On Sun, Oct 09, 2011 at 07:45:51AM 0700, Phil Steitz wrote:
> On 10/9/11 5:09 AM, Gilles Sadowski wrote:
> > On Sun, Oct 09, 2011 at 10:32:38AM +0200, Sébastien Brisard wrote:
> >>> Answering with a few examples:
> >>> x=1.000000000000000000e15 f=1.000000000000000000e+00 s=1.000000000000000000e+00
> >>> x=5.000000000000001000e15 f=1.000000000000000000e+00 s=1.000000000000000000e+00
> >>> x=2.500000000000000400e14 f=1.000000000000000000e+00 s=1.000000000000000000e+00
> >>> x=1.250000000000000200e13 f=1.000000000000000000e+00 s=1.000000000000000000e+00
> >>> x=6.250000000000001000e13 f=1.000000000000000000e+00 s=1.000000000000000000e+00
> >>> x=3.125000000000000500e12 f=1.000000000000000000e+00 s=1.000000000000000000e+00
> >>> x=1.562500000000000400e11 f=1.000000000000000000e+00 s=1.000000000000000000e+00
> >>> x=7.812500000000003000e11 f=1.000000000000000000e+00 s=1.000000000000000000e+00
> >>> x=3.906250000000001400e10 f=1.000000000000000000e+00 s=1.000000000000000000e+00
> >>> x=1.953125000000000700e09 f=1.000000000000000000e+00 s=1.000000000000000000e+00
> >>> x=9.765625000000004000e09 f=1.000000000000000000e+00 s=1.000000000000000000e+00
> >>> x=4.882812500000002000e08 f=9.999999999999996000e01 s=9.999999999999996000e01
> >>> x=2.441406250000001000e07 f=9.999999999999900000e01 s=9.999999999999900000e01
> >>> x=1.220703125000000700e06 f=9.999999999997516000e01 s=9.999999999997516000e01
> >>> x=6.103515625000004000e06 f=9.999999999937912000e01 s=9.999999999937912000e01
> >>> x=3.051757812500002000e05 f=9.999999998447796000e01 s=9.999999998447796000e01
> >>> x=1.525878906250001000e04 f=9.999999961194893000e01 s=9.999999961194893000e01
> >>> x=7.629394531250005000e04 f=9.999999029872346000e01 s=9.999999029872346000e01
> >>> x=3.814697265625002600e03 f=9.999975746825600000e01 s=9.999975746825600000e01
> >>> x=1.907348632812501400e02 f=9.999393681227797000e01 s=9.999393681227797000e01
> >>>
> >>> Thus: below 9.765625000000004E9, 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 "1e9").]
> >>>
> >>> 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 1E9 value is actually a mathematically provable value,
> >> since sin(x)/x is an alternating series, so the difference
> >> sin(x)/x1 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 < 1E16, which is actually less
> >> strong than x < 1E9.
> >> 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 floatingpoint
> >> (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?
>
> No, please at least point out that we define the value at 0 to be
> 1. There is no reason to link to Wikipedia,
This gives more than just the definition (e.g. inspiration for unit tests).
> when you can directly
> provide the formula, including treatment of 0.
Done in revision 1180642.
> The reason that I wanted to point out the top coding was that I
> thought it might be possible results that were equal, but had
> different internal representations could be returned over the
> interval (1E9, 1E9), where the actual value is close to but less
> than 1. This could effect computations involving the returned
> values. I have not been able to demonstrate this (and suspect that
> I may in fact be misreading the JLS on this issue  i.e., on method
> return you have to map to the standard value set, so returned values
> have to have the same internal representation), so I am fine leaving
> out the reference to the cut point.
Do you mean the numbers representation used inside the CPU? I think that
this is beyond the scope of a Java library.
As for the numbers returned from the above table they are _strictly_ the
same (see the unit test, where the tolerance is zero).
Gilles

To unsubscribe, email: devunsubscribe@commons.apache.org
For additional commands, email: devhelp@commons.apache.org
