On Sun, Oct 09, 2011 at 10:32:38AM +0200, Sébastien Brisard wrote: > > > > Answering with a few examples: > > x=1.000000000000000000e-15 f=1.000000000000000000e+00 > > s=1.000000000000000000e+00 > > x=5.000000000000001000e-15 f=1.000000000000000000e+00 > > s=1.000000000000000000e+00 > > x=2.500000000000000400e-14 f=1.000000000000000000e+00 > > s=1.000000000000000000e+00 > > x=1.250000000000000200e-13 f=1.000000000000000000e+00 > > s=1.000000000000000000e+00 > > x=6.250000000000001000e-13 f=1.000000000000000000e+00 > > s=1.000000000000000000e+00 > > x=3.125000000000000500e-12 f=1.000000000000000000e+00 > > s=1.000000000000000000e+00 > > x=1.562500000000000400e-11 f=1.000000000000000000e+00 > > s=1.000000000000000000e+00 > > x=7.812500000000003000e-11 f=1.000000000000000000e+00 > > s=1.000000000000000000e+00 > > x=3.906250000000001400e-10 f=1.000000000000000000e+00 > > s=1.000000000000000000e+00 > > x=1.953125000000000700e-09 f=1.000000000000000000e+00 > > s=1.000000000000000000e+00 > > x=9.765625000000004000e-09 f=1.000000000000000000e+00 > > s=1.000000000000000000e+00 > > x=4.882812500000002000e-08 f=9.999999999999996000e-01 > > s=9.999999999999996000e-01 > > x=2.441406250000001000e-07 f=9.999999999999900000e-01 > > s=9.999999999999900000e-01 > > x=1.220703125000000700e-06 f=9.999999999997516000e-01 > > s=9.999999999997516000e-01 > > x=6.103515625000004000e-06 f=9.999999999937912000e-01 > > s=9.999999999937912000e-01 > > x=3.051757812500002000e-05 f=9.999999998447796000e-01 > > s=9.999999998447796000e-01 > > x=1.525878906250001000e-04 f=9.999999961194893000e-01 > > s=9.999999961194893000e-01 > > x=7.629394531250005000e-04 f=9.999999029872346000e-01 > > s=9.999999029872346000e-01 > > x=3.814697265625002600e-03 f=9.999975746825600000e-01 > > s=9.999975746825600000e-01 > > x=1.907348632812501400e-02 f=9.999393681227797000e-01 > > s=9.999393681227797000e-01 > > > > Thus: below 9.765625000000004E-9, the value of the definitional formula > > (without shortcut, indicated by "f=" in the above table) is strictly the > > same as the CM implementation (with shortcut, indicated by "s=" in the above > > table) in that they are both the double value "1.0". > > > > [I still don't understand what you mean by "(despite all being equal to 1 > > under double equals)".] > > > > What the implementation returns is, within double precision, the result of > > the mathematical definition of sinc: sin(x) / x. Hence, there is no *need* > > to document any special case, not even x = 0: The fact that without the > > shortcut check, we'd get NaN does not mean that the implementation departs > > from the definition, but rather that it correctly simulates it (at double > > precision). > > [However, if we assume that the doc should integrate numerical > > considerations, it doesn't hurt to add the extra prose (in which case it > > becomes necessary, IMHO, to explain the "1e-9").] > > > > Maybe I should also add a "SincTest" class where a unit test would check the > > strict equality of returned values from the actual division and from the > > shortcut implementation? > > > > > > Gilles > > > I think the 1E-9 value is actually a mathematically provable value, > since sin(x)/x is an alternating series, so the difference > |sin(x)/x-1| is bounded from above by the next term in the taylor > expansion, namely x^2/6. Replacing sinc(x) by one is therefore > rigorous provided this error remains below one ulp of one. This leads > to the following condition x^2/6 < 1E-16, which is actually less > strong than |x| < 1E-9. > So I personally think that this is indeed an implementation detail. If > you look at the implementation of any basic functions, there are > *always* boundary cases, with different formulas for different ranges > of the argument. These cases are not advertised anywhere in the doc > (and we should be thankful of that, in my opinion). > As a user of sinc (and spherical Bessel functions in general, for > diffraction problems), the only thing I really care about is: > reliability near zero. How this reliability is enforced is really none > of my concerns. > One further point. If you were to try and implement the next spherical > Bessel function, you would find that the analytical floating-point > (without short cut, using Gilles' terminology) estimate near zero does > *not* behave well, because, I think, of catastrophic cancellations. So > in this case, you *have* to carry out a test on x, and if x is too > small, you have to replace the analytical expression of j1 by its > Taylor expansion, in order to remain within one ulp of the expected > value. The boundary is found using the preceding line of reasoning > (alternating series). In this case, this is still an implementation > detail, but *not* an optimization one-- it is a *requirement*!
I made some changes (revision 1180588), following the above argument which concurs with my original remark in this thread (implementation detail). Is this now satisfactory? Gilles --------------------------------------------------------------------- To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org For additional commands, e-mail: dev-h...@commons.apache.org
