A faster pow(x,1/12) function available in Fast12thRoot.jl

@MrBelette ,

Further analysis suggests my earlier proposal is seriously flawed.

In part due to issues that Jeffrey raised about instability with certain sets of coefficients, and also because I didn’t pay attention to being clear about what was required.

  • Is the desired interval [\frac{1}{2},1] or [\frac{1}{2},1) ?
    The semi-open interval does not need the special case of t=0, so its perhaps not sensible to use transform t=x-1 for the best approximation.
  • The interval [\frac{1}{2},1) is best for simple range/strength reduction, but is it best for simple polynomial approximation ? (I conjecture [\frac{3}{4},\frac{3}{2}) might be numerically more stable, but will be more expensive for range reduction.)
  • What is the actual required precision ?
    I was targeting faithful rounding (which is basically \pm1 ulp) but you seem to be looking for correct rounding (100% \pm\frac{1}{2} ulp). However, without resolving the Table Makers Dilemma for this function (which I am not volunteering to do) I don’t know how much intermediate precision is required to guarantee correctness.

But If I go back and look at the text of your question, your use case suggests you code will end up a line that looks like fast12throot(x)-1, (see for example converting annual interest rate to monthly when compounding frequency known) so I wonder whether what you really want is the function fast12throotm1(x)\equivfast12throot(x)-1. It is often undesirable to calculate (x+1)-1 for small x because it loses precision on the intermediate step, which is exactly why standard maths libraries include functions for log1p(x)\equiv\log(x)+1 and expm1(x)\equiv\exp(x)-1 specifically to avoid problems with precision near special cases (since the implementations include approximations then add 1 or subtract 1 respectively).

So I suppose there are two questions I should have asked earlier:

  1. Do you really want both functions ?
  2. What is actually the acceptable error for the desired function(s) ?

regards,
JonT.

2 Likes