gmpbench-0.1 and gexpr.c and Win32

[Sisyphus]

Hi,

gexpr.c now compiles on Win32 with the latest MinGW (gcc). Earlier 
versions of MinGW's math.h did not include asinh, acosh, or atanh 
functions, and gexpr.c needed modification before it would compile.

Attached is such a modified gexpr.c and fastmath.h (which is needed by 
the modifed gexpr.c). It still won't compile with MSVC 7.0, owing to 
issues with syntax (C99 standards and assembler code) in fastmath.h - 
which I don't know how to fix. The asinh(), acosh() and atanh() code was 
pinched directly from the MinGW files asinh.c, acosh.c and atanh.c, 
which can be viewed at 
http://cvs.sourceforge.net/viewcvs.py/mingw/runtime/mingwex/math/

'fastmath.h' was also pinched directly from the same location. I'm 
unsure of what's required here as regards copyright and license issues 
so I haven't provided a patch. (Chances are that if I *did* provide a 
patch, it would be inadequate for legal reasons :-)

First problem with runbench is that there's an assumption that gcc 
always produces an executable named 'a.out' by default - but on Win32 
the default is 'a.exe'. This, of course, poses a problem when, in 
runbench, we reach the line:

echo "Using `./a.out`"

I think the correct (portable) fix is to change the preceding line from:

$CC $CFLAGS gmpver.c $LIBS

to:

$CC $CFLAGS -o a.out gmpver.c $LIBS

At least that works fine with Win32.


multiply.c, divide.c and rsa.c each contain the following piece of code:

#if !defined (__sun) \
    && (defined (USG) || defined (__SVR4) || defined (_UNICOS) \
        || defined (__hpux))

For Win32, we need to tack on '|| defined (_WIN32)' so that it becomes:

#if !defined (__sun) \
    && (defined (USG) || defined (__SVR4) || defined (_UNICOS) \
        || defined (__hpux) || defined (_WIN32))

With those changes in place, and gexpr.exe in the path, './runbench' 
runs fine in the MSYS shell.

The only problem I have is that if I use the GMP library built using 
default configuration values, I get a segfault-type error when it comes 
to divide 8388608 4194304. If I use the GMP library built with malloc 
(instead of alloca) then that problem does not arise. I assume the error 
arises because I have only 64Mb of ram - which is a ridiculously small 
amount these days. I believe that issue can therefore be ignored by the 
GMP developers.

Cheers,
Rob





 /***********************************************************
 fastmath.h ************************************************/

 static __inline__ double __fast_sqrt (double x)
{
  double res;
  asm __volatile__ ("fsqrt" : "=t" (res) : "0" (x));
  return res;
}

static __inline__ long double __fast_sqrtl (long double x)
{
  long double res;
  asm __volatile__ ("fsqrt" : "=t" (res) : "0" (x));
  return res;
}

static __inline__ float __fast_sqrtf (float x)
{
  float res;
  asm __volatile__ ("fsqrt" : "=t" (res) : "0" (x));
  return res;
}


static __inline__ double __fast_log (double x)
{
   double res;
   asm __volatile__
     ("fldln2\n\t"
      "fxch\n\t"
      "fyl2x"
       : "=t" (res) : "0" (x) : "st(1)");
   return res;
}

static __inline__ long double __fast_logl (long double x)
{
  long double res;
   asm __volatile__
     ("fldln2\n\t"
      "fxch\n\t"
      "fyl2x"
       : "=t" (res) : "0" (x) : "st(1)");
   return res;
}


static __inline__ float __fast_logf (float x)
{
   float res;
   asm __volatile__
     ("fldln2\n\t"
      "fxch\n\t"
      "fyl2x"
       : "=t" (res) : "0" (x) : "st(1)");
   return res;
}

static __inline__ double __fast_log1p (double x)
{
  double res;
  /* fyl2xp1 accurate only for |x| <= 1.0 - 0.5 * sqrt (2.0) */
  if (fabs (x) >= 1.0 - 0.5 * 1.41421356237309504880)
    res = __fast_log (1.0 + x);
  else
    asm __volatile__
      ("fldln2\n\t"
       "fxch\n\t"
       "fyl2xp1"
       : "=t" (res) : "0" (x) : "st(1)");
   return res;
}

static __inline__ long double __fast_log1pl (long double x)
{
  long double res;
  /* fyl2xp1 accurate only for |x| <= 1.0 - 0.5 * sqrt (2.0) */
  if (fabsl (x) >= 1.0L - 0.5L * 1.41421356237309504880L)
    res = __fast_logl (1.0L + x);
  else
    asm __volatile__
      ("fldln2\n\t"
       "fxch\n\t"
       "fyl2xp1"
       : "=t" (res) : "0" (x) : "st(1)");
   return res;
}

static __inline__ float __fast_log1pf (float x)
{
  float res;
  /* fyl2xp1 accurate only for |x| <= 1.0 - 0.5 * sqrt (2.0) */
  if (fabsf (x) >= 1.0 - 0.5 * 1.41421356237309504880)
    res = __fast_logf (1.0 + x);
  else
    asm __volatile__
      ("fldln2\n\t"
       "fxch\n\t"
       "fyl2xp1"
       : "=t" (res) : "0" (x) : "st(1)");
   return res;
}

/* Expression evaluation using plain floating-point arithmetic.
   Copyright (C) 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc.

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with
this program; see the file COPYING.  If not, write to the Free Software
Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */

#include <string.h>
#include <stdio.h>
#include <setjmp.h>
#include <math.h>

#ifdef _WIN32
#include <errno.h>
#include "fastmath.h"
#endif

jmp_buf top;

static char *bp, *expr_start_p, *expr_end_p;
static int ch;
static int previous_result_valid_flag;
double previous_result;

double floor (), strtod ();
double term (), expo (), factor (), number ();

#define next skip ()

void
skip ()
{
  do
    ch = *bp++;
  while (ch == ' ' || ch == '\n');
}

void
error ()
{
  fprintf (stderr, " ^ syntax error\n");
  longjmp (top, 1);
}

#ifdef _WIN32

  /* asinh(x) = copysign(log(fabs(x) + sqrt(x * x + 1.0)), x) */
 double asinh(double x)
 {
   double z;
   if (!isfinite (x))
     return x;
   z = fabs (x);
 
   /* Avoid setting FPU underflow exception flag in x * x. */
 #if 0
   if ( z < 0x1p-32)
     return x;
 #endif
 
   /* Use log1p to avoid cancellation with small x. Put
      x * x in denom, so overflow is harmless. 
      asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0)
               = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0))  */
 
   z = __fast_log1p (z + z * z / (__fast_sqrt (z * z + 1.0) + 1.0));
 
   return ( x > 0.0 ? z : -z);
 }

 /* atanh (x) = 0.5 * log ((1.0 + x)/(1.0 - x)) */
 
 double atanh(double x)
 {
   double z;
   if isnan (x)
     return x;
   z = fabs (x);
   if (z == 1.0)
     {
       errno  = ERANGE;
       return (x > 0 ? INFINITY : -INFINITY);
     }
   if (z > 1.0)
     {
       errno = EDOM;
       return nan("");
     }
   /* Rearrange formula to avoid precision loss for small x.
 
   atanh(x) = 0.5 * log ((1.0 + x)/(1.0 - x))
           = 0.5 * log1p ((1.0 + x)/(1.0 - x) - 1.0)
            = 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x)) 
            = 0.5 * log1p ((2.0 * x ) / (1.0 - x))  */
   z = 0.5 * __fast_log1p ((z + z) / (1.0 - z));
   return x >= 0 ? z : -z;
 }
 
 /* acosh(x) = log (x + sqrt(x * x - 1)) */
 double acosh (double x)
 {
   if (isnan (x)) 
     return x;
 
   if (x < 1.0)
     {
       errno = EDOM;
       return nan("");
     }
 
   if (x > 0x1p32)
     /*  Avoid overflow (and unnecessary calculation when
         sqrt (x * x - 1) == x). GCC optimizes by replacing
         the long double M_LN2 const with a fldln2 insn.  */ 
     return __fast_log (x) + 6.9314718055994530941723E-1L;
 
   /* Since  x >= 1, the arg to log will always be greater than
      the fyl2xp1 limit (approx 0.29) so just use logl. */ 
   return __fast_log (x + __fast_sqrt((x + 1.0) * (x - 1.0)));
 }

#endif

/* END OF WIN32-SPECIFIC STUFF */

double
expr ()
{
  double e;
  if (ch == '+')
    {
      next;
      e = term ();
    }
  else if (ch == '-')
    {
      next;
      e = - term ();
    }
  else
    e = term ();
  while (ch == '+' || ch == '-')
    { char op;
      op = ch;
      next;
      if (op == '-')
        e -= term ();
      else
        e += term ();
    }
  return e;
}

double
term ()
{
  double t;
  t = expo ();
  for (;;)
    switch (ch)
      {
      case '*':
        next;
        t *= expo ();
        break;
      case '/':
        next;
        t /= expo ();
        break;
      case '%':
        next;
        t = fmod (t, expo ());
        break;
      default:
        return t;
      }
}

double
expo ()
{
  double e;
  e = factor ();
  if (ch == '^')
    {
      next;
      e = pow (e, expo ());
    }
  return e;
}

struct functions
{
  char *spelling;
  double (* evalfn) ();
};

double
log2 (x)
     double x;
{
  return log (x) * 1.4426950408889634073599246810019;
}

struct functions fns[] =
{
  {"log2", log2},
  {"log10", log10},
  {"log", log},
  {"exp", exp},
  {"sqrt", sqrt},
  {"floor", floor},
  {"ceil", ceil},
  {"sin", sin},
  {"cos", cos},
  {"tan", tan},
  {"asin", asin},
  {"acos", acos},
  {"atan", atan},
  {"sinh", sinh},
  {"cosh", cosh},
  {"tanh", tanh},
  {"asinh", asinh},
  {"acosh", acosh},
  {"atanh", atanh},
  {0, 0}
};


double
factor ()
{
  double f;
  int i;

  for (i = 0; fns[i].spelling != 0; i++)
    {
      char *spelling = fns[i].spelling;
      int len = strlen (spelling);
      if (strncmp (spelling, bp - 1, len) == 0 && ! isalnum (bp[-1 + len]))
        {
          bp += len - 1;
          next;
          if (ch != '(')
            error ();
          next;
          f = expr ();
          if (ch != ')')
            error ();
          next;
          return (fns[i].evalfn) (f);
        }
    }

  if (ch == '(')
    {
      next;
      f = expr ();
      if (ch == ')')
        next;
      else
        error ();
    }
  else
    f = number ();
  if (ch == '!')
    {
      unsigned long n;
      if (floor (f) != f)
        error ();
      for (n = f, f = 1; n > 1; n--)
        f *= n;
      next;
    }
  return f;
}

double
number ()
{
  double n;
  char *endp;

  if (strncmp ("pi", bp - 1, 2) == 0 && ! isalnum (bp[1]))
    {
      bp += 2 - 1;
      next;
      return 3.1415926535897932384626433832795;
    }
  if (ch == '$')
    {
      if (! previous_result_valid_flag)
        error ();
      next;
      return previous_result;
    }
  if (ch != '.' && (ch < '0' || ch > '9'))
    error ();
  n = strtod (bp - 1, &endp);
  if (endp == bp - 1)
    error ();
  bp = endp;
  next;
  return n;
}

main (argc, argv)
     int argc;
     char **argv;
{
  int nl_flag = 1;
  int hhmm_flag = 0;
  int dhhmm_flag = 0;
  int round_flag = 0;
  int prec = 5;

  while (argc >= 2)
    {
      if (!strcmp (argv[1], "-n"))
        {
          nl_flag = 0;
          argv++;
          argc--;
        }
      else if (!strcmp (argv[1], "-hhmm"))
        {
          hhmm_flag = 1;
          argv++;
          argc--;
        }
      else if (!strcmp (argv[1], "-dhhmm"))
        {
          dhhmm_flag = 1;
          argv++;
          argc--;
        }
      else if (!strcmp (argv[1], "-round"))
        {
          round_flag = 1;
          argv++;
          argc--;
        }
      else if (!strcmp (argv[1], "-prec"))
        {
          prec = atoi (argv[2]);
          argv += 2;
          argc -= 2;
        }
      else if (!strcmp (argv[1], "-help") || !strcmp (argv[1], "-h"))
        {
          printf ("usage: %s [options] expr [expr ... ]\n", argv[0]);
          printf ("   options:    -n      -- suppress newline\n");
          printf ("               -prec n -- print n digits\n");
          printf ("               -round  -- round to nearesr integer\n");
          printf ("               -hhmm   -- print in base 60 (time format)\n");
          printf ("               -dhhmm   -- print in base 24,60,60 (time 
format)\n");
          printf ("               -help   -- make demons drop from your 
nose\n");
          exit (0);
        }
      else
        break;
    }

  if (argc >= 2)
    {
      int i;
      double exp;

      for (i = 1; i < argc; i++)
        {
          expr_start_p = argv[i];
          expr_end_p = expr_end_p + strlen (expr_start_p);
          bp = expr_start_p;
          next;
          if (setjmp (top) == 0)
            {
              int h, m;
              exp = expr ();
              if (round_flag)
                exp = floor (exp + 0.5);
              if (hhmm_flag)
                {
                  h = (int) floor (exp);
                  m = (int) ((exp - floor (exp)) * 60.0 + 0.5);
                  h += m / 60;  m = m % 60;
                  printf ("%02d:%02d", h, m);
                }
              else if (dhhmm_flag)
                {
                  double tmp = (exp - floor (exp)) * 24.0;
                  h = (int) floor (tmp);
                  m = (int) ((tmp - floor (tmp)) * 60.0 + 0.5);
                  h += m / 60;  m = m % 60;
                  printf ("%dd %02d:%02d", (int) floor (exp), h, m);
                }
              else
                printf ("%.*g", prec, exp);
              if (nl_flag)
                puts ("");
              previous_result = exp;
              previous_result_valid_flag = 1;
            }
        }
    }
  else
    {
#define BUFSIZE 1024
      char buf[BUFSIZE];
      double exp;

      for (;;)
        {
          fputs ("eval> ", stdout);
          bp = fgets (buf, BUFSIZE, stdin);
          if (bp == NULL)
            break;
          next;
          if (setjmp (top) == 0)
            {
              int h, m;
              exp = expr ();
              if (round_flag)
                exp = floor (exp + 0.5);
              if (hhmm_flag)
                {
                  h = (int) floor (exp);
                  m = (int) ((exp - floor (exp)) * 60.0 + 0.5);
                  h += m / 60;  m = m % 60;
                  printf ("%02d:%02d", h, m);
                }
              else if (dhhmm_flag)
                {
                  double tmp = (exp - floor (exp)) * 24.0;
                  h = (int) floor (tmp);
                  m = (int) ((tmp - floor (tmp)) * 60.0 + 0.5);
                  h += m / 60;  m = m % 60;
                  printf ("%dd %02d:%02d", (int) floor (exp), h, m);
                }
              else
                printf ("%.*g", prec, exp);
              if (nl_flag)
                puts ("");
              previous_result = exp;
              previous_result_valid_flag = 1;
            }
        }
    }

  exit (0);
}