* math argument reduction
@ 2011-06-26 18:04 Szabolcs Nagy
2011-06-26 18:33 ` Rich Felker
0 siblings, 1 reply; 3+ messages in thread
From: Szabolcs Nagy @ 2011-06-26 18:04 UTC (permalink / raw)
To: musl
dalias, i thought about the simple argument reduction method
you proposed on irc, imho it won't be that good
the task is to calculate X mod M, so to get the fractional part of
x = X/M
you said you would treat X as a sum (of 2^n terms) let's say
X = A+B
and assume we have precalculated a=A/M and b=B/M in a table, then
x = a+b
is easy to get
but it well might be that the form of a+b is
[..]1011 . 0000000000000001101011[..]
so the fractional part has many leading zeros
(it can be as many as 61) so you need to store a and b
with large precision to get good floating point representation
of the fractional part
(53 + 61 bits precision + some extra bits for rounding)
so you need to store and do arithmetics with 120bit fixpoint values
at that point it's probably not much harder to implement
the standard method of argument reduction used in sun
and bsd libm:
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.67.5616
^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: math argument reduction
2011-06-26 18:04 math argument reduction Szabolcs Nagy
@ 2011-06-26 18:33 ` Rich Felker
2011-06-26 18:51 ` Szabolcs Nagy
0 siblings, 1 reply; 3+ messages in thread
From: Rich Felker @ 2011-06-26 18:33 UTC (permalink / raw)
To: musl
On Sun, Jun 26, 2011 at 08:04:02PM +0200, Szabolcs Nagy wrote:
> dalias, i thought about the simple argument reduction method
> you proposed on irc, imho it won't be that good
>
> the task is to calculate X mod M, so to get the fractional part of
> x = X/M
> you said you would treat X as a sum (of 2^n terms) let's say
> X = A+B
> and assume we have precalculated a=A/M and b=B/M in a table, then
> x = a+b
> is easy to get
>
> but it well might be that the form of a+b is
>
> [..]1011 . 0000000000000001101011[..]
>
> so the fractional part has many leading zeros
> (it can be as many as 61) so you need to store a and b
> with large precision to get good floating point representation
> of the fractional part
> (53 + 61 bits precision + some extra bits for rounding)
>
> so you need to store and do arithmetics with 120bit fixpoint values
I think you misunderstood my algorithm. You never need to compute X/M.
It's A%M and B%M, not A/M and B/M, which you have precalculated in
your tables. And therefore, (A+B)%M is congruent to A%M + B%M. In
general it could be outside the interval [0,M), but not by so much
that it's hard to get a good answer.
Rich
^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: math argument reduction
2011-06-26 18:33 ` Rich Felker
@ 2011-06-26 18:51 ` Szabolcs Nagy
0 siblings, 0 replies; 3+ messages in thread
From: Szabolcs Nagy @ 2011-06-26 18:51 UTC (permalink / raw)
To: musl
* Rich Felker <dalias@aerifal.cx> [2011-06-26 14:33:44 -0400]:
> It's A%M and B%M, not A/M and B/M, which you have precalculated in
> your tables. And therefore, (A+B)%M is congruent to A%M + B%M. In
> general it could be outside the interval [0,M), but not by so much
> that it's hard to get a good answer.
the hard part is when there is a lot of leading zeros in the result
A%M + B%M can be just above M by a little so (A%M+B%M)%M is almost 0
eg: 0.00000000000000011111
to represent this number properly as float you need all those 1s eg
1.1111 * 2^-n
if A%M and B%M is not precise enough then those trailing 1s will get
lost so the abs value is about right but as a floating point value
the result does not have 53 correct bits..
^ permalink raw reply [flat|nested] 3+ messages in thread
end of thread, other threads:[~2011-06-26 18:51 UTC | newest]
Thread overview: 3+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2011-06-26 18:04 math argument reduction Szabolcs Nagy
2011-06-26 18:33 ` Rich Felker
2011-06-26 18:51 ` Szabolcs Nagy
Code repositories for project(s) associated with this public inbox
https://git.vuxu.org/mirror/musl/
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for NNTP newsgroup(s).