Mailing-List: contact dev-help@commons.apache.org; run by ezmlm
Precedence: bulk
Reply-To: "Commons Developers List" <dev@commons.apache.org>
Received-SPF: pass (nike.apache.org: local policy)
Date: Sun, 9 Oct 2011 18:38:04 +0200
From: Gilles Sadowski <gilles@harfang.homelinux.org>
To: dev@commons.apache.org
Subject: Re: svn commit: r1179928 - in [...]
Message-ID: <20111009163804.GJ17021@dusk.harfang.homelinux.org>
Mail-Followup-To: dev@commons.apache.org
References: <20111007131745.GN17388@dusk.harfang.homelinux.org>
 <4E8FDC32.3000407@gmail.com>
 <20111008120547.GF17021@dusk.harfang.homelinux.org>
 <4E905C63.20306@gmail.com>
 <20111008211224.GG17021@dusk.harfang.homelinux.org>
 <4E90C0B1.6050109@gmail.com>
 <20111008230628.GH17021@dusk.harfang.homelinux.org>
 <CAGRH7HoWf78B0iLOCUgMF-Oq_u=q3SQYPnzrwVdRyX+iHqh86A@mail.gmail.com>
 <20111009120907.GI17021@dusk.harfang.homelinux.org>
 <4E91B39F.2030404@gmail.com>
MIME-Version: 1.0
Content-Type: text/plain; charset=iso-8859-1
Content-Disposition: inline
Content-Transfer-Encoding: 8bit
In-Reply-To: <4E91B39F.2030404@gmail.com>
User-Agent: Mutt/1.5.21 (2010-09-15)
X-DCC-scarlet.be-Metrics: eir 20001; Body=1 Fuz1=1 Fuz2=1
X-Virus-Scanned: clamav-milter 0.97.1-exp at eir
X-Virus-Status: Clean
X-Virus-Checked: Checked by ClamAV on apache.org

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.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?
> 
> 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 (-1E-9, 1E-9), 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, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org