gmp_printf bug?
Vincent Lefevre
vincent at vinc17.net
Mon Aug 8 03:46:32 CEST 2011
On 2011-08-04 22:17:45 +0200, Niels Möller wrote:
> The idea is, to compute log(b)/log(2), first compute sqrt(2),
> sqrt(sqrt(2)), .... Then compute the logrithm starting from the most
> significant bit by multiplying selected roots to get a value just below
> b.
>
> (And then the integer part of the logarithm must be done upfront, or,
> equivalently, you need to include repeated squares as well).
Wouldn't it be better to use a shift-and-add algorithm, so that only
additions are done instead of multiplications? I'm not sure about
the number of iterations, but it seems to be similar in both cases.
I've attached a small program I've written that shows how this can
be done. It uses double-precision arithmetic, but of course, the
same algorithm works in fixed point (so that mpz can be used). With
the included parameters, it can compute log2(x) with x between 1 and
2^64 at least.
--
Vincent Lefèvre <vincent at vinc17.net> - Web: <http://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / Arénaire project (LIP, ENS-Lyon)
-------------- next part --------------
/* Shift-and-add log2 computation inspired by Henry Briggs' algorithm.
No error analysis. Can also be done in fixed point. */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define K0 6
#define KN 53
static int s[K0+KN];
static double w[K0+KN];
void compsw (void)
{
int k;
for (k = K0; k > -KN; k--)
{
s[K0-k] = k <= 0 ? k : 1 << (k - 1);
w[K0-k] = log1p (ldexp (1.0, s[K0-k])) / log (2.0);
}
}
double my_log (double x)
{
double xk = 1.0, yk = 0.0, xt;
int k;
for (k = K0; k > -KN; k--)
{
xt = xk + ldexp (xk, s[K0-k]);
if (xt <= x)
{
xk = xt;
yk += w[K0-k];
}
}
return yk;
}
int main (int argc, char *argv[])
{
double x, y, z;
if (argc != 2)
{
fprintf (stderr, "Usage: briggs-log <double>\n");
exit (1);
}
compsw ();
x = atof (argv[1]);
y = my_log (x);
z = log2 (x);
printf ("x = %.20g\n", x);
printf ("[mine] log2(x) = %.20g\n", y);
printf ("[libc] log2(x) = %.20g\n", z);
return 0;
}
/* $Id: briggs-log.c 45502 2011-08-08 01:33:37Z vinc17/xvii $ */
More information about the gmp-devel
mailing list