This is the mail archive of the
`libc-alpha@sourceware.org`
mailing list for the glibc project.

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |

Other format: | [Raw text] |

*From*: Siddhesh Poyarekar <siddhesh at redhat dot com>*To*: libc-alpha at sourceware dot org*Date*: Fri, 21 Dec 2012 00:43:20 +0530*Subject*: [PATCH] Convert mantissa storage in mp_no to int

Hi, The mantissa in mp_no is currently stored as an array of doubles. This is unnecessary and slow since: 1) The largest number stored is 2^24 2) Only integers are stored. The side-effect of the storage being double is that the generated code is all fp arithmetic, which is a lost opportunity for optimization. Attached patch converts this double array to int and adjusts mp code accordingly. This patch is not exhaustive in terms of converting all possible fp operations to integer operations, but it ensures that correctness of operations is maintained. I'll be following up with patches for individual functions to tighten them up and eliminate fp operations wherever possible. I have verified that the testsuite runs clean after this change. Like my last patch, I tested the patch for performance improvement with the help of the old trusty pow function test in the following blog post: http://entropymine.com/imageworsener/slowpow/ and the command line: time ./powtest 1000000 1.0000000000000020 1.5000000050000000 The patch takes a whole second off the execution time for the above command. OK for 2.18? Siddhesh ChangeLog: 2012-12-21 Siddhesh Poyarekar <siddhesh@redhat.com> * sysdeps/ieee754/dbl-64/mpa.c (I_RADIX): New macro. (TWOPOW): Likewise. (__acr): Juggle logic a bit to make it faster and use int instead of doubles. (__mp_dbl): Use int values instead of double. (__dbl_mp): Likewise. (add_magnitudes): Likewise. (sub_magnitudes): Likewise. (__add): Likewise. (__sub): Likewise. (__mul): Likewise. * sysdeps/ieee754/dbl-64/mpa.h (mp_no): Change type of D to int array. * sysdeps/ieee754/dbl-64/mpa.c (__mpatan): Use int values instead of double. * sysdeps/ieee754/dbl-64/mpatan2.c (__mpatan2): Likewise. * sysdeps/ieee754/dbl-64/mpexp.c (__mpexp): Likewise. * sysdeps/ieee754/dbl-64/mpsqrt.c (__mpsqrt): Likewise. * sysdeps/ieee754/dbl-64/mptan.c (__mptan): Likewise. * sysdeps/ieee754/dbl-64/sincos32.c (__c32): Likewise. (__mpranred): Likewise. * sysdeps/ieee754/dbl-64/sincos32.h (oofac27): Likewise. (pi): Likewise. (hp): Likewise. * sysdeps/ieee754/dbl-64/slowpow.c (__slowpow): Likewise. diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c index b5d25ed..5ea13b8 100644 --- a/sysdeps/ieee754/dbl-64/mpa.c +++ b/sysdeps/ieee754/dbl-64/mpa.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2011 Free Software Foundation + * Copyright (C) 2001-2012 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -47,6 +47,9 @@ #include "mpa2.h" #include <sys/param.h> /* For MIN() */ +#define I_RADIX (1 << 24) +#define TWOPOW(n) (1 << n) + #ifndef SECTION # define SECTION #endif @@ -76,16 +79,12 @@ int __acr(const mp_no *x, const mp_no *y, int p) { int i; - if (X[0] == ZERO) { - if (Y[0] == ZERO) i= 0; - else i=-1; - } - else if (Y[0] == ZERO) i= 1; - else { - if (EX > EY) i= 1; - else if (EX < EY) i=-1; - else i= mcr(x,y,p); - } + if (X[0] == 0 || Y[0] == 0) + return X[0] - Y[0]; + + if (EX > EY) i= 1; + else if (EX < EY) i=-1; + else i= mcr(x,y,p); return i; } @@ -245,10 +244,10 @@ void __mp_dbl(const mp_no *x, double *y, int p) { double a,c,u,v,z[5]; #endif - if (X[0] == ZERO) {*y = ZERO; return; } + if (X[0] == 0) {*y = ZERO; return; } if (EX> -42) norm(x,y,p); - else if (EX==-42 && X[1]>=TWO10) norm(x,y,p); + else if (EX==-42 && X[1]>=TWOPOW(10)) norm(x,y,p); else denorm(x,y,p); } #endif @@ -263,24 +262,21 @@ SECTION __dbl_mp(double x, mp_no *y, int p) { int i,n; - double u; /* Sign */ - if (x == ZERO) {Y[0] = ZERO; return; } - else if (x > ZERO) Y[0] = ONE; - else {Y[0] = MONE; x=-x; } + if (x == ZERO) {Y[0] = 0; return; } + else if (x > ZERO) Y[0] = 1; + else {Y[0] = -1; x=-x; } /* Exponent */ - for (EY=ONE; x >= RADIX; EY += ONE) x *= RADIXI; - for ( ; x < ONE; EY -= ONE) x *= RADIX; + for (EY=1; x >= RADIX; EY++) x *= RADIXI; + for ( ; x < ONE; EY--) x *= RADIX; /* Digits */ n=MIN(p,4); for (i=1; i<=n; i++) { - u = (x + TWO52) - TWO52; - if (u>x) u -= ONE; - Y[i] = u; x -= u; x *= RADIX; } - for ( ; i<=p; i++) Y[i] = ZERO; + Y[i] = x; x -= Y[i]; x *= RADIX; } + for ( ; i<=p; i++) Y[i] = 0; } @@ -302,29 +298,29 @@ add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { if (j<1) {__cpy(x,z,p); return; } - else Z[k] = ZERO; + else Z[k] = 0; for (; j>0; i--,j--) { Z[k] += X[i] + Y[j]; - if (Z[k] >= RADIX) { - Z[k] -= RADIX; - Z[--k] = ONE; } + if (Z[k] >= I_RADIX) { + Z[k] -= I_RADIX; + Z[--k] = 1; } else - Z[--k] = ZERO; + Z[--k] = 0; } for (; i>0; i--) { Z[k] += X[i]; - if (Z[k] >= RADIX) { - Z[k] -= RADIX; - Z[--k] = ONE; } + if (Z[k] >= I_RADIX) { + Z[k] -= I_RADIX; + Z[--k] = 1; } else - Z[--k] = ZERO; + Z[--k] = 0; } - if (Z[1] == ZERO) { + if (Z[1] == 0) { for (i=1; i<=p; i++) Z[i] = Z[i+1]; } - else EZ += ONE; + else EZ++; } @@ -344,45 +340,45 @@ sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { if (EX == EY) { i=j=k=p; - Z[k] = Z[k+1] = ZERO; } + Z[k] = Z[k+1] = 0; } else { j= EX - EY; if (j > p) {__cpy(x,z,p); return; } else { i=p; j=p+1-j; k=p; - if (Y[j] > ZERO) { - Z[k+1] = RADIX - Y[j--]; - Z[k] = MONE; } + if (Y[j] > 0) { + Z[k+1] = I_RADIX - Y[j--]; + Z[k] = -1; } else { - Z[k+1] = ZERO; - Z[k] = ZERO; j--;} + Z[k+1] = 0; + Z[k] = 0; j--;} } } for (; j>0; i--,j--) { Z[k] += (X[i] - Y[j]); - if (Z[k] < ZERO) { - Z[k] += RADIX; - Z[--k] = MONE; } + if (Z[k] < 0) { + Z[k] += I_RADIX; + Z[--k] = -1; } else - Z[--k] = ZERO; + Z[--k] = 0; } for (; i>0; i--) { Z[k] += X[i]; - if (Z[k] < ZERO) { - Z[k] += RADIX; - Z[--k] = MONE; } + if (Z[k] < 0) { + Z[k] += I_RADIX; + Z[--k] = -1; } else - Z[--k] = ZERO; + Z[--k] = 0; } - for (i=1; Z[i] == ZERO; i++) ; + for (i=1; Z[i] == 0; i++) ; EZ = EZ - i + 1; for (k=1; i <= p+1; ) Z[k++] = Z[i++]; for (; k <= p; ) - Z[k++] = ZERO; + Z[k++] = 0; } @@ -396,8 +392,8 @@ __add(const mp_no *x, const mp_no *y, mp_no *z, int p) { int n; - if (X[0] == ZERO) {__cpy(y,z,p); return; } - else if (Y[0] == ZERO) {__cpy(x,z,p); return; } + if (X[0] == 0) {__cpy(y,z,p); return; } + else if (Y[0] == 0) {__cpy(x,z,p); return; } if (X[0] == Y[0]) { if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; } @@ -406,7 +402,7 @@ __add(const mp_no *x, const mp_no *y, mp_no *z, int p) { else { if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; } else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; } - else Z[0] = ZERO; + else Z[0] = 0; } } @@ -421,8 +417,8 @@ __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { int n; - if (X[0] == ZERO) {__cpy(y,z,p); Z[0] = -Z[0]; return; } - else if (Y[0] == ZERO) {__cpy(x,z,p); return; } + if (X[0] == 0) {__cpy(y,z,p); Z[0] = -Z[0]; return; } + else if (Y[0] == 0) {__cpy(x,z,p); return; } if (X[0] != Y[0]) { if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; } @@ -431,7 +427,7 @@ __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { else { if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; } else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = -Y[0]; } - else Z[0] = ZERO; + else Z[0] = 0; } } @@ -446,28 +442,26 @@ SECTION __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { int i, i1, i2, j, k, k2; - double u; /* Is z=0? */ - if (X[0]*Y[0]==ZERO) - { Z[0]=ZERO; return; } + if (X[0]*Y[0] == 0) + { Z[0]=0; return; } /* Multiply, add and carry */ k2 = (p<3) ? p+p : p+3; - Z[k2]=ZERO; + Z[k2]=0; for (k=k2; k>1; ) { if (k > p) {i1=k-p; i2=p+1; } else {i1=1; i2=k; } - for (i=i1,j=i2-1; i<i2; i++,j--) Z[k] += X[i]*Y[j]; + long tmp = Z[k]; + for (i=i1,j=i2-1; i<i2; i++,j--) tmp += (long) X[i]*Y[j]; - u = (Z[k] + CUTTER)-CUTTER; - if (u > Z[k]) u -= RADIX; - Z[k] -= u; - Z[--k] = u*RADIXI; + Z[k] = (int) (tmp % (1 << 24)); + Z[--k] = (int) (tmp / (1 << 24)); } /* Is there a carry beyond the most significant digit? */ - if (Z[1] == ZERO) { + if (Z[1] == 0) { for (i=1; i<=p; i++) Z[i]=Z[i+1]; EZ = EX + EY - 1; } else diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h index 4fdecb6..c09689a 100644 --- a/sysdeps/ieee754/dbl-64/mpa.h +++ b/sysdeps/ieee754/dbl-64/mpa.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001, 2011 Free Software Foundation, Inc. + * Copyright (C) 2001-2012 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -41,7 +41,7 @@ typedef struct {/* This structure holds the details of a multi-precision */ int e; /* floating point number, x: d[0] holds its sign (-1,0 or 1) */ - double d[40]; /* e holds its exponent (...,-2,-1,0,1,2,...) and */ + int d[40]; /* e holds its exponent (...,-2,-1,0,1,2,...) and */ } mp_no; /* d[1]...d[p] hold its mantissa digits. The value of x is, */ /* x = d[1]*r**(e-1) + d[2]*r**(e-2) + ... + d[p]*r**(e-p). */ /* Here r = 2**24, 0 <= d[i] < r and 1 <= p <= 32. */ diff --git a/sysdeps/ieee754/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c index d897bbb..9119e58 100644 --- a/sysdeps/ieee754/dbl-64/mpatan.c +++ b/sysdeps/ieee754/dbl-64/mpatan.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2011 Free Software Foundation + * Copyright (C) 2001-2012 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -48,12 +48,12 @@ __mpatan(mp_no *x, mp_no *y, int p) { int i,m,n; double dx; mp_no - mptwo = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, - mptwoim1 = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + mptwo = {0,{0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0}}, + mptwoim1 = {0,{0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0}}; mp_no mps,mpsm,mpt,mpt1,mpt2,mpt3; @@ -66,8 +66,8 @@ __mpatan(mp_no *x, mp_no *y, int p) { {if (dx>__atan_xm[m].d) break;} } mptwo.e = mptwoim1.e = 1; - mptwo.d[0] = mptwoim1.d[0] = ONE; - mptwo.d[1] = TWO; + mptwo.d[0] = mptwoim1.d[0] = 1; + mptwo.d[1] = 2; /* Reduce x m times */ __mul(x,x,&mpsm,p); @@ -89,7 +89,7 @@ __mpatan(mp_no *x, mp_no *y, int p) { n=__atan_np[p]; mptwoim1.d[1] = __atan_twonm1[p].d; __dvd(&mpsm,&mptwoim1,&mpt,p); for (i=n-1; i>1; i--) { - mptwoim1.d[1] -= TWO; + mptwoim1.d[1] -= 2; __dvd(&mpsm,&mptwoim1,&mpt1,p); __mul(&mpsm,&mpt,&mpt2,p); __sub(&mpt1,&mpt2,&mpt,p); diff --git a/sysdeps/ieee754/dbl-64/mpatan2.c b/sysdeps/ieee754/dbl-64/mpatan2.c index 42cb6a7..3e3252b 100644 --- a/sysdeps/ieee754/dbl-64/mpatan2.c +++ b/sysdeps/ieee754/dbl-64/mpatan2.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2011 Free Software Foundation + * Copyright (C) 2001-2012 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -49,14 +49,12 @@ void SECTION __mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) { - static const double ZERO = 0.0, ONE = 1.0; - mp_no mpt1,mpt2,mpt3; - if (X[0] <= ZERO) { + if (X[0] <= 0) { __dvd(x,y,&mpt1,p); __mul(&mpt1,&mpt1,&mpt2,p); - if (mpt1.d[0] != ZERO) mpt1.d[0] = ONE; + if (mpt1.d[0] != 0) mpt1.d[0] = 1; __add(&mpt2,&mpone,&mpt3,p); __mpsqrt(&mpt3,&mpt2,p); __add(&mpt1,&mpt2,&mpt3,p); mpt3.d[0]=Y[0]; __mpatan(&mpt3,&mpt1,p); __add(&mpt1,&mpt1,z,p); diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c index 2762ba9..87391dc 100644 --- a/sysdeps/ieee754/dbl-64/mpexp.c +++ b/sysdeps/ieee754/dbl-64/mpexp.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2011 Free Software Foundation + * Copyright (C) 2001-2012 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -56,9 +56,9 @@ __mpexp(mp_no *x, mp_no *y, int p) { { 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63}, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}}; - mp_no mpk = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + mp_no mpk = {0,{0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0}}; mp_no mps,mpak,mpt1,mpt2; /* Choose m,n and compute a=2**(-m) */ @@ -81,7 +81,7 @@ __mpexp(mp_no *x, mp_no *y, int p) { __mul(x,&mpt1,&mps,p); /* Evaluate the polynomial. Put result in mpt2 */ - mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=__mpexp_nn[n].d; + mpk.e = 1; mpk.d[0] = 1; mpk.d[1]=__mpexp_nn[n].d; __dvd(&mps,&mpk,&mpt1,p); __add(&mpone,&mpt1,&mpak,p); for (k=n-1; k>1; k--) { diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c index 92bf5ef..e29105b 100644 --- a/sysdeps/ieee754/dbl-64/mpsqrt.c +++ b/sysdeps/ieee754/dbl-64/mpsqrt.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2011 Free Software Foundation + * Copyright (C) 2001-2012 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -53,17 +53,17 @@ __mpsqrt(mp_no *x, mp_no *y, int p) { int i,m,ey; double dx,dy; mp_no - mphalf = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, - mp3halfs = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + mphalf = {0,{0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0}}, + mp3halfs = {0,{0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0}}; mp_no mpxn,mpz,mpu,mpt1,mpt2; /* Prepare multi-precision 1/2 and 3/2 */ - mphalf.e =0; mphalf.d[0] =ONE; mphalf.d[1] =HALFRAD; - mp3halfs.e=1; mp3halfs.d[0]=ONE; mp3halfs.d[1]=ONE; mp3halfs.d[2]=HALFRAD; + mphalf.e =0; mphalf.d[0] =1; mphalf.d[1] =(1 << 23); + mp3halfs.e=1; mp3halfs.d[0]=1; mp3halfs.d[1]=1; mp3halfs.d[2]=(1 << 23); ey=EX/2; __cpy(x,&mpxn,p); mpxn.e -= (ey+ey); __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p); diff --git a/sysdeps/ieee754/dbl-64/mptan.c b/sysdeps/ieee754/dbl-64/mptan.c index e75b2da..356bc04 100644 --- a/sysdeps/ieee754/dbl-64/mptan.c +++ b/sysdeps/ieee754/dbl-64/mptan.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2011 Free Software Foundation + * Copyright (C) 2001-2012 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -47,8 +47,6 @@ void SECTION __mptan(double x, mp_no *mpy, int p) { - static const double MONE = -1.0; - int n; mp_no mpw, mpc, mps; @@ -56,7 +54,7 @@ __mptan(double x, mp_no *mpy, int p) { __c32(&mpw, &mpc, &mps, p); /* computing sin(x) and cos(x) */ if (n) /* second or fourth quarter of unit circle */ { __dvd(&mpc,&mps,mpy,p); - mpy->d[0] *= MONE; + mpy->d[0] *= -1; } /* tan is negative in this area */ else __dvd(&mps,&mpc,mpy,p); diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c index 038d1bd..3890815 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.c +++ b/sysdeps/ieee754/dbl-64/sincos32.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2011 Free Software Foundation + * Copyright (C) 2001-2012 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -119,7 +119,7 @@ cc32(mp_no *x, mp_no *y, int p) { void SECTION __c32(mp_no *x, mp_no *y, mp_no *z, int p) { - static const mp_no mpt={1,{1.0,2.0}}, one={1,{1.0,1.0}}; + static const mp_no mpt={1,{1,2}}, one={1,{1,1}}; mp_no u,t,t1,t2,c,s; int i; __cpy(x,&u,p); @@ -251,7 +251,7 @@ __mpranred(double x, mp_no *y, int p) number v; double t,xn; int i,k,n; - static const mp_no one = {1,{1.0,1.0}}; + static const mp_no one = {1,{1,1}}; mp_no a,b,c; if (ABS(x) < 2.8e14) { @@ -267,18 +267,18 @@ __mpranred(double x, mp_no *y, int p) } else { /* if x is very big more precision required */ __dbl_mp(x,&a,p); - a.d[0]=1.0; + a.d[0]=1; k = a.e-5; if (k < 0) k=0; b.e = -k; - b.d[0] = 1.0; + b.d[0] = 1; for (i=0;i<p;i++) b.d[i+1] = toverp[i+k]; __mul(&a,&b,&c,p); t = c.d[c.e]; for (i=1;i<=p-c.e;i++) c.d[i]=c.d[i+c.e]; for (i=p+1-c.e;i<=p;i++) c.d[i]=0; c.e=0; - if (c.d[1] >= 8388608.0) + if (c.d[1] >= 8388608) { t +=1.0; __sub(&c,&one,&b,p); __mul(&b,&hp,y,p); diff --git a/sysdeps/ieee754/dbl-64/sincos32.h b/sysdeps/ieee754/dbl-64/sincos32.h index 6efe4d4..14562fb 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.h +++ b/sysdeps/ieee754/dbl-64/sincos32.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001-2012 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -42,24 +42,24 @@ static const number #endif static const mp_no - oofac27 = {-3,{1.0,7.0,4631664.0,12006312.0,13118056.0,6538613.0,646354.0, - 8508025.0,9131256.0,7548776.0,2529842.0,8864927.0,660489.0,15595125.0,12777885.0, - 11618489.0,13348664.0,5486686.0,514518.0,11275535.0,4727621.0,3575562.0, - 13579710.0,5829745.0,7531862.0,9507898.0,6915060.0,4079264.0,1907586.0, - 6078398.0,13789314.0,5504104.0,14136.0}}, - pi = {1,{1.0,3.0, - 2375530.0,8947107.0,578323.0,1673774.0,225395.0,4498441.0,3678761.0, - 10432976.0,536314.0,10021966.0,7113029.0,2630118.0,3723283.0,7847508.0, - 6737716.0,15273068.0,12626985.0,12044668.0,5299519.0,8705461.0,11880201.0, - 1544726.0,14014857.0,7994139.0,13709579.0,10918111.0,11906095.0,16610011.0, - 13638367.0,12040417.0,11529578.0,2522774.0}}, - hp = {1,{1.0, 1.0, - 9576373.0,4473553.0,8677769.0,9225495.0,112697.0,10637828.0, - 10227988.0,13605096.0,268157.0,5010983.0,3556514.0,9703667.0, - 1861641.0,12312362.0,3368858.0,7636534.0,6313492.0,14410942.0, - 2649759.0,12741338.0,14328708.0,9160971.0,7007428.0,12385677.0, - 15243397.0,13847663.0,14341655.0,16693613.0,15207791.0,14408816.0, - 14153397.0,1261387.0,6110792.0,2291862.0,4181138.0,5295267.0}}; + oofac27 = {-3,{1,7,4631664,12006312,13118056,6538613,646354, + 8508025,9131256,7548776,2529842,8864927,660489,15595125,12777885, + 11618489,13348664,5486686,514518,11275535,4727621,3575562, + 13579710,5829745,7531862,9507898,6915060,4079264,1907586, + 6078398,13789314,5504104,14136}}, + pi = {1,{1,3, + 2375530,8947107,578323,1673774,225395,4498441,3678761, + 10432976,536314,10021966,7113029,2630118,3723283,7847508, + 6737716,15273068,12626985,12044668,5299519,8705461,11880201, + 1544726,14014857,7994139,13709579,10918111,11906095,16610011, + 13638367,12040417,11529578,2522774}}, + hp = {1,{1, 1, + 9576373,4473553,8677769,9225495,112697,10637828, + 10227988,13605096,268157,5010983,3556514,9703667, + 1861641,12312362,3368858,7636534,6313492,14410942, + 2649759,12741338,14328708,9160971,7007428,12385677, + 15243397,13847663,14341655,16693613,15207791,14408816, + 14153397,1261387,6110792,2291862,4181138,5295267}}; static const double toverp[75] = { 10680707.0, 7228996.0, 1387004.0, 2578385.0, 16069853.0, diff --git a/sysdeps/ieee754/dbl-64/slowpow.c b/sysdeps/ieee754/dbl-64/slowpow.c index 7c829c7..be2808f 100644 --- a/sysdeps/ieee754/dbl-64/slowpow.c +++ b/sysdeps/ieee754/dbl-64/slowpow.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2011 Free Software Foundation + * Copyright (C) 2001-2012 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -48,7 +48,7 @@ SECTION __slowpow(double x, double y, double z) { double res,res1; mp_no mpx, mpy, mpz,mpw,mpp,mpr,mpr1; - static const mp_no eps = {-3,{1.0,4.0}}; + static const mp_no eps = {-3,{1,4}}; int p; res = __halfulp(x,y); /* halfulp() returns -10 or x^y */

**Follow-Ups**:**Re: [PATCH] Convert mantissa storage in mp_no to int***From:*Joseph S. Myers

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |