[<prev] [next>] [<thread-prev] [day] [month] [year] [list]
Message-ID: <20171208004203.GF15263@port70.net>
Date: Fri, 8 Dec 2017 01:42:03 +0100
From: Szabolcs Nagy <nsz@...t70.net>
To: musl@...ts.openwall.com
Subject: Re: remquo - underlying logic
* Damian McGuckin <damianm@....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.
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
