On Tue, 20 Apr 2021 at 13:51, Gilles Sadowski <gillese...@gmail.com> wrote:
> Hi Alex.
>
> Le mar. 20 avr. 2021 à 11:17, Alex Herbert <alex.d.herb...@gmail.com> a
> écrit :
> >
> > [...]
>
> I'm a bit lost in all these bits... ;-)
>
I also had to remind myself of the work since it was a long time ago.
>
> > Any opinions on how changing LinearComination may affect Geometry? Either
> > we clean up the implementation to use the fast dot2s algorithm with
> correct
> > support for splitting 53-bit mantissas, or switch to the extended
> precision
> > version. But I do not think we should leave it as the current
> > implementation which has disadvantages against either of the
> alternatives.
>
> What is your suggestion?
>
Some background...
The dot2s is a specialisation of the DotK algorithm where the K represents
the k-fold increase in precision over a standard scalar product and denotes
qualitatively the level of robustness for a sum with massive cancellation
(e.g. 1e-100 + 1e200 - 2e-100 - 1e200 = -1e-100). Any floating-point sum is
limited by the representation of a 64-bit double as a binary floating point
number with only 53-bits of precision and an exponent. When adding numbers
the result is computed and any extra bits that cannot fit into the mantissa
are lost. Whether this occurs is dependent on the input numbers, but in
general as the difference in magnitude becomes greater then the result will
contain progressively less of the true result due to the 53-bit mantissa.
If two numbers are different by more than 2^53 then adding them makes no
difference to the bigger number. In the example above a standard precision
result is 0. Only with extended precision can you hold the result of 1e-100
+ 1e200, etc.
In a sum of magnitude this only really matters when cancellation (opposite
signs) may occur. Take the recently discussed code to compute the 3D vector
length:
len = sqrt(x^2 + y^2 + z^2)
This uses a sum where the sign of the terms is the same. There will be no
cancellation and thus the answer will be close to the true result (within a
few ulps).
But if the sum of multiple values has different signs then the answer may
be in between the magnitudes of the terms. So the intermediate sum should
store knowledge of more bits of the intermediate result than the 64-bits of
a double.
The LinearCombination is actually providing two operations, multiplication
and addition:
sum_i {a_i * b_i}
The multiplication has a potential 106-bit mantissa for the result with a
single exponent. Any addition may require a longer mantissa with a single
exponent which is how BigDecimal would store the result. An extended
precision double number would represent the result as an array of two
doubles of very different magnitudes. Adding or multiplying the extended
precision result can increase the length of the result further. The use of
arrays of double to represent a number in extended precision is well
described in Shewchuk (1997): Arbitrary Precision Floating-Point Arithmetic
[1], on which a lot of the work in Numbers-142 is based.
There could be a case for a class to implement this generically:
public class ExtendedDouble {
public static ExtendedDouble of(double a);
public ExtendedDouble add(double a);
public ExtendedDouble add(ExtendedDouble a);
public ExtendedDouble multiply(double a);
public ExtendedDouble multiply(ExtendedDouble a);
public double toDouble();
}
That could be in another module in numbers. Shewchuk does not describe
division but it may be in references he provides. His paper is based upon
developing exact arithmetic for geometry applications so addition and
multiplication are all that are used.
As for LinearCombination then I would suggest that the use case is for any
dot product where cancellation may be expected. If handling cancellation is
more important than speed then using the most accurate method would seem to
be a better choice.
I do recall seeing LinearCombination used in the Geometry code inside a
loop with the result stored as a double on each iteration. This would
result in loss of the extended precision intermediate sum and so could
cause error due to cancellation. So there is a case for going through
Geometry to find use cases for LinearCombination and potentially find
places where changes should be made to compute all the terms in advance and
call LinearCombination once with the two input arrays.
> My impression is that [Geometry] emphasizes accuracy over ultimate
> speed (as opposed to libraries used for real-time rendering, I guess).
>
> However, could it be possible to leave this as a user's decision?
> Quoting from Matt's tutorial:
> ---CUT---
> Typically, users of Commons Geometry will construct a single instance
> of this type for use by multiple objects throughout an entire
> operation, or even application. Since we don't want our class to
> assume such a heavy responsibility, we will simply accept a
> DoublePrecisionContext in the constructor.
> ---CUT---
> Would it be conceivable that the choice of the implementation
> activated by a call to the "LinearCombination" facility is also
> encapsulated in the "DoublePrecisionContext"?
>
Numbers-142 did discuss the change of LinearCombination to an interface.
Then we already have several implementations for the interface where a user
can choose the balance between speed and robustness to cancellation.
It would mean reworking Geometry at all places that use LinearCombination
to use an instance rather than the current static method call. There would
be a default instance or a user supplied one.
I would not opt for using a singleton to define the behaviour of the static
method. This would work for a standalone application using Geometry but
would have issues for any application where different code requires a
different precision level.
[1]
http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
