{ int p;
mp_no mpx,mpy,mpz,mperr,mpz2,mpt1;
p=6;
- dbl_mp(ax,&mpx,p); dbl_mp(ay,&mpy,p); dvd(&mpy,&mpx,&mpz,p);
- dbl_mp(ue.d,&mpt1,p); mul(&mpz,&mpt1,&mperr,p);
+ __dbl_mp(ax,&mpx,p); __dbl_mp(ay,&mpy,p); dvd(&mpy,&mpx,&mpz,p);
+ __dbl_mp(ue.d,&mpt1,p); mul(&mpz,&mpt1,&mperr,p);
sub(&mpz,&mperr,&mpz2,p); __mp_dbl(&mpz2,&z,p);
return signArctan2(y,z);
}
mp_no mpx,mpy,mpz,mpz1,mpz2,mperr,mpt1;
for (i=0; i<MM; i++) {
p = pr[i];
- dbl_mp(x,&mpx,p); dbl_mp(y,&mpy,p);
+ __dbl_mp(x,&mpx,p); __dbl_mp(y,&mpy,p);
__mpatan2(&mpy,&mpx,&mpz,p);
- dbl_mp(ud[i].d,&mpt1,p); mul(&mpz,&mpt1,&mperr,p);
+ __dbl_mp(ud[i].d,&mpt1,p); mul(&mpz,&mpt1,&mperr,p);
add(&mpz,&mperr,&mpz1,p); sub(&mpz,&mperr,&mpz2,p);
__mp_dbl(&mpz1,&z1,p); __mp_dbl(&mpz2,&z2,p);
if (z1==z2) return z1;
for (i=0; i<M; i++) {
p = pr[i];
- dbl_mp(x,&mpx,p); dbl_mp(y,&mpy,p);
+ __dbl_mp(x,&mpx,p); __dbl_mp(y,&mpy,p);
__mplog(&mpx,&mpy,p);
- dbl_mp(e[i].d,&mperr,p);
+ __dbl_mp(e[i].d,&mperr,p);
add(&mpy,&mperr,&mpy1,p); sub(&mpy,&mperr,&mpy2,p);
__mp_dbl(&mpy1,&y1,p); __mp_dbl(&mpy2,&y2,p);
if (y1==y2) return y1;
/* number *y. If the precision p is too small the result is truncated. x is */
/* left unchanged. */
-void dbl_mp(double x, mp_no *y, int p) {
+void __dbl_mp(double x, mp_no *y, int p) {
int i,n;
double u;
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
cpy(x,&z,p); z.e=0; __mp_dbl(&z,&t,p);
- t=ONE/t; dbl_mp(t,y,p); EY -= EX;
+ t=ONE/t; __dbl_mp(t,y,p); EY -= EX;
for (i=0; i<np1[p]; i++) {
cpy(y,&w,p);
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define ABS(x) ((x) < 0 ? -(x) : (x))
-int acr(const mp_no *, const mp_no *, int);
-int cr(const mp_no *, const mp_no *, int);
-void cpy(const mp_no *, mp_no *, int);
-void cpymn(const mp_no *, int, mp_no *, int);
+int __acr(const mp_no *, const mp_no *, int);
+int __cr(const mp_no *, const mp_no *, int);
+void __cpy(const mp_no *, mp_no *, int);
+void __cpymn(const mp_no *, int, mp_no *, int);
void __mp_dbl(const mp_no *, double *, int);
-void dbl_mp(double, mp_no *, int);
-void add(const mp_no *, const mp_no *, mp_no *, int);
-void sub(const mp_no *, const mp_no *, mp_no *, int);
-void mul(const mp_no *, const mp_no *, mp_no *, int);
-void inv(const mp_no *, mp_no *, int);
-void dvd(const mp_no *, const mp_no *, mp_no *, int);
+void __dbl_mp(double, mp_no *, int);
+void __add(const mp_no *, const mp_no *, mp_no *, int);
+void __sub(const mp_no *, const mp_no *, mp_no *, int);
+void __mul(const mp_no *, const mp_no *, mp_no *, int);
+void __inv(const mp_no *, mp_no *, int);
+void __dvd(const mp_no *, const mp_no *, mp_no *, int);
#include "endian.h"
#include "mpa.h"
-void mpsqrt(mp_no *, mp_no *, int);
+void __mpsqrt(mp_no *, mp_no *, int);
-void mpatan(mp_no *x, mp_no *y, int p) {
+void __mpatan(mp_no *x, mp_no *y, int p) {
#include "mpatan.h"
int i,m,n;
else {
for (i=0; i<m; i++) {
add(&mpone,&mpsm,&mpt1,p);
- mpsqrt(&mpt1,&mpt2,p);
+ __mpsqrt(&mpt1,&mpt2,p);
add(&mpt2,&mpt2,&mpt1,p);
add(&mptwo,&mpsm,&mpt2,p);
add(&mpt1,&mpt2,&mpt3,p);
dvd(&mpsm,&mpt3,&mpt1,p);
cpy(&mpt1,&mpsm,p);
}
- mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0];
+ __mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0];
}
/* Evaluate a truncated power series for Atan(s) */
#include "mpa.h"
-void mpsqrt(mp_no *, mp_no *, int);
-void mpatan(mp_no *, mp_no *, int);
+void __mpsqrt(mp_no *, mp_no *, int);
+void __mpatan(mp_no *, mp_no *, int);
/* Multi-Precision Atan2(y,x) function subroutine, for p >= 4. */
/* y=0 is not permitted if x<=0. No error messages are given. */
-void mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) {
+void __mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) {
static const double ZERO = 0.0, ONE = 1.0;
mpone.e = 1; mpone.d[0] = mpone.d[1] = ONE;
dvd(x,y,&mpt1,p); mul(&mpt1,&mpt1,&mpt2,p);
if (mpt1.d[0] != ZERO) mpt1.d[0] = ONE;
- add(&mpt2,&mpone,&mpt3,p); mpsqrt(&mpt3,&mpt2,p);
+ 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);
+ __mpatan(&mpt3,&mpt1,p); add(&mpt1,&mpt1,z,p);
}
else
{ dvd(y,x,&mpt1,p);
- mpatan(&mpt1,z,p);
+ __mpatan(&mpt1,z,p);
}
return;
}
/* Compute s=x*2**(-m). Put result in mps */
- dbl_mp(a,&mpt1,p);
+ __dbl_mp(a,&mpt1,p);
mul(x,&mpt1,&mps,p);
/* Evaluate the polynomial. Put result in mpt2 */
#include "endian.h"
#include "mpa.h"
-void mpexp(mp_no *, mp_no *, int);
+void __mpexp(mp_no *, mp_no *, int);
-void mplog(mp_no *x, mp_no *y, int p) {
+void __mplog(mp_no *x, mp_no *y, int p) {
#include "mplog.h"
int i,m;
#if 0
cpy(y,&mpt1,p);
for (i=0; i<m; i++) {
mpt1.d[0]=-mpt1.d[0];
- mpexp(&mpt1,&mpt2,p);
+ __mpexp(&mpt1,&mpt2,p);
mul(x,&mpt2,&mpt1,p);
sub(&mpt1,&mpone,&mpt2,p);
add(y,&mpt2,&mpt1,p);
double fastiroot(double);
-void mpsqrt(mp_no *x, mp_no *y, int p) {
+void __mpsqrt(mp_no *x, mp_no *y, int p) {
#include "mpsqrt.h"
int i,m,ex,ey;
mp3halfs.e=1; mp3halfs.d[0]=ONE; mp3halfs.d[1]=ONE; mp3halfs.d[2]=HALFRAD;
ex=EX; ey=EX/2; cpy(x,&mpxn,p); mpxn.e -= (ey+ey);
- __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); dbl_mp(dy,&mpu,p);
+ __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p);
mul(&mpxn,&mphalf,&mpz,p);
m=mp[p];
for (i=0; i<M; i++) {
p = pr[i];
- dbl_mp(x,&mpx,p); __mpatan(&mpx,&mpy,p);
- dbl_mp(u9[i].d,&mpt1,p); mul(&mpy,&mpt1,&mperr,p);
- add(&mpy,&mperr,&mpy1,p); sub(&mpy,&mperr,&mpy2,p);
+ __dbl_mp(x,&mpx,p); __mpatan(&mpx,&mpy,p);
+ __dbl_mp(u9[i].d,&mpt1,p); __mul(&mpy,&mpt1,&mperr,p);
+ __add(&mpy,&mperr,&mpy1,p); __sub(&mpy,&mperr,&mpy2,p);
__mp_dbl(&mpy1,&y1,p); __mp_dbl(&mpy2,&y2,p);
if (y1==y2) return y1;
}
/* Second stage */
/* Reduction by algorithm iv */
p=10; n = (mpranred(x,&mpa,p)) & 0x00000001;
- __mp_dbl(&mpa,&a,p); dbl_mp(a,&mpt1,p);
+ __mp_dbl(&mpa,&a,p); __dbl_mp(a,&mpt1,p);
sub(&mpa,&mpt1,&mpt2,p); __mp_dbl(&mpt2,&da,p);
MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
int p;
mp_no a,b,c;
p=32;
- dbl_mp(res,&a,p);
- dbl_mp(0.5*(res1-res),&b,p);
+ __dbl_mp(res,&a,p);
+ __dbl_mp(0.5*(res1-res),&b,p);
add(&a,&b,&c,p);
if (x>0.8)
{ sub(&hp,&c,&a,p);
__c32(&a,&b,&c,p);
}
else __c32(&c,&a,&b,p); /* b=sin(0.5*(res+res1)) */
- dbl_mp(x,&c,p); /* c = x */
+ __dbl_mp(x,&c,p); /* c = x */
sub(&b,&c,&a,p);
/* if a>0 return min(res,res1), otherwise return max(res,res1) */
if (a.d[0]>0) return (res<res1)?res:res1;
int p;
mp_no a,b,c;
p=32;
- dbl_mp(res,&a,p);
- dbl_mp(0.5*(res1-res),&b,p);
+ __dbl_mp(res,&a,p);
+ __dbl_mp(0.5*(res1-res),&b,p);
add(&a,&b,&c,p);
if (x>2.4)
{ sub(&pi,&c,&a,p);
__c32(&a,&c,&b,p);
}
else __c32(&c,&b,&a,p); /* b=cos(0.5*(res+res1)) */
- dbl_mp(x,&c,p); /* c = x */
+ __dbl_mp(x,&c,p); /* c = x */
sub(&b,&c,&a,p);
/* if a>0 return max(res,res1), otherwise return min(res,res1) */
if (a.d[0]>0) return (res>res1)?res:res1;
double y;
mp_no a,b,c;
p=32;
- dbl_mp(x,&a,p);
- dbl_mp(dx,&b,p);
+ __dbl_mp(x,&a,p);
+ __dbl_mp(dx,&b,p);
add(&a,&b,&c,p);
if (x>0.8) { sub(&hp,&c,&a,p); __c32(&a,&b,&c,p); }
else __c32(&c,&a,&b,p); /* b = sin(x+dx) */
double y;
mp_no a,b,c;
p=32;
- dbl_mp(x,&a,p);
- dbl_mp(dx,&b,p);
+ __dbl_mp(x,&a,p);
+ __dbl_mp(dx,&b,p);
add(&a,&b,&c,p);
if (x>0.8)
{ sub(&hp,&c,&b,p);
xn = t - toint.d;
v.d = t;
n =v.i[LOW_HALF]&3;
- dbl_mp(xn,&a,p);
+ __dbl_mp(xn,&a,p);
mul(&a,&hp,&b,p);
- dbl_mp(x,&c,p);
+ __dbl_mp(x,&c,p);
sub(&c,&b,y,p);
return n;
}
else { /* if x is very big more precision required */
- dbl_mp(x,&a,p);
+ __dbl_mp(x,&a,p);
a.d[0]=1.0;
k = a.e-5;
if (k < 0) k=0;
/**************************************************************************/
#include "mpa.h"
-void mpexp(mp_no *x, mp_no *y, int p);
+void __mpexp(mp_no *x, mp_no *y, int p);
/*Converting from double precision to Multi-precision and calculating e^x */
-double slowexp(double x) {
+double __slowexp(double x) {
double w,z,res,eps=3.0e-26;
#if 0
double y;
mp_no mpx, mpy, mpz,mpw,mpeps,mpcor;
p=6;
- dbl_mp(x,&mpx,p); /* Convert a double precision number x */
+ __dbl_mp(x,&mpx,p); /* Convert a double precision number x */
/* into a multiple precision number mpx with prec. p. */
- mpexp(&mpx, &mpy, p); /* Multi-Precision exponential function */
- dbl_mp(eps,&mpeps,p);
+ __mpexp(&mpx, &mpy, p); /* Multi-Precision exponential function */
+ __dbl_mp(eps,&mpeps,p);
mul(&mpeps,&mpy,&mpcor,p);
add(&mpy,&mpcor,&mpw,p);
sub(&mpy,&mpcor,&mpz,p);
if (w == z) return w;
else { /* if calculating is not exactly */
p = 32;
- dbl_mp(x,&mpx,p);
- mpexp(&mpx, &mpy, p);
+ __dbl_mp(x,&mpx,p);
+ __mpexp(&mpx, &mpy, p);
__mp_dbl(&mpy, &res, p);
return res;
}
#include "mpa.h"
-void mpexp(mp_no *x, mp_no *y, int p);
-void mplog(mp_no *x, mp_no *y, int p);
+void __mpexp(mp_no *x, mp_no *y, int p);
+void __mplog(mp_no *x, mp_no *y, int p);
double ulog(double);
double halfulp(double x,double y);
if (res >= 0) return res; /* if result was really computed by halfulp */
/* else, if result was not really computed by halfulp */
p = 10; /* p=precision */
- dbl_mp(x,&mpx,p);
- dbl_mp(y,&mpy,p);
- dbl_mp(z,&mpz,p);
- mplog(&mpx, &mpz, p); /* log(x) = z */
+ __dbl_mp(x,&mpx,p);
+ __dbl_mp(y,&mpy,p);
+ __dbl_mp(z,&mpz,p);
+ __mplog(&mpx, &mpz, p); /* log(x) = z */
mul(&mpy,&mpz,&mpw,p); /* y * z =w */
- mpexp(&mpw, &mpp, p); /* e^w =pp */
+ __mpexp(&mpw, &mpp, p); /* e^w =pp */
add(&mpp,&eps,&mpr,p); /* pp+eps =r */
__mp_dbl(&mpr, &res, p);
sub(&mpp,&eps,&mpr1,p); /* pp -eps =r1 */
if (res == res1) return res;
p = 32; /* if we get here result wasn't calculated exactly, continue */
- dbl_mp(x,&mpx,p); /* for more exact calculation */
- dbl_mp(y,&mpy,p);
- dbl_mp(z,&mpz,p);
- mplog(&mpx, &mpz, p); /* log(c)=z */
+ __dbl_mp(x,&mpx,p); /* for more exact calculation */
+ __dbl_mp(y,&mpy,p);
+ __dbl_mp(z,&mpz,p);
+ __mplog(&mpx, &mpz, p); /* log(c)=z */
mul(&mpy,&mpz,&mpw,p); /* y*z =w */
- mpexp(&mpw, &mpp, p); /* e^w=pp */
+ __mpexp(&mpw, &mpp, p); /* e^w=pp */
__mp_dbl(&mpp, &res, p); /* converting into double precision */
return res;
}