gmp_printf bug?

[Vincent Lefevre]

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 $ */