From mboxrd@z Thu Jan  1 00:00:00 1970
X-Msuck: nntp://news.gmane.org/gmane.linux.lib.musl.general/12229
Path: news.gmane.org!.POSTED!not-for-mail
From: Szabolcs Nagy <nsz@port70.net>
Newsgroups: gmane.linux.lib.musl.general
Subject: Re: remquo - underlying logic
Date: Fri, 8 Dec 2017 01:42:03 +0100
Message-ID: <20171208004203.GF15263@port70.net>
References: <alpine.LRH.2.02.1712010059520.23842@key0.esi.com.au>
 <20171130185956.GS15263@port70.net>
 <alpine.LRH.2.02.1712010723070.4154@key0.esi.com.au>
 <20171130211713.GT15263@port70.net>
 <alpine.LRH.2.02.1712061141550.18308@key0.esi.com.au>
 <20171206103708.GA15263@port70.net>
 <alpine.LRH.2.02.1712071157010.3671@key0.esi.com.au>
Reply-To: musl@lists.openwall.com
NNTP-Posting-Host: blaine.gmane.org
Mime-Version: 1.0
Content-Type: text/plain; charset=us-ascii
X-Trace: blaine.gmane.org 1512693734 14513 195.159.176.226 (8 Dec 2017 00:42:14 GMT)
X-Complaints-To: usenet@blaine.gmane.org
NNTP-Posting-Date: Fri, 8 Dec 2017 00:42:14 +0000 (UTC)
User-Agent: Mutt/1.6.0 (2016-04-01)
To: musl@lists.openwall.com
Original-X-From: musl-return-12245-gllmg-musl=m.gmane.org@lists.openwall.com Fri Dec 08 01:42:10 2017
Return-path: <musl-return-12245-gllmg-musl=m.gmane.org@lists.openwall.com>
Envelope-to: gllmg-musl@m.gmane.org
Original-Received: from mother.openwall.net ([195.42.179.200])
	by blaine.gmane.org with smtp (Exim 4.84_2)
	(envelope-from <musl-return-12245-gllmg-musl=m.gmane.org@lists.openwall.com>)
	id 1eN6jy-0003e1-Dg
	for gllmg-musl@m.gmane.org; Fri, 08 Dec 2017 01:42:10 +0100
Original-Received: (qmail 32307 invoked by uid 550); 8 Dec 2017 00:42:15 -0000
Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm
Precedence: bulk
List-Post: <mailto:musl@lists.openwall.com>
List-Help: <mailto:musl-help@lists.openwall.com>
List-Unsubscribe: <mailto:musl-unsubscribe@lists.openwall.com>
List-Subscribe: <mailto:musl-subscribe@lists.openwall.com>
List-ID: <musl.lists.openwall.com>
Original-Received: (qmail 32283 invoked from network); 8 Dec 2017 00:42:14 -0000
Mail-Followup-To: musl@lists.openwall.com
Content-Disposition: inline
In-Reply-To: <alpine.LRH.2.02.1712071157010.3671@key0.esi.com.au>
Xref: news.gmane.org gmane.linux.lib.musl.general:12229
Archived-At: <http://permalink.gmane.org/gmane.linux.lib.musl.general/12229>

* 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.