[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
gmpbench-0.1 and gexpr.c and Win32
From: |
Sisyphus |
Subject: |
gmpbench-0.1 and gexpr.c and Win32 |
Date: |
Tue, 08 Mar 2005 05:54:35 +0000 |
User-agent: |
Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.3) Gecko/20030313 |
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);
}
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- gmpbench-0.1 and gexpr.c and Win32,
Sisyphus <=