Re: remquo - underlying logic

[Szabolcs Nagy <nsz@port70.net>]

* Damian McGuckin <damianm@esi.com.au> [2017-12-07 12:09:23 +1100]:
> That said, over that range, I am experimenting using a simplistic form of
> double-double arithmetic for that calculation.
> 
> Would you agree that when
> 
> 	(int) (x / y) < 2^52
> 
> the computation (int) (x / y) is accurate to within epsilon, i.e. if it
> should be at most be incorrect by +/- 1.?
> 

this part is not the issue:

if the exact mathematical (int)(x/y) is k, i.e.

k <= x/y < k+1

then the rounded double prec x/y is

k <= x/y <= k+1

because close to k+1 some values get rounded up,
so sometimes you would compute x-(k+1)*y instead
of x-k*y, but this is easy to correct: just add +y
in the end if the result is out of range.
(same is true when other round-to-int method is
used instead of trunc, assumes k and k+1 are
representable)

> If so, and using the same sort of logic that log.c uses to split the
> calculation of
> 
> 	k * log(2.0)
> 
> into a high and low component, or maybe into 4 components, would you agree
> that it is possible to come up with an accurate computation of
> 
> 	x - y * (int) (x / y)
> 
> It should be much quicker than long division.

the real problem is doing x-k*y exactly (it is
representable as double).

when evaluated without fma, there is a rounding
error on the mul (the sub is exact).

one possibility is if k < 2^26 and the bottom 26
bits of y are zeroed out in yz then

x - k*yz - k*(y-yz)

is exact, another way is to use veltkamp/dekker
exact multiplication algorithm.