diff options
author | Ulrich Drepper <drepper@redhat.com> | 2001-03-12 00:04:52 +0000 |
---|---|---|
committer | Ulrich Drepper <drepper@redhat.com> | 2001-03-12 00:04:52 +0000 |
commit | e4d8276142b9c07b23043ef44b0fe8fa7bcc3121 (patch) | |
tree | f153a80b6ce0fdd3261ff18a16fd80bd965231c3 /sysdeps/ieee754/dbl-64/dosincos.c | |
parent | d3c8723f6415af59a6ec14fcb918ad0e4d1fb588 (diff) | |
download | glibc-e4d8276142b9c07b23043ef44b0fe8fa7bcc3121.tar.gz |
Update.
2001-03-11 Ulrich Drepper <drepper@redhat.com>
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
Anna Reitman <reitman@il.ibm.com>.
* math/Makefile (dbl-only-routines): New variable.
(libm-routines): Add $(dbl-only-routines).
* sysdeps/ieee754/dbl-64/e_acos.c: Empty, definition is in e_asin.c.
* sysdeps/ieee754/dbl-64/e_asin.c: Replaced with accurate asin
implementation.
* sysdeps/ieee754/dbl-64/e_atan2.c: Replaced with accurate atan2
implementation.
* sysdeps/ieee754/dbl-64/e_exp.c: Replaced with accurate exp
implementation.
* sysdeps/ieee754/dbl-64/e_lgamma_r.c: Don't use __kernel_sin and
__kernel_cos.
* sysdeps/ieee754/dbl-64/e_log.c: Replaced with accurate log
implementation.
* sysdeps/ieee754/dbl-64/e_remainder.c: Replaced with accurate
remainder implementation.
* sysdeps/ieee754/dbl-64/e_pow.c: Replaced with accurate pow
implementation.
* sysdeps/ieee754/dbl-64/e_sqrt.c: Replaced with accurate sqrt
implementation.
* sysdeps/ieee754/dbl-64/k_cos.c: Empty, definition is in s_sin.c.
* sysdeps/ieee754/dbl-64/k_sin.c: Empty, definition is in s_sin.c.
* sysdeps/ieee754/dbl-64/s_atan.c: Replaced with accurate atan
implementation.
* sysdeps/ieee754/dbl-64/s_cos.c: Empty, definition is in s_sin.c.
* sysdeps/ieee754/dbl-64/s_sin.c: Replaced with accurate sin/cos
implementation.
* sysdeps/ieee754/dbl-64/s_sincos.c: Rewritten to not use __kernel_sin
and __kernel_cos.
* sysdeps/ieee754/dbl-64/s_tan.c: Replaced with accurate tan
implementation.
* sysdeps/ieee754/dbl-64/Dist: Add new non-code files.
* sysdeps/ieee754/dbl-64/MathLib.h: New file.
* sysdeps/ieee754/dbl-64/asincos.tbl: New file.
* sysdeps/ieee754/dbl-64/atnat.h: New file.
* sysdeps/ieee754/dbl-64/atnat2.h: New file.
* sysdeps/ieee754/dbl-64/branred.c: New file.
* sysdeps/ieee754/dbl-64/branred.h: New file.
* sysdeps/ieee754/dbl-64/dla.h: New file.
* sysdeps/ieee754/dbl-64/doasin.c: New file.
* sysdeps/ieee754/dbl-64/doasin.h: New file.
* sysdeps/ieee754/dbl-64/dosincos.c: New file.
* sysdeps/ieee754/dbl-64/dosincos.h: New file.
* sysdeps/ieee754/dbl-64/endian.h: New file.
* sysdeps/ieee754/dbl-64/halfulp.c: New file.
* sysdeps/ieee754/dbl-64/mpa.c: New file.
* sysdeps/ieee754/dbl-64/mpa.h: New file.
* sysdeps/ieee754/dbl-64/mpa2.h: New file.
* sysdeps/ieee754/dbl-64/mpatan.c: New file.
* sysdeps/ieee754/dbl-64/mpatan.h: New file.
* sysdeps/ieee754/dbl-64/mpatan2.c: New file.
* sysdeps/ieee754/dbl-64/mpexp.c: New file.
* sysdeps/ieee754/dbl-64/mpexp.h: New file.
* sysdeps/ieee754/dbl-64/mplog.c: New file.
* sysdeps/ieee754/dbl-64/mplog.h: New file.
* sysdeps/ieee754/dbl-64/mpsqrt.c: New file.
* sysdeps/ieee754/dbl-64/mpsqrt.h: New file.
* sysdeps/ieee754/dbl-64/mptan.c: New file.
* sysdeps/ieee754/dbl-64/mydefs.h: New file.
* sysdeps/ieee754/dbl-64/powtwo.tbl: New file.
* sysdeps/ieee754/dbl-64/root.tbl: New file.
* sysdeps/ieee754/dbl-64/sincos.tbl: New file.
* sysdeps/ieee754/dbl-64/sincos32.c: New file.
* sysdeps/ieee754/dbl-64/sincos32.h: New file.
* sysdeps/ieee754/dbl-64/slowexp.c: New file.
* sysdeps/ieee754/dbl-64/slowpow.c: New file.
* sysdeps/ieee754/dbl-64/uasncs.h: New file.
* sysdeps/ieee754/dbl-64/uatan.tbl: New file.
* sysdeps/ieee754/dbl-64/uexp.h: New file.
* sysdeps/ieee754/dbl-64/uexp.tbl: New file.
* sysdeps/ieee754/dbl-64/ulog.h: New file.
* sysdeps/ieee754/dbl-64/ulog.tbl: New file.
* sysdeps/ieee754/dbl-64/upow.h: New file.
* sysdeps/ieee754/dbl-64/upow.tbl: New file.
* sysdeps/ieee754/dbl-64/urem.h: New file.
* sysdeps/ieee754/dbl-64/uroot.h: New file.
* sysdeps/ieee754/dbl-64/usncs.h: New file.
* sysdeps/ieee754/dbl-64/utan.h: New file.
* sysdeps/ieee754/dbl-64/utan.tbl: New file.
* sysdeps/i386/fpu/branred.c: New file.
* sysdeps/i386/fpu/doasin.c: New file.
* sysdeps/i386/fpu/dosincos.c: New file.
* sysdeps/i386/fpu/halfulp.c: New file.
* sysdeps/i386/fpu/mpa.c: New file.
* sysdeps/i386/fpu/mpatan.c: New file.
* sysdeps/i386/fpu/mpatan2.c: New file.
* sysdeps/i386/fpu/mpexp.c: New file.
* sysdeps/i386/fpu/mplog.c: New file.
* sysdeps/i386/fpu/mpsqrt.c: New file.
* sysdeps/i386/fpu/mptan.c: New file.
* sysdeps/i386/fpu/sincos32.c: New file.
* sysdeps/i386/fpu/slowexp.c: New file.
* sysdeps/i386/fpu/slowpow.c: New file.
* sysdeps/ia64/fpu/branred.c: New file.
* sysdeps/ia64/fpu/doasin.c: New file.
* sysdeps/ia64/fpu/dosincos.c: New file.
* sysdeps/ia64/fpu/halfulp.c: New file.
* sysdeps/ia64/fpu/mpa.c: New file.
* sysdeps/ia64/fpu/mpatan.c: New file.
* sysdeps/ia64/fpu/mpatan2.c: New file.
* sysdeps/ia64/fpu/mpexp.c: New file.
* sysdeps/ia64/fpu/mplog.c: New file.
* sysdeps/ia64/fpu/mpsqrt.c: New file.
* sysdeps/ia64/fpu/mptan.c: New file.
* sysdeps/ia64/fpu/sincos32.c: New file.
* sysdeps/ia64/fpu/slowexp.c: New file.
* sysdeps/ia64/fpu/slowpow.c: New file.
* sysdeps/m68k/fpu/branred.c: New file.
* sysdeps/m68k/fpu/doasin.c: New file.
* sysdeps/m68k/fpu/dosincos.c: New file.
* sysdeps/m68k/fpu/halfulp.c: New file.
* sysdeps/m68k/fpu/mpa.c: New file.
* sysdeps/m68k/fpu/mpatan.c: New file.
* sysdeps/m68k/fpu/mpatan2.c: New file.
* sysdeps/m68k/fpu/mpexp.c: New file.
* sysdeps/m68k/fpu/mplog.c: New file.
* sysdeps/m68k/fpu/mpsqrt.c: New file.
* sysdeps/m68k/fpu/mptan.c: New file.
* sysdeps/m68k/fpu/sincos32.c: New file.
* sysdeps/m68k/fpu/slowexp.c: New file.
* sysdeps/m68k/fpu/slowpow.c: New file.
* iconvdata/gconv-modules: Add a number of alias, mostly for IBM
codepages.
Diffstat (limited to 'sysdeps/ieee754/dbl-64/dosincos.c')
-rw-r--r-- | sysdeps/ieee754/dbl-64/dosincos.c | 181 |
1 files changed, 181 insertions, 0 deletions
diff --git a/sysdeps/ieee754/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c new file mode 100644 index 0000000000..2a5f8fcc34 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/dosincos.c @@ -0,0 +1,181 @@ + +/* + * IBM Accurate Mathematical Library + * Copyright (c) International Business Machines Corp., 2001 + * + * 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 + * (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 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. + */ +/********************************************************************/ +/* */ +/* MODULE_NAME: dosincos.c */ +/* */ +/* */ +/* FUNCTIONS: dubsin */ +/* dubcos */ +/* docos */ +/* FILES NEEDED: endian.h mydefs.h dla.h dosincos.h */ +/* sincos.tbl */ +/* */ +/* Routines compute sin() and cos() as Double-Length numbers */ +/********************************************************************/ + + + +#include "endian.h" +#include "mydefs.h" +#include "sincos.tbl" +#include "dla.h" +#include "dosincos.h" +/***********************************************************************/ +/* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */ +/* as Double-Length number and store it at array v .It computes it by */ +/* arithmetic action on Double-Length numbers */ +/*(x+dx) between 0 and PI/4 */ +/***********************************************************************/ + +void dubsin(double x, double dx, double v[]) { + double xx,y,yy,z,zz,r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee, + sn,ssn,cs,ccs,ds,dss,dc,dcc; + mynumber u; + int4 k; + + u.x=x+big.x; + k = u.i[LOW_HALF]<<2; + x=x-(u.x-big.x); + d=x+dx; + dd=(x-d)+dx; + /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */ + MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc); + sn=sincos.x[k]; /* */ + ssn=sincos.x[k+1]; /* sin(Xi) and cos(Xi) */ + cs=sincos.x[k+2]; /* */ + ccs=sincos.x[k+3]; /* */ + MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* Taylor */ + ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); + MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* series */ + ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); + MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* for sin */ + MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); + ADD2(ds,dss,d,dd,ds,dss,r,s); /* ds=sin(t) */ + + MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor */ + ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); + MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* series */ + ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); + MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* for cos */ + ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); + MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* dc=cos(t) */ + + MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); + MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc); + SUB2(e,ee,dc,dcc,e,ee,r,s); + ADD2(e,ee,sn,ssn,e,ee,r,s); /* e+ee=sin(x+dx) */ + + v[0]=e; + v[1]=ee; +} +/**********************************************************************/ +/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */ +/* as Double-Length number and store it in array v .It computes it by */ +/* arithmetic action on Double-Length numbers */ +/*(x+dx) between 0 and PI/4 */ +/**********************************************************************/ + +void dubcos(double x, double dx, double v[]) { + double xx,y,yy,z,zz,r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee, + sn,ssn,cs,ccs,ds,dss,dc,dcc; + mynumber u; + int4 k; + u.x=x+big.x; + k = u.i[LOW_HALF]<<2; + x=x-(u.x-big.x); + d=x+dx; + dd=(x-d)+dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */ + MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc); + sn=sincos.x[k]; /* */ + ssn=sincos.x[k+1]; /* sin(Xi) and cos(Xi) */ + cs=sincos.x[k+2]; /* */ + ccs=sincos.x[k+3]; /* */ + MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); + ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); + MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); + ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); + MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); + MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); + ADD2(ds,dss,d,dd,ds,dss,r,s); + + MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); + ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); + MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); + ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); + MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); + ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); + MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); + + MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); + MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc); + + MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); + ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); + MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); + ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); + MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); + MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); + ADD2(ds,dss,d,dd,ds,dss,r,s); + MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); + ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); + MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); + ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); + MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); + ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); + MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); + MUL2(sn,ssn,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); + MUL2(dc,dcc,cs,ccs,dc,dcc,p,hx,tx,hy,ty,q,c,cc); + ADD2(e,ee,dc,dcc,e,ee,r,s); + SUB2(cs,ccs,e,ee,e,ee,r,s); + + v[0]=e; + v[1]=ee; +} +/**********************************************************************/ +/* 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[]) { + 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];} + 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];} + /* cos(x) = sin ( 90 - x ) */ + else {dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1]; + } + } + else { /* y>= 3/4 * PI */ + p=2.0*hp0.x-y; /* p = PI- y */ + yy=2.0*hp1.x-yy; + y=p+yy; + yy=(p-y)+yy; + dubcos(y,yy,w); + v[0]=-w[0]; + v[1]=-w[1]; + } +} |