mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-25 22:40:05 +00:00
Update.
2001-03-12 Ulrich Drepper <drepper@redhat.com> * sysdeps/ieee754/dbl-64/e_remainder.c: Fix handling of boundary conditions. * sysdeps/ieee754/dbl-64/e_pow.c: Fix handling of boundary conditions. * sysdeps/ieee754/dbl-64/s_sin.c (__sin): Handle Inf and NaN correctly. (__cos): Likewise. * sysdeps/ieee754/dbl-64/e_asin.c (__ieee754_asin): Handle NaN correctly. (__ieee754_acos): Likewise. redefinition. * sysdeps/ieee754/dbl-64/endian.h: Define also one of BIG_ENDI and LITTLE_ENDI. * sysdeps/ieee754/dbl-64/MathLib.h (Init_Lib): Use void as parameter list.
This commit is contained in:
Ulrich Drepper
parent
9a656848ea
commit
ca58f1dbeb
24
ChangeLog
24
ChangeLog
@ -1,10 +1,26 @@
|
||||
2001-03-12 Ulrich Drepper <drepper@redhat.com>
|
||||
|
||||
* sysdeps/ieee754/dbl-64/e_remainder.c: Fix handling of boundary
|
||||
conditions.
|
||||
|
||||
* sysdeps/ieee754/dbl-64/e_pow.c: Fix handling of boundary
|
||||
conditions.
|
||||
|
||||
* sysdeps/ieee754/dbl-64/s_sin.c (__sin): Handle Inf and NaN
|
||||
correctly.
|
||||
(__cos): Likewise.
|
||||
|
||||
* sysdeps/ieee754/dbl-64/e_asin.c (__ieee754_asin): Handle NaN
|
||||
correctly.
|
||||
(__ieee754_acos): Likewise.
|
||||
|
||||
2001-03-12 Andreas Jaeger <aj@suse.de>
|
||||
|
||||
* sysdeps/unix/sysv/linux/s390/sysdep.h (_LINUX_S390_SYSDEP_H):
|
||||
Fix typo. Patch by Martin Schwidefsky <schwidefsky@de.ibm.com>.
|
||||
|
||||
* sysdeps/s390/bits/string.h: Protect __STRING_INLINE against
|
||||
redefinition.
|
||||
redefinition.
|
||||
|
||||
2001-03-11 Roland McGrath <roland@frob.com>
|
||||
|
||||
@ -12,6 +28,12 @@
|
||||
|
||||
2001-03-11 Ulrich Drepper <drepper@redhat.com>
|
||||
|
||||
* sysdeps/ieee754/dbl-64/endian.h: Define also one of BIG_ENDI and
|
||||
LITTLE_ENDI.
|
||||
|
||||
* sysdeps/ieee754/dbl-64/MathLib.h (Init_Lib): Use void as
|
||||
parameter list.
|
||||
|
||||
Last-bit accurate math library implementation by IBM Haifa.
|
||||
Contributed by Abraham Ziv <ziv@il.ibm.com>, Moshe Olshansky
|
||||
<olshansk@il.ibm.com>, Ealan Henis <ealan@il.ibm.com>, and
|
||||
|
||||
@ -45,7 +45,7 @@
|
||||
/* x=n*pi/2+(a+aa), abs(a+aa)<pi/4, n=0,+-1,+-2,.... */
|
||||
/* Routine return integer (n mod 4) */
|
||||
/*******************************************************************/
|
||||
int branred(double x, double *a, double *aa)
|
||||
int __branred(double x, double *a, double *aa)
|
||||
{
|
||||
int i,k;
|
||||
#if 0
|
||||
|
||||
@ -5,9 +5,9 @@
|
||||
*
|
||||
* 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
|
||||
* the Free Software Foundation; either version 2 of the License, or
|
||||
* the Free Software Foundation; either version 2 of the License, 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
|
||||
@ -15,12 +15,12 @@
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License
|
||||
* along with this program; if not, write to the Free Software
|
||||
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
|
||||
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
|
||||
*/
|
||||
/**********************************************************************/
|
||||
/* MODULE_NAME: doasin.c */
|
||||
/* */
|
||||
/* FUNCTION: doasin */
|
||||
/* FUNCTION: doasin */
|
||||
/* */
|
||||
/* FILES NEEDED:endian.h mydefs.h dla.h doasin.h */
|
||||
/* mpa.c */
|
||||
@ -31,17 +31,17 @@
|
||||
|
||||
#include "endian.h"
|
||||
#include "mydefs.h"
|
||||
#include "dla.h"
|
||||
#include "dla.h"
|
||||
|
||||
/********************************************************************/
|
||||
/* Compute arcsin(x,dx,v) of double-length number (x+dx) the result */
|
||||
/* stored in v where v= v[0]+v[1] =arcsin(x+dx) */
|
||||
/********************************************************************/
|
||||
void doasin(double x, double dx, double v[]) {
|
||||
|
||||
void __doasin(double x, double dx, double v[]) {
|
||||
|
||||
#include "doasin.h"
|
||||
|
||||
static const double
|
||||
static const double
|
||||
d5 = 0.22372159090911789889975459505194491E-01,
|
||||
d6 = 0.17352764422456822913014975683014622E-01,
|
||||
d7 = 0.13964843843786693521653681033981614E-01,
|
||||
@ -49,10 +49,10 @@ void doasin(double x, double dx, double v[]) {
|
||||
d9 = 0.97622386568166960207425666787248914E-02,
|
||||
d10 = 0.83638737193775788576092749009744976E-02,
|
||||
d11 = 0.79470250400727425881446981833568758E-02;
|
||||
|
||||
|
||||
double xx,p,pp,u,uu,r,s;
|
||||
double hx,tx,hy,ty,tp,tq,tc,tcc;
|
||||
|
||||
|
||||
|
||||
/* Taylor series for arcsin for Double-Length numbers */
|
||||
xx = x*x+2.0*x*dx;
|
||||
@ -73,9 +73,3 @@ void doasin(double x, double dx, double v[]) {
|
||||
v[0]=p;
|
||||
v[1]=pp; /* arcsin(x+dx)=v[0]+v[1] */
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@ -45,7 +45,7 @@
|
||||
/*(x+dx) between 0 and PI/4 */
|
||||
/***********************************************************************/
|
||||
|
||||
void dubsin(double x, double dx, double v[]) {
|
||||
void __dubsin(double x, double dx, double v[]) {
|
||||
double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
|
||||
sn,ssn,cs,ccs,ds,dss,dc,dcc;
|
||||
#if 0
|
||||
@ -96,7 +96,7 @@ void dubsin(double x, double dx, double v[]) {
|
||||
/*(x+dx) between 0 and PI/4 */
|
||||
/**********************************************************************/
|
||||
|
||||
void dubcos(double x, double dx, double v[]) {
|
||||
void __dubcos(double x, double dx, double v[]) {
|
||||
double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
|
||||
sn,ssn,cs,ccs,ds,dss,dc,dcc;
|
||||
#if 0
|
||||
@ -159,20 +159,20 @@ void dubcos(double x, double dx, double v[]) {
|
||||
/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
|
||||
/* as Double-Length number and store it in array v */
|
||||
/**********************************************************************/
|
||||
void docos(double x, double dx, double v[]) {
|
||||
void __docos(double x, double dx, double v[]) {
|
||||
double y,yy,p,w[2];
|
||||
if (x>0) {y=x; yy=dx;}
|
||||
else {y=-x; yy=-dx;}
|
||||
if (y<0.5*hp0.x) /* y< PI/4 */
|
||||
{dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
|
||||
{__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
|
||||
else if (y<1.5*hp0.x) { /* y< 3/4 * PI */
|
||||
p=hp0.x-y; /* p = PI/2 - y */
|
||||
yy=hp1.x-yy;
|
||||
y=p+yy;
|
||||
yy=(p-y)+yy;
|
||||
if (y>0) {dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
|
||||
if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
|
||||
/* cos(x) = sin ( 90 - x ) */
|
||||
else {dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
|
||||
else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
|
||||
}
|
||||
}
|
||||
else { /* y>= 3/4 * PI */
|
||||
@ -180,7 +180,7 @@ void docos(double x, double dx, double v[]) {
|
||||
yy=2.0*hp1.x-yy;
|
||||
y=p+yy;
|
||||
yy=(p-y)+yy;
|
||||
dubcos(y,yy,w);
|
||||
__dubcos(y,yy,w);
|
||||
v[0]=-w[0];
|
||||
v[1]=-w[1];
|
||||
}
|
||||
|
||||
@ -312,6 +312,8 @@ double __ieee754_asin(double x){
|
||||
} /* else if (k < 0x3ff00000) */
|
||||
/*---------------------------- |x|>=1 -------------------------------*/
|
||||
else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x;
|
||||
else
|
||||
if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
|
||||
else {
|
||||
u.i[HIGH_HALF]=0x7ff00000;
|
||||
v.i[HIGH_HALF]=0x7ff00000;
|
||||
@ -621,6 +623,8 @@ double __ieee754_acos(double x)
|
||||
/*---------------------------- |x|>=1 -----------------------*/
|
||||
else
|
||||
if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x;
|
||||
else
|
||||
if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
|
||||
else {
|
||||
u.i[HIGH_HALF]=0x7ff00000;
|
||||
v.i[HIGH_HALF]=0x7ff00000;
|
||||
|
||||
@ -374,9 +374,9 @@ static double normalized(double ax,double ay,double y, double z)
|
||||
{ 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);
|
||||
sub(&mpz,&mperr,&mpz2,p); __mp_dbl(&mpz2,&z,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);
|
||||
}
|
||||
/* Fix the sign and return after stage 1 or stage 2 */
|
||||
@ -392,10 +392,10 @@ static double atan2Mp(double x,double y,const int pr[])
|
||||
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);
|
||||
add(&mpz,&mperr,&mpz1,p); sub(&mpz,&mperr,&mpz2,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;
|
||||
}
|
||||
|
||||
@ -175,7 +175,7 @@ static double zero= 0.00000000000000000000e+00;
|
||||
GET_HIGH_WORD(ix,x);
|
||||
ix &= 0x7fffffff;
|
||||
|
||||
if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0);
|
||||
if(ix<0x3fd00000) return sin(pi*x);
|
||||
y = -x; /* x is assume negative */
|
||||
|
||||
/*
|
||||
|
||||
@ -189,10 +189,10 @@ double __ieee754_log(double x) {
|
||||
|
||||
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);
|
||||
add(&mpy,&mperr,&mpy1,p); sub(&mpy,&mperr,&mpy2,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;
|
||||
}
|
||||
|
||||
@ -45,7 +45,7 @@
|
||||
double __exp1(double x, double xx, double error);
|
||||
static double log1(double x, double *delta, double *error);
|
||||
static double log2(double x, double *delta, double *error);
|
||||
double slowpow(double x, double y,double z);
|
||||
double __slowpow(double x, double y,double z);
|
||||
static double power1(double x, double y);
|
||||
static int checkint(double x);
|
||||
|
||||
@ -69,8 +69,8 @@ double __ieee754_pow(double x, double y) {
|
||||
if (((qx==0x7ff00000)&&(u.i[LOW_HALF]!=0))||(qx>0x7ff00000)) return NaNQ.x;
|
||||
if (y == 1.0) return x;
|
||||
if (y == 2.0) return x*x;
|
||||
if (y == -1.0) return (x!=0)?1.0/x:NaNQ.x;
|
||||
if (y == 0) return ((x>0)&&(qx<0x7ff00000))?1.0:NaNQ.x;
|
||||
if (y == -1.0) return 1.0/x;
|
||||
if (y == 0) return 1.0;
|
||||
}
|
||||
/* else */
|
||||
if(((u.i[HIGH_HALF]>0 && u.i[HIGH_HALF]<0x7ff00000)|| /* x>0 and not x->0 */
|
||||
@ -94,16 +94,37 @@ double __ieee754_pow(double x, double y) {
|
||||
}
|
||||
|
||||
if (x == 0) {
|
||||
if (ABS(y) > 1.0e20) return (y>0)?0:NaNQ.x;
|
||||
if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
|
||||
|| (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)
|
||||
return y;
|
||||
if (ABS(y) > 1.0e20) return (y>0)?0:INF.x;
|
||||
k = checkint(y);
|
||||
if (k == 0 || y < 0) return NaNQ.x;
|
||||
else return (k==1)?0:x; /* return 0 */
|
||||
if (k == -1)
|
||||
return y < 0 ? 1.0/x : x;
|
||||
else
|
||||
return y < 0 ? 1.0/ABS(x) : 0.0; /* return 0 */
|
||||
}
|
||||
/* if x<0 */
|
||||
if (u.i[HIGH_HALF] < 0) {
|
||||
k = checkint(y);
|
||||
if (k==0) return NaNQ.x; /* y not integer and x<0 */
|
||||
return (k==1)?upow(-x,y):-upow(-x,y); /* if y even or odd */
|
||||
if (k==0) {
|
||||
if ((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] == 0) {
|
||||
if (x == -1.0) return 1.0;
|
||||
else if (x > -1.0) return v.i[HIGH_HALF] < 0 ? INF.x : 0.0;
|
||||
else return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x;
|
||||
}
|
||||
else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0)
|
||||
return y < 0 ? 0.0 : INF.x;
|
||||
return NaNQ.x; /* y not integer and x<0 */
|
||||
}
|
||||
else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0)
|
||||
{
|
||||
if (k < 0)
|
||||
return y < 0 ? nZERO.x : nINF.x;
|
||||
else
|
||||
return y < 0 ? 0.0 : INF.x;
|
||||
}
|
||||
return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */
|
||||
}
|
||||
/* x>0 */
|
||||
qx = u.i[HIGH_HALF]&0x7fffffff; /* no sign */
|
||||
@ -111,7 +132,8 @@ double __ieee754_pow(double x, double y) {
|
||||
|
||||
if (qx > 0x7ff00000 || (qx == 0x7ff00000 && u.i[LOW_HALF] != 0)) return NaNQ.x;
|
||||
/* if 0<x<2^-0x7fe */
|
||||
if (qy > 0x7ff00000 || (qy == 0x7ff00000 && v.i[LOW_HALF] != 0)) return NaNQ.x;
|
||||
if (qy > 0x7ff00000 || (qy == 0x7ff00000 && v.i[LOW_HALF] != 0))
|
||||
return x == 1.0 ? 1.0 : NaNQ.x;
|
||||
/* if y<2^-0x7fe */
|
||||
|
||||
if (qx == 0x7ff00000) /* x= 2^-0x3ff */
|
||||
@ -124,7 +146,7 @@ double __ieee754_pow(double x, double y) {
|
||||
if (y<0) return (x<1.0)?INF.x:0;
|
||||
}
|
||||
|
||||
if (x == 1.0) return NaNQ.x;
|
||||
if (x == 1.0) return 1.0;
|
||||
if (y>0) return (x>1.0)?INF.x:0;
|
||||
if (y<0) return (x<1.0)?INF.x:0;
|
||||
return 0; /* unreachable, to make the compiler happy */
|
||||
@ -148,7 +170,7 @@ static double power1(double x, double y) {
|
||||
a2 = (a-a1)+aa;
|
||||
error = error*ABS(y);
|
||||
t = __exp1(a1,a2,1.9e16*error);
|
||||
return (t >= 0)?t:slowpow(x,y,z);
|
||||
return (t >= 0)?t:__slowpow(x,y,z);
|
||||
}
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
@ -103,11 +103,13 @@ double __ieee754_remainder(double x, double y)
|
||||
else {
|
||||
if (kx<0x7ff00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
|
||||
y=ABS(y)*t128.x;
|
||||
z=uremainder(x,y)*t128.x;
|
||||
z=uremainder(z,y)*tm128.x;
|
||||
z=__ieee754_remainder(x,y)*t128.x;
|
||||
z=__ieee754_remainder(z,y)*tm128.x;
|
||||
return z;
|
||||
}
|
||||
else { /* if x is too big */
|
||||
if (kx == 0x7ff00000 && u.i[LOW_HALF] == 0 && y == 1.0)
|
||||
return x / x;
|
||||
if (kx>=0x7ff00000||(ky==0&&t.i[LOW_HALF]==0)||ky>0x7ff00000||
|
||||
(ky==0x7ff00000&&t.i[LOW_HALF]!=0))
|
||||
return (u.i[HIGH_HALF]&0x80000000)?nNAN.x:NAN.x;
|
||||
|
||||
@ -5,9 +5,9 @@
|
||||
*
|
||||
* 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
|
||||
* the Free Software Foundation; either version 2 of the License, or
|
||||
* the Free Software Foundation; either version 2 of the License, 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
|
||||
@ -15,12 +15,12 @@
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License
|
||||
* along with this program; if not, write to the Free Software
|
||||
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
|
||||
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
|
||||
*/
|
||||
/************************************************************************/
|
||||
/* */
|
||||
/* MODULE_NAME:halfulp.c */
|
||||
/* */
|
||||
/* MODULE_NAME:halfulp.c */
|
||||
/* */
|
||||
/* FUNCTIONS:halfulp */
|
||||
/* FILES NEEDED: mydefs.h dla.h endian.h */
|
||||
/* uroot.c */
|
||||
@ -35,11 +35,11 @@
|
||||
/*3. if x can be represented by x=2**n for some integer n. */
|
||||
/************************************************************************/
|
||||
|
||||
#include "endian.h"
|
||||
#include "endian.h"
|
||||
#include "mydefs.h"
|
||||
#include "dla.h"
|
||||
|
||||
double usqrt(double x);
|
||||
double __ieee754_sqrt(double x);
|
||||
|
||||
int4 tab54[32] = {
|
||||
262143, 11585, 1782, 511, 210, 107, 63, 42,
|
||||
@ -48,17 +48,17 @@ int4 tab54[32] = {
|
||||
3, 3, 3, 3, 3, 3, 3, 3 };
|
||||
|
||||
|
||||
double halfulp(double x, double y)
|
||||
double __halfulp(double x, double y)
|
||||
{
|
||||
mynumber v;
|
||||
double z,u,uu,j1,j2,j3,j4,j5;
|
||||
int4 k,l,m,n;
|
||||
if (y <= 0) { /*if power is negative or zero */
|
||||
v.x = y;
|
||||
if (v.i[LOW_HALF] != 0) return -10.0;
|
||||
if (v.i[LOW_HALF] != 0) return -10.0;
|
||||
v.x = x;
|
||||
if (v.i[LOW_HALF] != 0) return -10.0;
|
||||
if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */
|
||||
if (v.i[LOW_HALF] != 0) return -10.0;
|
||||
if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */
|
||||
k = ((v.i[HIGH_HALF]&0x7fffffff)>>20)-1023; /* find this n */
|
||||
z = (double) k;
|
||||
return (z*y == -1075.0)?0: -10.0;
|
||||
@ -66,19 +66,19 @@ double halfulp(double x, double y)
|
||||
/* if y > 0 */
|
||||
v.x = y;
|
||||
if (v.i[LOW_HALF] != 0) return -10.0;
|
||||
|
||||
|
||||
v.x=x;
|
||||
/* case where x = 2**n for some integer n */
|
||||
/* case where x = 2**n for some integer n */
|
||||
if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) {
|
||||
k=(v.i[HIGH_HALF]>>20)-1023;
|
||||
return (((double) k)*y == -1075.0)?0:-10.0;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
v.x = y;
|
||||
k = v.i[HIGH_HALF];
|
||||
m = k<<12;
|
||||
l = 0;
|
||||
while (m)
|
||||
while (m)
|
||||
{m = m<<1; l++; }
|
||||
n = (k&0x000fffff)|0x00100000;
|
||||
n = n>>(20-l); /* n is the odd integer of y */
|
||||
@ -88,19 +88,19 @@ double halfulp(double x, double y)
|
||||
if (n > 34) return -10.0;
|
||||
k = -k;
|
||||
if (k>5) return -10.0;
|
||||
|
||||
|
||||
/* now treat x */
|
||||
while (k>0) {
|
||||
z = usqrt(x);
|
||||
z = __ieee754_sqrt(x);
|
||||
EMULV(z,z,u,uu,j1,j2,j3,j4,j5);
|
||||
if (((u-x)+uu) != 0) break;
|
||||
x = z;
|
||||
k--;
|
||||
}
|
||||
if (k) return -10.0;
|
||||
|
||||
if (k) return -10.0;
|
||||
|
||||
/* it is impossible that n == 2, so the mantissa of x must be short */
|
||||
|
||||
|
||||
v.x = x;
|
||||
if (v.i[LOW_HALF]) return -10.0;
|
||||
k = v.i[HIGH_HALF];
|
||||
@ -109,16 +109,14 @@ double halfulp(double x, double y)
|
||||
while (m) {m = m<<1; l++; }
|
||||
m = (k&0x000fffff)|0x00100000;
|
||||
m = m>>(20-l); /* m is the odd integer of x */
|
||||
|
||||
|
||||
/* now check whether the length of m**n is at most 54 bits */
|
||||
|
||||
|
||||
if (m > tab54[n-3]) return -10.0;
|
||||
|
||||
|
||||
/* yes, it is - now compute x**n by simple multiplications */
|
||||
|
||||
|
||||
u = x;
|
||||
for (k=1;k<n;k++) u = u*x;
|
||||
return u;
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -62,7 +62,7 @@ static int mcr(const mp_no *x, const mp_no *y, int p) {
|
||||
|
||||
|
||||
/* acr() compares the absolute values of two multiple precision numbers */
|
||||
int acr(const mp_no *x, const mp_no *y, int p) {
|
||||
int __acr(const mp_no *x, const mp_no *y, int p) {
|
||||
int i;
|
||||
|
||||
if (X[0] == ZERO) {
|
||||
@ -81,20 +81,20 @@ int acr(const mp_no *x, const mp_no *y, int p) {
|
||||
|
||||
|
||||
/* cr90 compares the values of two multiple precision numbers */
|
||||
int cr(const mp_no *x, const mp_no *y, int p) {
|
||||
int __cr(const mp_no *x, const mp_no *y, int p) {
|
||||
int i;
|
||||
|
||||
if (X[0] > Y[0]) i= 1;
|
||||
else if (X[0] < Y[0]) i=-1;
|
||||
else if (X[0] < ZERO ) i= acr(y,x,p);
|
||||
else i= acr(x,y,p);
|
||||
else if (X[0] < ZERO ) i= __acr(y,x,p);
|
||||
else i= __acr(x,y,p);
|
||||
|
||||
return i;
|
||||
}
|
||||
|
||||
|
||||
/* Copy a multiple precision number. Set *y=*x. x=y is permissible. */
|
||||
void cpy(const mp_no *x, mp_no *y, int p) {
|
||||
void __cpy(const mp_no *x, mp_no *y, int p) {
|
||||
int i;
|
||||
|
||||
EY = EX;
|
||||
@ -110,7 +110,7 @@ void cpy(const mp_no *x, mp_no *y, int p) {
|
||||
/* n<m, the digits of x beyond the n'th are ignored. */
|
||||
/* x=y is permissible. */
|
||||
|
||||
void cpymn(const mp_no *x, int m, mp_no *y, int n) {
|
||||
void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
|
||||
|
||||
int i,k;
|
||||
|
||||
@ -246,7 +246,7 @@ void __mp_dbl(const mp_no *x, double *y, int p) {
|
||||
/* 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;
|
||||
@ -286,7 +286,7 @@ static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
i=p; j=p+ EY - EX; k=p+1;
|
||||
|
||||
if (j<1)
|
||||
{cpy(x,z,p); return; }
|
||||
{__cpy(x,z,p); return; }
|
||||
else Z[k] = ZERO;
|
||||
|
||||
for (; j>0; i--,j--) {
|
||||
@ -330,7 +330,7 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
Z[k] = Z[k+1] = ZERO; }
|
||||
else {
|
||||
j= EX - EY;
|
||||
if (j > p) {cpy(x,z,p); return; }
|
||||
if (j > p) {__cpy(x,z,p); return; }
|
||||
else {
|
||||
i=p; j=p+1-j; k=p;
|
||||
if (Y[j] > ZERO) {
|
||||
@ -375,19 +375,19 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
/* but not x&z or y&z. One guard digit is used. The error is less than */
|
||||
/* one ulp. *x & *y are left unchanged. */
|
||||
|
||||
void add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
void __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] == ZERO) {__cpy(y,z,p); return; }
|
||||
else if (Y[0] == ZERO) {__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]; }
|
||||
if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; }
|
||||
else {add_magnitudes(y,x,z,p); Z[0] = Y[0]; }
|
||||
}
|
||||
else {
|
||||
if ((n=acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; }
|
||||
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;
|
||||
}
|
||||
@ -399,19 +399,19 @@ void add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
/* overlap but not x&z or y&z. One guard digit is used. The error is */
|
||||
/* less than one ulp. *x & *y are left unchanged. */
|
||||
|
||||
void sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
void __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] == ZERO) {__cpy(y,z,p); Z[0] = -Z[0]; return; }
|
||||
else if (Y[0] == ZERO) {__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]; }
|
||||
if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; }
|
||||
else {add_magnitudes(y,x,z,p); Z[0] = -Y[0]; }
|
||||
}
|
||||
else {
|
||||
if ((n=acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; }
|
||||
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;
|
||||
}
|
||||
@ -424,7 +424,7 @@ void sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
/* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp. */
|
||||
/* *x & *y are left unchanged. */
|
||||
|
||||
void mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
|
||||
int i, i1, i2, j, k, k2;
|
||||
double u;
|
||||
@ -464,7 +464,7 @@ void mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
/* 2.001*r**(1-p) for p>3. */
|
||||
/* *x=0 is not permissible. *x is left unchanged. */
|
||||
|
||||
void inv(const mp_no *x, mp_no *y, int p) {
|
||||
void __inv(const mp_no *x, mp_no *y, int p) {
|
||||
int i;
|
||||
#if 0
|
||||
int l;
|
||||
@ -478,14 +478,14 @@ void inv(const mp_no *x, mp_no *y, int p) {
|
||||
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
|
||||
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;
|
||||
__cpy(x,&z,p); z.e=0; __mp_dbl(&z,&t,p);
|
||||
t=ONE/t; __dbl_mp(t,y,p); EY -= EX;
|
||||
|
||||
for (i=0; i<np1[p]; i++) {
|
||||
cpy(y,&w,p);
|
||||
mul(x,&w,y,p);
|
||||
sub(&mptwo,y,&z,p);
|
||||
mul(&w,&z,y,p);
|
||||
__cpy(y,&w,p);
|
||||
__mul(x,&w,y,p);
|
||||
__sub(&mptwo,y,&z,p);
|
||||
__mul(&w,&z,y,p);
|
||||
}
|
||||
return;
|
||||
}
|
||||
@ -496,11 +496,11 @@ void inv(const mp_no *x, mp_no *y, int p) {
|
||||
/* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3 */
|
||||
/* and 3.001*r**(1-p) for p>3. *y=0 is not permissible. */
|
||||
|
||||
void dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
|
||||
|
||||
mp_no w;
|
||||
|
||||
if (X[0] == ZERO) Z[0] = ZERO;
|
||||
else {inv(y,&w,p); mul(x,&w,z,p);}
|
||||
else {__inv(y,&w,p); __mul(x,&w,z,p);}
|
||||
return;
|
||||
}
|
||||
|
||||
@ -65,14 +65,14 @@ typedef union { int i[2]; double d; } number;
|
||||
#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);
|
||||
|
||||
@ -33,9 +33,9 @@
|
||||
|
||||
#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;
|
||||
@ -66,36 +66,36 @@ void mpatan(mp_no *x, mp_no *y, int p) {
|
||||
mptwo.d[1] = TWO;
|
||||
|
||||
/* Reduce x m times */
|
||||
mul(x,x,&mpsm,p);
|
||||
if (m==0) cpy(x,&mps,p);
|
||||
__mul(x,x,&mpsm,p);
|
||||
if (m==0) __cpy(x,&mps,p);
|
||||
else {
|
||||
for (i=0; i<m; i++) {
|
||||
add(&mpone,&mpsm,&mpt1,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);
|
||||
__add(&mpone,&mpsm,&mpt1,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) */
|
||||
n=np[p]; mptwoim1.d[1] = twonm1[p].d;
|
||||
dvd(&mpsm,&mptwoim1,&mpt,p);
|
||||
__dvd(&mpsm,&mptwoim1,&mpt,p);
|
||||
for (i=n-1; i>1; i--) {
|
||||
mptwoim1.d[1] -= TWO;
|
||||
dvd(&mpsm,&mptwoim1,&mpt1,p);
|
||||
mul(&mpsm,&mpt,&mpt2,p);
|
||||
sub(&mpt1,&mpt2,&mpt,p);
|
||||
__dvd(&mpsm,&mptwoim1,&mpt1,p);
|
||||
__mul(&mpsm,&mpt,&mpt2,p);
|
||||
__sub(&mpt1,&mpt2,&mpt,p);
|
||||
}
|
||||
mul(&mps,&mpt,&mpt1,p);
|
||||
sub(&mps,&mpt1,&mpt,p);
|
||||
__mul(&mps,&mpt,&mpt1,p);
|
||||
__sub(&mps,&mpt1,&mpt,p);
|
||||
|
||||
/* Compute Atan(x) */
|
||||
mptwoim1.d[1] = twom[m].d;
|
||||
mul(&mptwoim1,&mpt,y,p);
|
||||
__mul(&mptwoim1,&mpt,y,p);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
@ -37,12 +37,12 @@
|
||||
|
||||
#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;
|
||||
|
||||
@ -54,15 +54,15 @@ void mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) {
|
||||
|
||||
if (X[0] <= ZERO) {
|
||||
mpone.e = 1; mpone.d[0] = mpone.d[1] = ONE;
|
||||
dvd(x,y,&mpt1,p); mul(&mpt1,&mpt1,&mpt2,p);
|
||||
__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(&mpt1,&mpt2,&mpt3,p); mpt3.d[0]=Y[0];
|
||||
mpatan(&mpt3,&mpt1,p); add(&mpt1,&mpt1,z,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);
|
||||
}
|
||||
else
|
||||
{ dvd(y,x,&mpt1,p);
|
||||
mpatan(&mpt1,z,p);
|
||||
{ __dvd(y,x,&mpt1,p);
|
||||
__mpatan(&mpt1,z,p);
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
@ -35,7 +35,7 @@
|
||||
|
||||
/* Multi-Precision exponential function subroutine (for p >= 4, */
|
||||
/* 2**(-55) <= abs(x) <= 1024). */
|
||||
void mpexp(mp_no *x, mp_no *y, int p) {
|
||||
void __mpexp(mp_no *x, mp_no *y, int p) {
|
||||
|
||||
int i,j,k,m,m1,m2,n;
|
||||
double a,b;
|
||||
@ -75,30 +75,30 @@ void mpexp(mp_no *x, mp_no *y, int p) {
|
||||
}
|
||||
|
||||
/* Compute s=x*2**(-m). Put result in mps */
|
||||
dbl_mp(a,&mpt1,p);
|
||||
mul(x,&mpt1,&mps,p);
|
||||
__dbl_mp(a,&mpt1,p);
|
||||
__mul(x,&mpt1,&mps,p);
|
||||
|
||||
/* Evaluate the polynomial. Put result in mpt2 */
|
||||
mpone.e=1; mpone.d[0]=ONE; mpone.d[1]=ONE;
|
||||
mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=nn[n].d;
|
||||
dvd(&mps,&mpk,&mpt1,p);
|
||||
add(&mpone,&mpt1,&mpak,p);
|
||||
__dvd(&mps,&mpk,&mpt1,p);
|
||||
__add(&mpone,&mpt1,&mpak,p);
|
||||
for (k=n-1; k>1; k--) {
|
||||
mul(&mps,&mpak,&mpt1,p);
|
||||
__mul(&mps,&mpak,&mpt1,p);
|
||||
mpk.d[1]=nn[k].d;
|
||||
dvd(&mpt1,&mpk,&mpt2,p);
|
||||
add(&mpone,&mpt2,&mpak,p);
|
||||
__dvd(&mpt1,&mpk,&mpt2,p);
|
||||
__add(&mpone,&mpt2,&mpak,p);
|
||||
}
|
||||
mul(&mps,&mpak,&mpt1,p);
|
||||
add(&mpone,&mpt1,&mpt2,p);
|
||||
__mul(&mps,&mpak,&mpt1,p);
|
||||
__add(&mpone,&mpt1,&mpt2,p);
|
||||
|
||||
/* Raise polynomial value to the power of 2**m. Put result in y */
|
||||
for (k=0,j=0; k<m; ) {
|
||||
mul(&mpt2,&mpt2,&mpt1,p); k++;
|
||||
__mul(&mpt2,&mpt2,&mpt1,p); k++;
|
||||
if (k==m) { j=1; break; }
|
||||
mul(&mpt1,&mpt1,&mpt2,p); k++;
|
||||
__mul(&mpt1,&mpt1,&mpt2,p); k++;
|
||||
}
|
||||
if (j) cpy(&mpt1,y,p);
|
||||
else cpy(&mpt2,y,p);
|
||||
if (j) __cpy(&mpt1,y,p);
|
||||
else __cpy(&mpt2,y,p);
|
||||
return;
|
||||
}
|
||||
|
||||
@ -37,9 +37,9 @@
|
||||
#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
|
||||
@ -58,14 +58,14 @@ void mplog(mp_no *x, mp_no *y, int p) {
|
||||
|
||||
/* Perform m newton iterations to solve for y: exp(y)-x=0. */
|
||||
/* The iterations formula is: y(n+1)=y(n)+(x*exp(-y(n))-1). */
|
||||
cpy(y,&mpt1,p);
|
||||
__cpy(y,&mpt1,p);
|
||||
for (i=0; i<m; i++) {
|
||||
mpt1.d[0]=-mpt1.d[0];
|
||||
mpexp(&mpt1,&mpt2,p);
|
||||
mul(x,&mpt2,&mpt1,p);
|
||||
sub(&mpt1,&mpone,&mpt2,p);
|
||||
add(y,&mpt2,&mpt1,p);
|
||||
cpy(&mpt1,y,p);
|
||||
__mpexp(&mpt1,&mpt2,p);
|
||||
__mul(x,&mpt2,&mpt1,p);
|
||||
__sub(&mpt1,&mpone,&mpt2,p);
|
||||
__add(y,&mpt2,&mpt1,p);
|
||||
__cpy(&mpt1,y,p);
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
@ -42,7 +42,7 @@
|
||||
|
||||
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;
|
||||
@ -60,19 +60,19 @@ void mpsqrt(mp_no *x, mp_no *y, int p) {
|
||||
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;
|
||||
|
||||
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);
|
||||
mul(&mpxn,&mphalf,&mpz,p);
|
||||
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);
|
||||
__mul(&mpxn,&mphalf,&mpz,p);
|
||||
|
||||
m=mp[p];
|
||||
for (i=0; i<m; i++) {
|
||||
mul(&mpu,&mpu,&mpt1,p);
|
||||
mul(&mpt1,&mpz,&mpt2,p);
|
||||
sub(&mp3halfs,&mpt2,&mpt1,p);
|
||||
mul(&mpu,&mpt1,&mpt2,p);
|
||||
cpy(&mpt2,&mpu,p);
|
||||
__mul(&mpu,&mpu,&mpt1,p);
|
||||
__mul(&mpt1,&mpz,&mpt2,p);
|
||||
__sub(&mp3halfs,&mpt2,&mpt1,p);
|
||||
__mul(&mpu,&mpt1,&mpt2,p);
|
||||
__cpy(&mpt2,&mpu,p);
|
||||
}
|
||||
mul(&mpxn,&mpu,y,p); EY += ey;
|
||||
__mul(&mpxn,&mpu,y,p); EY += ey;
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
@ -37,23 +37,23 @@
|
||||
#include "endian.h"
|
||||
#include "mpa.h"
|
||||
|
||||
int mpranred(double, mp_no *, int);
|
||||
int __mpranred(double, mp_no *, int);
|
||||
void __c32(mp_no *, mp_no *, mp_no *, int);
|
||||
|
||||
void mptan(double x, mp_no *mpy, int p) {
|
||||
void __mptan(double x, mp_no *mpy, int p) {
|
||||
|
||||
static const double MONE = -1.0;
|
||||
|
||||
int n;
|
||||
mp_no mpw, mpc, mps;
|
||||
|
||||
n = mpranred(x, &mpw, p) & 0x00000001; /* negative or positive result */
|
||||
n = __mpranred(x, &mpw, p) & 0x00000001; /* negative or positive result */
|
||||
__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);
|
||||
{ __dvd(&mpc,&mps,mpy,p);
|
||||
mpy->d[0] *= MONE;
|
||||
} /* tan is negative in this area */
|
||||
else dvd(&mps,&mpc,mpy,p);
|
||||
else __dvd(&mps,&mpc,mpy,p);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
@ -214,9 +214,9 @@ static double atanMp(double x,const int pr[]){
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
@ -60,8 +60,8 @@ static const double
|
||||
cs4 = -4.16666666666664434524222570944589E-02,
|
||||
cs6 = 1.38888874007937613028114285595617E-03;
|
||||
|
||||
void dubsin(double x, double dx, double w[]);
|
||||
void docos(double x, double dx, double w[]);
|
||||
void __dubsin(double x, double dx, double w[]);
|
||||
void __docos(double x, double dx, double w[]);
|
||||
double __mpsin(double x, double dx);
|
||||
double __mpcos(double x, double dx);
|
||||
double __mpsin1(double x);
|
||||
@ -75,7 +75,7 @@ static double sloww2(double x, double dx, double orig, int n);
|
||||
static double bsloww(double x, double dx, double orig, int n);
|
||||
static double bsloww1(double x, double dx, double orig, int n);
|
||||
static double bsloww2(double x, double dx, double orig, int n);
|
||||
int branred(double x, double *a, double *aa);
|
||||
int __branred(double x, double *a, double *aa);
|
||||
static double cslow2(double x);
|
||||
static double csloww(double x, double dx, double orig);
|
||||
static double csloww1(double x, double dx, double orig);
|
||||
@ -84,7 +84,7 @@ static double csloww2(double x, double dx, double orig, int n);
|
||||
/* An ultimate sin routine. Given an IEEE double machine number x */
|
||||
/* it computes the correctly rounded (to nearest) value of sin(x) */
|
||||
/*******************************************************************/
|
||||
double sin(double x){
|
||||
double __sin(double x){
|
||||
double xx,res,t,cor,y,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2;
|
||||
#if 0
|
||||
double w[2];
|
||||
@ -307,7 +307,7 @@ double sin(double x){
|
||||
/* -----------------281474976710656 <|x| <2^1024----------------------------*/
|
||||
else if (k < 0x7ff00000) {
|
||||
|
||||
n = branred(x,&a,&da);
|
||||
n = __branred(x,&a,&da);
|
||||
switch (n) {
|
||||
case 0:
|
||||
if (a*a < 0.01588) return bsloww(a,da,x,n);
|
||||
@ -327,7 +327,7 @@ double sin(double x){
|
||||
} /* else if (k < 0x7ff00000 ) */
|
||||
|
||||
/*--------------------- |x| > 2^1024 ----------------------------------*/
|
||||
else return NAN.x;
|
||||
else return x / x;
|
||||
return 0; /* unreachable */
|
||||
}
|
||||
|
||||
@ -337,7 +337,7 @@ double sin(double x){
|
||||
/* it computes the correctly rounded (to nearest) value of cos(x) */
|
||||
/*******************************************************************/
|
||||
|
||||
double cos(double x)
|
||||
double __cos(double x)
|
||||
{
|
||||
double y,xx,res,t,cor,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2;
|
||||
mynumber u,v;
|
||||
@ -548,7 +548,7 @@ double cos(double x)
|
||||
|
||||
else if (k < 0x7ff00000) {/* 281474976710656 <|x| <2^1024 */
|
||||
|
||||
n = branred(x,&a,&da);
|
||||
n = __branred(x,&a,&da);
|
||||
switch (n) {
|
||||
case 1:
|
||||
if (a*a < 0.01588) return bsloww(-a,-da,x,n);
|
||||
@ -570,7 +570,7 @@ double cos(double x)
|
||||
|
||||
|
||||
|
||||
else return NAN.x; /* |x| > 2^1024 */
|
||||
else return x / x; /* |x| > 2^1024 */
|
||||
return 0;
|
||||
|
||||
}
|
||||
@ -594,7 +594,7 @@ static const double th2_36 = 206158430208.0; /* 1.5*2**37 */
|
||||
cor = (r-res)+t;
|
||||
if (res == res + 1.0007*cor) return res;
|
||||
else {
|
||||
dubsin(ABS(x),0,w);
|
||||
__dubsin(ABS(x),0,w);
|
||||
if (w[0] == w[0]+1.000000001*w[1]) return (x>0)?w[0]:-w[0];
|
||||
else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
|
||||
}
|
||||
@ -631,7 +631,7 @@ static double slow1(double x) {
|
||||
cor=(y-res)+cor;
|
||||
if (res == res+1.0005*cor) return (x>0)?res:-res;
|
||||
else {
|
||||
dubsin(ABS(x),0,w);
|
||||
__dubsin(ABS(x),0,w);
|
||||
if (w[0] == w[0]+1.000000005*w[1]) return (x>0)?w[0]:-w[0];
|
||||
else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
|
||||
}
|
||||
@ -679,7 +679,7 @@ static double slow2(double x) {
|
||||
y=ABS(x)-hp0.x;
|
||||
y1=y-hp1.x;
|
||||
y2=(y-y1)-hp1.x;
|
||||
docos(y1,y2,w);
|
||||
__docos(y1,y2,w);
|
||||
if (w[0] == w[0]+1.000000005*w[1]) return (x>0)?w[0]:-w[0];
|
||||
else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
|
||||
}
|
||||
@ -709,7 +709,7 @@ static double sloww(double x,double dx, double orig) {
|
||||
cor = (cor>0)? 1.0005*cor+ABS(orig)*3.1e-30 : 1.0005*cor-ABS(orig)*3.1e-30;
|
||||
if (res == res + cor) return res;
|
||||
else {
|
||||
(x>0)? dubsin(x,dx,w) : dubsin(-x,-dx,w);
|
||||
(x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w);
|
||||
cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-30 : 1.000000001*w[1] - ABS(orig)*1.1e-30;
|
||||
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
|
||||
else {
|
||||
@ -725,7 +725,7 @@ static double sloww(double x,double dx, double orig) {
|
||||
a = t - y;
|
||||
da = ((t-a)-y)+da;
|
||||
if (n&2) {a=-a; da=-da;}
|
||||
(a>0)? dubsin(a,da,w) : dubsin(-a,-da,w);
|
||||
(a>0)? __dubsin(a,da,w) : __dubsin(-a,-da,w);
|
||||
cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-40 : 1.000000001*w[1] - ABS(orig)*1.1e-40;
|
||||
if (w[0] == w[0]+cor) return (a>0)?w[0]:-w[0];
|
||||
else return __mpsin1(orig);
|
||||
@ -768,7 +768,7 @@ static double sloww1(double x, double dx, double orig) {
|
||||
cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
|
||||
if (res == res + cor) return (x>0)?res:-res;
|
||||
else {
|
||||
dubsin(ABS(x),dx,w);
|
||||
__dubsin(ABS(x),dx,w);
|
||||
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
|
||||
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
|
||||
else return __mpsin1(orig);
|
||||
@ -811,7 +811,7 @@ static double sloww2(double x, double dx, double orig, int n) {
|
||||
cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
|
||||
if (res == res + cor) return (n&2)?-res:res;
|
||||
else {
|
||||
docos(ABS(x),dx,w);
|
||||
__docos(ABS(x),dx,w);
|
||||
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
|
||||
if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0];
|
||||
else return __mpsin1(orig);
|
||||
@ -844,7 +844,7 @@ static double bsloww(double x,double dx, double orig,int n) {
|
||||
cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24;
|
||||
if (res == res + cor) return res;
|
||||
else {
|
||||
(x>0)? dubsin(x,dx,w) : dubsin(-x,-dx,w);
|
||||
(x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w);
|
||||
cor = (w[1]>0)? 1.000000001*w[1] + 1.1e-24 : 1.000000001*w[1] - 1.1e-24;
|
||||
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
|
||||
else return (n&1)?__mpcos1(orig):__mpsin1(orig);
|
||||
@ -887,7 +887,7 @@ mynumber u;
|
||||
cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24;
|
||||
if (res == res + cor) return (x>0)?res:-res;
|
||||
else {
|
||||
dubsin(ABS(x),dx,w);
|
||||
__dubsin(ABS(x),dx,w);
|
||||
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24: 1.000000005*w[1]-1.1e-24;
|
||||
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
|
||||
else return (n&1)?__mpcos1(orig):__mpsin1(orig);
|
||||
@ -931,7 +931,7 @@ mynumber u;
|
||||
cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24;
|
||||
if (res == res + cor) return (n&2)?-res:res;
|
||||
else {
|
||||
docos(ABS(x),dx,w);
|
||||
__docos(ABS(x),dx,w);
|
||||
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24 : 1.000000005*w[1]-1.1e-24;
|
||||
if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0];
|
||||
else return (n&1)?__mpsin1(orig):__mpcos1(orig);
|
||||
@ -972,7 +972,7 @@ static double cslow2(double x) {
|
||||
return res;
|
||||
else {
|
||||
y=ABS(x);
|
||||
docos(y,0,w);
|
||||
__docos(y,0,w);
|
||||
if (w[0] == w[0]+1.000000005*w[1]) return w[0];
|
||||
else return __mpcos(x,0);
|
||||
}
|
||||
@ -1004,7 +1004,7 @@ static double csloww(double x,double dx, double orig) {
|
||||
cor = (cor>0)? 1.0005*cor+ABS(orig)*3.1e-30 : 1.0005*cor-ABS(orig)*3.1e-30;
|
||||
if (res == res + cor) return res;
|
||||
else {
|
||||
(x>0)? dubsin(x,dx,w) : dubsin(-x,-dx,w);
|
||||
(x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w);
|
||||
cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-30 : 1.000000001*w[1] - ABS(orig)*1.1e-30;
|
||||
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
|
||||
else {
|
||||
@ -1020,7 +1020,7 @@ static double csloww(double x,double dx, double orig) {
|
||||
a = t - y;
|
||||
da = ((t-a)-y)+da;
|
||||
if (n==1) {a=-a; da=-da;}
|
||||
(a>0)? dubsin(a,da,w) : dubsin(-a,-da,w);
|
||||
(a>0)? __dubsin(a,da,w) : __dubsin(-a,-da,w);
|
||||
cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-40 : 1.000000001*w[1] - ABS(orig)*1.1e-40;
|
||||
if (w[0] == w[0]+cor) return (a>0)?w[0]:-w[0];
|
||||
else return __mpcos1(orig);
|
||||
@ -1064,7 +1064,7 @@ static double csloww1(double x, double dx, double orig) {
|
||||
cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
|
||||
if (res == res + cor) return (x>0)?res:-res;
|
||||
else {
|
||||
dubsin(ABS(x),dx,w);
|
||||
__dubsin(ABS(x),dx,w);
|
||||
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
|
||||
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
|
||||
else return __mpcos1(orig);
|
||||
@ -1109,14 +1109,19 @@ static double csloww2(double x, double dx, double orig, int n) {
|
||||
cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
|
||||
if (res == res + cor) return (n)?-res:res;
|
||||
else {
|
||||
docos(ABS(x),dx,w);
|
||||
__docos(ABS(x),dx,w);
|
||||
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
|
||||
if (w[0] == w[0]+cor) return (n)?-w[0]:w[0];
|
||||
else return __mpcos1(orig);
|
||||
}
|
||||
}
|
||||
|
||||
weak_alias (__cos, cos)
|
||||
weak_alias (__sin, sin)
|
||||
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
weak_alias (sin, sinl)
|
||||
weak_alias (cos, cosl)
|
||||
strong_alias (__sin, __sinl)
|
||||
weak_alias (__sin, sinl)
|
||||
strong_alias (__cos, __cosl)
|
||||
weak_alias (__cos, cosl)
|
||||
#endif
|
||||
|
||||
@ -53,8 +53,8 @@ double tan(double x) {
|
||||
mp_no mpy;
|
||||
#endif
|
||||
|
||||
int branred(double, double *, double *);
|
||||
int mpranred(double, mp_no *, int);
|
||||
int __branred(double, double *, double *);
|
||||
int __mpranred(double, mp_no *, int);
|
||||
|
||||
/* x=+-INF, x=NaN */
|
||||
num.d = x; ux = num.i[HIGH_HALF];
|
||||
@ -361,7 +361,7 @@ double tan(double x) {
|
||||
|
||||
/* (---) The case 1e8 < abs(x) < 2**1024 */
|
||||
/* Range reduction by algorithm iii */
|
||||
n = (branred(x,&a,&da)) & 0x00000001;
|
||||
n = (__branred(x,&a,&da)) & 0x00000001;
|
||||
EADD(a,da,t1,t2) a=t1; da=t2;
|
||||
if (a<ZERO) {ya=-a; yya=-da; sy=MONE;}
|
||||
else {ya= a; yya= da; sy= ONE;}
|
||||
@ -384,9 +384,9 @@ double tan(double x) {
|
||||
|
||||
/* Second stage */
|
||||
/* Reduction by algorithm iv */
|
||||
p=10; n = (mpranred(x,&mpa,p)) & 0x00000001;
|
||||
__mp_dbl(&mpa,&a,p); dbl_mp(a,&mpt1,p);
|
||||
sub(&mpa,&mpt1,&mpt2,p); __mp_dbl(&mpt2,&da,p);
|
||||
p=10; n = (__mpranred(x,&mpa,p)) & 0x00000001;
|
||||
__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)
|
||||
c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
|
||||
|
||||
@ -61,17 +61,17 @@ static void ss32(mp_no *x, mp_no *y, int p) {
|
||||
#endif
|
||||
for (i=1;i<=p;i++) mpk.d[i]=0;
|
||||
|
||||
mul(x,x,&x2,p);
|
||||
cpy(&oofac27,&gor,p);
|
||||
cpy(&gor,&sum,p);
|
||||
__mul(x,x,&x2,p);
|
||||
__cpy(&oofac27,&gor,p);
|
||||
__cpy(&gor,&sum,p);
|
||||
for (a=27.0;a>1.0;a-=2.0) {
|
||||
mpk.d[1]=a*(a-1.0);
|
||||
mul(&gor,&mpk,&mpt1,p);
|
||||
cpy(&mpt1,&gor,p);
|
||||
mul(&x2,&sum,&mpt1,p);
|
||||
sub(&gor,&mpt1,&sum,p);
|
||||
__mul(&gor,&mpk,&mpt1,p);
|
||||
__cpy(&mpt1,&gor,p);
|
||||
__mul(&x2,&sum,&mpt1,p);
|
||||
__sub(&gor,&mpt1,&sum,p);
|
||||
}
|
||||
mul(x,&sum,y,p);
|
||||
__mul(x,&sum,y,p);
|
||||
}
|
||||
|
||||
/**********************************************************************/
|
||||
@ -91,18 +91,18 @@ static void cc32(mp_no *x, mp_no *y, int p) {
|
||||
#endif
|
||||
for (i=1;i<=p;i++) mpk.d[i]=0;
|
||||
|
||||
mul(x,x,&x2,p);
|
||||
__mul(x,x,&x2,p);
|
||||
mpk.d[1]=27.0;
|
||||
mul(&oofac27,&mpk,&gor,p);
|
||||
cpy(&gor,&sum,p);
|
||||
__mul(&oofac27,&mpk,&gor,p);
|
||||
__cpy(&gor,&sum,p);
|
||||
for (a=26.0;a>2.0;a-=2.0) {
|
||||
mpk.d[1]=a*(a-1.0);
|
||||
mul(&gor,&mpk,&mpt1,p);
|
||||
cpy(&mpt1,&gor,p);
|
||||
mul(&x2,&sum,&mpt1,p);
|
||||
sub(&gor,&mpt1,&sum,p);
|
||||
__mul(&gor,&mpk,&mpt1,p);
|
||||
__cpy(&mpt1,&gor,p);
|
||||
__mul(&x2,&sum,&mpt1,p);
|
||||
__sub(&gor,&mpt1,&sum,p);
|
||||
}
|
||||
mul(&x2,&sum,y,p);
|
||||
__mul(&x2,&sum,y,p);
|
||||
}
|
||||
|
||||
/***************************************************************************/
|
||||
@ -112,20 +112,20 @@ void __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}};
|
||||
mp_no u,t,t1,t2,c,s;
|
||||
int i;
|
||||
cpy(x,&u,p);
|
||||
__cpy(x,&u,p);
|
||||
u.e=u.e-1;
|
||||
cc32(&u,&c,p);
|
||||
ss32(&u,&s,p);
|
||||
for (i=0;i<24;i++) {
|
||||
mul(&c,&s,&t,p);
|
||||
sub(&s,&t,&t1,p);
|
||||
add(&t1,&t1,&s,p);
|
||||
sub(&mpt,&c,&t1,p);
|
||||
mul(&t1,&c,&t2,p);
|
||||
add(&t2,&t2,&c,p);
|
||||
__mul(&c,&s,&t,p);
|
||||
__sub(&s,&t,&t1,p);
|
||||
__add(&t1,&t1,&s,p);
|
||||
__sub(&mpt,&c,&t1,p);
|
||||
__mul(&t1,&c,&t2,p);
|
||||
__add(&t2,&t2,&c,p);
|
||||
}
|
||||
sub(&one,&c,y,p);
|
||||
cpy(&s,z,p);
|
||||
__sub(&one,&c,y,p);
|
||||
__cpy(&s,z,p);
|
||||
}
|
||||
|
||||
/************************************************************************/
|
||||
@ -137,16 +137,16 @@ double __sin32(double x, double res, double res1) {
|
||||
int p;
|
||||
mp_no a,b,c;
|
||||
p=32;
|
||||
dbl_mp(res,&a,p);
|
||||
dbl_mp(0.5*(res1-res),&b,p);
|
||||
add(&a,&b,&c,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);
|
||||
{ __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 */
|
||||
sub(&b,&c,&a,p);
|
||||
__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;
|
||||
else return (res>res1)?res:res1;
|
||||
@ -161,21 +161,21 @@ double __cos32(double x, double res, double res1) {
|
||||
int p;
|
||||
mp_no a,b,c;
|
||||
p=32;
|
||||
dbl_mp(res,&a,p);
|
||||
dbl_mp(0.5*(res1-res),&b,p);
|
||||
add(&a,&b,&c,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);
|
||||
{ __sub(&pi,&c,&a,p);
|
||||
__c32(&a,&b,&c,p);
|
||||
b.d[0]=-b.d[0];
|
||||
}
|
||||
else if (x>0.8)
|
||||
{ sub(&hp,&c,&a,p);
|
||||
{ __sub(&hp,&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 */
|
||||
sub(&b,&c,&a,p);
|
||||
__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;
|
||||
else return (res<res1)?res:res1;
|
||||
@ -190,10 +190,10 @@ double __mpsin(double x, double dx) {
|
||||
double y;
|
||||
mp_no a,b,c;
|
||||
p=32;
|
||||
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); }
|
||||
__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) */
|
||||
__mp_dbl(&b,&y,p);
|
||||
return y;
|
||||
@ -208,11 +208,11 @@ double __mpcos(double x, double dx) {
|
||||
double y;
|
||||
mp_no a,b,c;
|
||||
p=32;
|
||||
dbl_mp(x,&a,p);
|
||||
dbl_mp(dx,&b,p);
|
||||
add(&a,&b,&c,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);
|
||||
{ __sub(&hp,&c,&b,p);
|
||||
__c32(&b,&a,&c,p);
|
||||
}
|
||||
else __c32(&c,&a,&b,p); /* a = cos(x+dx) */
|
||||
@ -226,7 +226,7 @@ double __mpcos(double x, double dx) {
|
||||
/* n=0,+-1,+-2,.... */
|
||||
/* Return int which indicates in which quarter of circle x is */
|
||||
/******************************************************************/
|
||||
int mpranred(double x, mp_no *y, int p)
|
||||
int __mpranred(double x, mp_no *y, int p)
|
||||
{
|
||||
number v;
|
||||
double t,xn;
|
||||
@ -239,31 +239,31 @@ int mpranred(double x, mp_no *y, int p)
|
||||
xn = t - toint.d;
|
||||
v.d = t;
|
||||
n =v.i[LOW_HALF]&3;
|
||||
dbl_mp(xn,&a,p);
|
||||
mul(&a,&hp,&b,p);
|
||||
dbl_mp(x,&c,p);
|
||||
sub(&c,&b,y,p);
|
||||
__dbl_mp(xn,&a,p);
|
||||
__mul(&a,&hp,&b,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;
|
||||
b.e = -k;
|
||||
b.d[0] = 1.0;
|
||||
for (i=0;i<p;i++) b.d[i+1] = toverp[i+k];
|
||||
mul(&a,&b,&c,p);
|
||||
__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)
|
||||
{ t +=1.0;
|
||||
sub(&c,&one,&b,p);
|
||||
mul(&b,&hp,y,p);
|
||||
__sub(&c,&one,&b,p);
|
||||
__mul(&b,&hp,y,p);
|
||||
}
|
||||
else mul(&c,&hp,y,p);
|
||||
else __mul(&c,&hp,y,p);
|
||||
n = (int) t;
|
||||
if (x < 0) { y->d[0] = - y->d[0]; n = -n; }
|
||||
return (n&3);
|
||||
@ -274,14 +274,14 @@ int mpranred(double x, mp_no *y, int p)
|
||||
/* Multi-Precision sin() function subroutine, for p=32. It is */
|
||||
/* based on the routines mpranred() and c32(). */
|
||||
/*******************************************************************/
|
||||
double mpsin1(double x)
|
||||
double __mpsin1(double x)
|
||||
{
|
||||
int p;
|
||||
int n;
|
||||
mp_no u,s,c;
|
||||
double y;
|
||||
p=32;
|
||||
n=mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */
|
||||
n=__mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */
|
||||
__c32(&u,&c,&s,p);
|
||||
switch (n) { /* in which quarter of unit circle y is*/
|
||||
case 0:
|
||||
@ -313,7 +313,7 @@ double mpsin1(double x)
|
||||
/* based on the routines mpranred() and c32(). */
|
||||
/*****************************************************************/
|
||||
|
||||
double mpcos1(double x)
|
||||
double __mpcos1(double x)
|
||||
{
|
||||
int p;
|
||||
int n;
|
||||
@ -321,7 +321,7 @@ double mpcos1(double x)
|
||||
double y;
|
||||
|
||||
p=32;
|
||||
n=mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */
|
||||
n=__mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */
|
||||
__c32(&u,&c,&s,p);
|
||||
switch (n) { /* in what quarter of unit circle y is*/
|
||||
|
||||
|
||||
@ -30,10 +30,10 @@
|
||||
/**************************************************************************/
|
||||
#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;
|
||||
@ -45,20 +45,20 @@ double slowexp(double x) {
|
||||
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);
|
||||
mul(&mpeps,&mpy,&mpcor,p);
|
||||
add(&mpy,&mpcor,&mpw,p);
|
||||
sub(&mpy,&mpcor,&mpz,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);
|
||||
__mp_dbl(&mpw, &w, p);
|
||||
__mp_dbl(&mpz, &z, 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;
|
||||
}
|
||||
|
||||
@ -34,40 +34,40 @@
|
||||
|
||||
#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);
|
||||
double __halfulp(double x,double y);
|
||||
|
||||
double slowpow(double x, double y, double z) {
|
||||
double __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}};
|
||||
int p;
|
||||
|
||||
res = halfulp(x,y); /* halfulp() returns -10 or x^y */
|
||||
res = __halfulp(x,y); /* halfulp() returns -10 or x^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 */
|
||||
mul(&mpy,&mpz,&mpw,p); /* y * z =w */
|
||||
mpexp(&mpw, &mpp, p); /* e^w =pp */
|
||||
add(&mpp,&eps,&mpr,p); /* pp+eps =r */
|
||||
__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 */
|
||||
__add(&mpp,&eps,&mpr,p); /* pp+eps =r */
|
||||
__mp_dbl(&mpr, &res, p);
|
||||
sub(&mpp,&eps,&mpr1,p); /* pp -eps =r1 */
|
||||
__sub(&mpp,&eps,&mpr1,p); /* pp -eps =r1 */
|
||||
__mp_dbl(&mpr1, &res1, p); /* converting into double precision */
|
||||
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 */
|
||||
mul(&mpy,&mpz,&mpw,p); /* y*z =w */
|
||||
mpexp(&mpw, &mpp, p); /* e^w=pp */
|
||||
__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 */
|
||||
__mp_dbl(&mpp, &res, p); /* converting into double precision */
|
||||
return res;
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user