miguelraz / OrthogonalPolynomials.jl

Other
6 stars 1 forks source link

Accuracy error in the readme? #16

Open MikaelSlevinsky opened 6 years ago

MikaelSlevinsky commented 6 years ago

In the context of numerical evaluation of L_{12}(1/2), the readme states

Note as well that because we are reducing dramatically the number of floating point operations, we are losing less precision than other methods.

Is this in reference to the method named goal12, whose result appears to be -0.23164963886602852?

The evaluation of L_{12}(1/2) can be done in Rational{BigInt} that avoids floating-point rounding errors. The result is:

julia> L12(x) = (1//479001600)*(x^12 - 144x^11 + 8_712x^10 - 290_400x^9 + 5_880_600x^8 - 75_271_680x^7 + 614_718_720x^6 - 3_161_410_560x^5 + 9_879_408_000x^4 - 17_563_392_000x^3 + 15_807_052_800x^2 - 5_748_019_200x + 479_001_600)
L12 (generic function with 1 method)

julia> L12(big(1//2))
-454494403199//1961990553600

julia> BigFloat(ans)
-2.316496388655191535372741969964192186414408636630858853081075303297525519747738e-01

julia> Float64(ans)
-0.23164963886551915

Perhaps someone may clarify to me how the goal12 method is more accurate than the other methods, notably the ApproxFun implementation that results in -0.2316496388655192, which is two units of least precision away from the correctly rounded rational result.

miguelraz commented 6 years ago

@MikaelSlevinsky ! Welcome and thanks for bringing this up. I don't think that I can clarify this, because assuming that less floating point operations leading to a more precise result is wrong in the general case, as you have suitably pointed it out. I will strike that from the README.md. Thanks!

MikaelSlevinsky commented 6 years ago

P.S. if you want more info about stable recurrences for orthogonal polynomials, check out Section 4.2 here, and here.