diff options
| -rw-r--r-- | sysdeps/ieee754/dbl-64/doasin.c | 81 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/doasin.h | 63 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/dosincos.c | 217 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/dosincos.h | 80 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpa-arch.h | 47 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpa.c | 913 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpa.h | 123 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpatan.c | 116 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpatan.h | 145 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpatan2.c | 67 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpsqrt.c | 111 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpsqrt.h | 38 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mptan.c | 63 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/sincos32.c | 307 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/sincos32.h | 81 |
15 files changed, 0 insertions, 2452 deletions
diff --git a/sysdeps/ieee754/dbl-64/doasin.c b/sysdeps/ieee754/dbl-64/doasin.c deleted file mode 100644 index a65ef11477..0000000000 --- a/sysdeps/ieee754/dbl-64/doasin.c +++ /dev/null @@ -1,81 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ -/**********************************************************************/ -/* MODULE_NAME: doasin.c */ -/* */ -/* FUNCTION: doasin */ -/* */ -/* FILES NEEDED:endian.h mydefs.h dla.h doasin.h */ -/* mpa.c */ -/* */ -/* 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) */ -/**********************************************************************/ - -#include "endian.h" -#include "mydefs.h" -#include <dla.h> -#include <math_private.h> - -#ifndef SECTION -# define SECTION -#endif - -/********************************************************************/ -/* 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 -SECTION -__doasin(double x, double dx, double v[]) { - -#include "doasin.h" - - static const double - d5 = 0.22372159090911789889975459505194491E-01, - d6 = 0.17352764422456822913014975683014622E-01, - d7 = 0.13964843843786693521653681033981614E-01, - d8 = 0.11551791438485242609036067259086589E-01, - d9 = 0.97622386568166960207425666787248914E-02, - d10 = 0.83638737193775788576092749009744976E-02, - d11 = 0.79470250400727425881446981833568758E-02; - - double xx,p,pp,u,uu,r,s; - double tc,tcc; - - -/* Taylor series for arcsin for Double-Length numbers */ - xx = x*x+2.0*x*dx; - p = ((((((d11*xx+d10)*xx+d9)*xx+d8)*xx+d7)*xx+d6)*xx+d5)*xx; - pp = 0; - - MUL2(x,dx,x,dx,u,uu,tc,tcc); - ADD2(p,pp,c4.x,cc4.x,p,pp,r,s); - MUL2(p,pp,u,uu,p,pp,tc,tcc); - ADD2(p,pp,c3.x,cc3.x,p,pp,r,s); - MUL2(p,pp,u,uu,p,pp,tc,tcc); - ADD2(p,pp,c2.x,cc2.x,p,pp,r,s); - MUL2(p,pp,u,uu,p,pp,tc,tcc); - ADD2(p,pp,c1.x,cc1.x,p,pp,r,s); - MUL2(p,pp,u,uu,p,pp,tc,tcc); - MUL2(p,pp,x,dx,p,pp,tc,tcc); - ADD2(p,pp,x,dx,p,pp,r,s); - v[0]=p; - v[1]=pp; /* arcsin(x+dx)=v[0]+v[1] */ -} diff --git a/sysdeps/ieee754/dbl-64/doasin.h b/sysdeps/ieee754/dbl-64/doasin.h deleted file mode 100644 index bfabca70ae..0000000000 --- a/sysdeps/ieee754/dbl-64/doasin.h +++ /dev/null @@ -1,63 +0,0 @@ - -/* - * IBM Accurate Mathematical Library - * Written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ - -/************************************************************************/ -/* MODULE_NAME: doasin.h */ -/* */ -/* */ -/* common data and variables definition for BIG or LITTLE ENDIAN */ -/************************************************************************/ - - - -#ifndef DOASIN_H -#define DOASIN_H - -#ifdef BIG_ENDI - - static const mynumber -/**/ c1 = {{0x3FC55555, 0x55555555}}, /* 0.16666666666666666 */ -/**/ cc1 = {{0x3C655555, 0x55775389}}, /* 9.2518585419753846e-18 */ -/**/ c2 = {{0x3FB33333, 0x33333333}}, /* 0.074999999999999997 */ -/**/ cc2 = {{0x3C499993, 0x63F1A115}}, /* 2.7755472886508899e-18 */ -/**/ c3 = {{0x3FA6DB6D, 0xB6DB6DB7}}, /* 0.044642857142857144 */ -/**/ cc3 = {{0xBC320FC0, 0x3D5CF0C5}}, /* -9.7911734574147224e-19 */ -/**/ c4 = {{0x3F9F1C71, 0xC71C71C5}}, /* 0.030381944444444437 */ -/**/ cc4 = {{0xBC02B240, 0xFF23ED1E}}; /* -1.2669108566898312e-19 */ - -#else -#ifdef LITTLE_ENDI - - static const mynumber -/**/ c1 = {{0x55555555, 0x3FC55555}}, /* 0.16666666666666666 */ -/**/ cc1 = {{0x55775389, 0x3C655555}}, /* 9.2518585419753846e-18 */ -/**/ c2 = {{0x33333333, 0x3FB33333}}, /* 0.074999999999999997 */ -/**/ cc2 = {{0x63F1A115, 0x3C499993}}, /* 2.7755472886508899e-18 */ -/**/ c3 = {{0xB6DB6DB7, 0x3FA6DB6D}}, /* 0.044642857142857144 */ -/**/ cc3 = {{0x3D5CF0C5, 0xBC320FC0}}, /* -9.7911734574147224e-19 */ -/**/ c4 = {{0xC71C71C5, 0x3F9F1C71}}, /* 0.030381944444444437 */ -/**/ cc4 = {{0xFF23ED1E, 0xBC02B240}}; /* -1.2669108566898312e-19 */ - - -#endif -#endif - - -#endif diff --git a/sysdeps/ieee754/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c deleted file mode 100644 index 68e3a11401..0000000000 --- a/sysdeps/ieee754/dbl-64/dosincos.c +++ /dev/null @@ -1,217 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ -/********************************************************************/ -/* */ -/* 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 <dla.h> -#include "dosincos.h" -#include <math_private.h> - -#ifndef SECTION -# define SECTION -#endif - -extern const union -{ - int4 i[880]; - double x[440]; -} __sincostab attribute_hidden; - -/***********************************************************************/ -/* 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 -SECTION -__dubsin (double x, double dx, double v[]) -{ - double r, s, 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, c, cc); - sn = __sincostab.x[k]; /* */ - ssn = __sincostab.x[k + 1]; /* sin(Xi) and cos(Xi) */ - cs = __sincostab.x[k + 2]; /* */ - ccs = __sincostab.x[k + 3]; /* */ - /* Taylor series for sin ds=sin(t) */ - MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, c, cc); - ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s); - MUL2 (d2, dd2, ds, dss, ds, dss, c, cc); - ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s); - MUL2 (d2, dd2, ds, dss, ds, dss, c, cc); - MUL2 (d, dd, ds, dss, ds, dss, c, cc); - ADD2 (ds, dss, d, dd, ds, dss, r, s); - - /* Taylor series for cos dc=cos(t) */ - MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, c, cc); - ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s); - MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc); - ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s); - MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc); - ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s); - MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc); - - MUL2 (cs, ccs, ds, dss, e, ee, c, cc); - MUL2 (dc, dcc, sn, ssn, dc, dcc, 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 -SECTION -__dubcos (double x, double dx, double v[]) -{ - double r, s, 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, c, cc); - sn = __sincostab.x[k]; /* */ - ssn = __sincostab.x[k + 1]; /* sin(Xi) and cos(Xi) */ - cs = __sincostab.x[k + 2]; /* */ - ccs = __sincostab.x[k + 3]; /* */ - MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, c, cc); - ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s); - MUL2 (d2, dd2, ds, dss, ds, dss, c, cc); - ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s); - MUL2 (d2, dd2, ds, dss, ds, dss, c, cc); - MUL2 (d, dd, ds, dss, ds, dss, c, cc); - ADD2 (ds, dss, d, dd, ds, dss, r, s); - - MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, c, cc); - ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s); - MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc); - ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s); - MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc); - ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s); - MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc); - - MUL2 (cs, ccs, ds, dss, e, ee, c, cc); - MUL2 (dc, dcc, sn, ssn, dc, dcc, c, cc); - - MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, c, cc); - ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s); - MUL2 (d2, dd2, ds, dss, ds, dss, c, cc); - ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s); - MUL2 (d2, dd2, ds, dss, ds, dss, c, cc); - MUL2 (d, dd, ds, dss, ds, dss, c, cc); - ADD2 (ds, dss, d, dd, ds, dss, r, s); - MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, c, cc); - ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s); - MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc); - ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s); - MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc); - ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s); - MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc); - MUL2 (sn, ssn, ds, dss, e, ee, c, cc); - MUL2 (dc, dcc, cs, ccs, dc, dcc, 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 -SECTION -__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]; - } -} diff --git a/sysdeps/ieee754/dbl-64/dosincos.h b/sysdeps/ieee754/dbl-64/dosincos.h deleted file mode 100644 index 9f34339063..0000000000 --- a/sysdeps/ieee754/dbl-64/dosincos.h +++ /dev/null @@ -1,80 +0,0 @@ - -/* - * IBM Accurate Mathematical Library - * Written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ - -/************************************************************************/ -/* MODULE_NAME: dosincos.h */ -/* */ -/* */ -/* common data and variables definition for BIG or LITTLE ENDIAN */ -/************************************************************************/ - - - -#ifndef DOSINCOS_H -#define DOSINCOS_H - - -#ifdef BIG_ENDI -static const mynumber -/**/ s3 = {{0xBFC55555, 0x55555555}},/* -0.16666666666666666 */ -/**/ ss3 = {{0xBC6553AA, 0xE77EE482}},/* -9.2490366677784492e-18 */ -/**/ s5 = {{0x3F811111, 0x11110F15}},/* 0.008333333333332452 */ -/**/ ss5 = {{0xBC21AC06, 0xDA488820}},/* -4.7899996586987931e-19 */ -/**/ s7 = {{0xBF2A019F, 0x5816C78D}},/* -0.00019841261022928957 */ -/**/ ss7 = {{0x3BCDCEC9, 0x6A18BF2A}},/* 1.2624077757871259e-20 */ -/**/ c2 = {{0x3FE00000, 0x00000000}},/* 0.5 */ -/**/ cc2 = {{0xBA282FD8, 0x00000000}},/* -1.5264073330037701e-28 */ -/**/ c4 = {{0xBFA55555, 0x55555555}},/* -0.041666666666666664 */ -/**/ cc4 = {{0xBC4554BC, 0x2FFF257E}},/* -2.312711276085743e-18 */ -/**/ c6 = {{0x3F56C16C, 0x16C16A96}},/* 0.0013888888888888055 */ -/**/ cc6 = {{0xBBD2E846, 0xE6346F14}},/* -1.6015133010194884e-20 */ -/**/ c8 = {{0xBEFA019F, 0x821D5987}},/* -2.480157866754367e-05 */ -/**/ cc8 = {{0x3B7AB71E, 0x72FFE5CC}},/* 3.5357416224857556e-22 */ - -/**/ big = {{0x42c80000, 0x00000000}}, /* 52776558133248 */ - -/**/ hp0 = {{0x3FF921FB, 0x54442D18}}, /* PI / 2 */ -/**/ hp1 = {{0x3C91A626, 0x33145C07}}; /* 6.123233995736766e-17 */ -#else -#ifdef LITTLE_ENDI -static const mynumber -/**/ s3 = {{0x55555555, 0xBFC55555}},/* -0.16666666666666666 */ -/**/ ss3 = {{0xE77EE482, 0xBC6553AA}},/* -9.2490366677784492e-18 */ -/**/ s5 = {{0x11110F15, 0x3F811111}},/* 0.008333333333332452 */ -/**/ ss5 = {{0xDA488820, 0xBC21AC06}},/* -4.7899996586987931e-19 */ -/**/ s7 = {{0x5816C78D, 0xBF2A019F}},/* -0.00019841261022928957 */ -/**/ ss7 = {{0x6A18BF2A, 0x3BCDCEC9}},/* 1.2624077757871259e-20 */ -/**/ c2 = {{0x00000000, 0x3FE00000}},/* 0.5 */ -/**/ cc2 = {{0x00000000, 0xBA282FD8}},/* -1.5264073330037701e-28 */ -/**/ c4 = {{0x55555555, 0xBFA55555}},/* -0.041666666666666664 */ -/**/ cc4 = {{0x2FFF257E, 0xBC4554BC}},/* -2.312711276085743e-18 */ -/**/ c6 = {{0x16C16A96, 0x3F56C16C}},/* 0.0013888888888888055 */ -/**/ cc6 = {{0xE6346F14, 0xBBD2E846}},/* -1.6015133010194884e-20 */ -/**/ c8 = {{0x821D5987, 0xBEFA019F}},/* -2.480157866754367e-05 */ -/**/ cc8 = {{0x72FFE5CC, 0x3B7AB71E}},/* 3.5357416224857556e-22 */ - -/**/ big = {{0x00000000, 0x42c80000}}, /* 52776558133248 */ - -/**/ hp0 = {{0x54442D18, 0x3FF921FB}}, /* PI / 2 */ -/**/ hp1 = {{0x33145C07, 0x3C91A626}}; /* 6.123233995736766e-17 */ -#endif -#endif - -#endif diff --git a/sysdeps/ieee754/dbl-64/mpa-arch.h b/sysdeps/ieee754/dbl-64/mpa-arch.h deleted file mode 100644 index fbe296d8b5..0000000000 --- a/sysdeps/ieee754/dbl-64/mpa-arch.h +++ /dev/null @@ -1,47 +0,0 @@ -/* Overridable constants and operations. - Copyright (C) 2013-2021 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 - the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. */ - -#include <stdint.h> - -typedef long mantissa_t; -typedef int64_t mantissa_store_t; - -#define TWOPOW(i) (1L << i) - -#define RADIX_EXP 24 -#define RADIX TWOPOW (RADIX_EXP) /* 2^24 */ - -/* Divide D by RADIX and put the remainder in R. D must be a non-negative - integral value. */ -#define DIV_RADIX(d, r) \ - ({ \ - r = d & (RADIX - 1); \ - d >>= RADIX_EXP; \ - }) - -/* Put the integer component of a double X in R and retain the fraction in - X. This is used in extracting mantissa digits for MP_NO by using the - integer portion of the current value of the number as the current mantissa - digit and then scaling by RADIX to get the next mantissa digit in the same - manner. */ -#define INTEGER_OF(x, i) \ - ({ \ - i = (mantissa_t) x; \ - x -= i; \ - }) - -/* Align IN down to F. The code assumes that F is a power of two. */ -#define ALIGN_DOWN_TO(in, f) ((in) & - (f)) diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c deleted file mode 100644 index eb5d8e8e89..0000000000 --- a/sysdeps/ieee754/dbl-64/mpa.c +++ /dev/null @@ -1,913 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ -/************************************************************************/ -/* MODULE_NAME: mpa.c */ -/* */ -/* FUNCTIONS: */ -/* mcr */ -/* acr */ -/* cpy */ -/* norm */ -/* denorm */ -/* mp_dbl */ -/* dbl_mp */ -/* add_magnitudes */ -/* sub_magnitudes */ -/* add */ -/* sub */ -/* mul */ -/* inv */ -/* dvd */ -/* */ -/* Arithmetic functions for multiple precision numbers. */ -/* Relative errors are bounded */ -/************************************************************************/ - - -#include "endian.h" -#include "mpa.h" -#include <sys/param.h> -#include <alloca.h> - -#ifndef SECTION -# define SECTION -#endif - -#ifndef NO__CONST -const mp_no __mpone = { 1, { 1.0, 1.0 } }; -const mp_no __mptwo = { 1, { 1.0, 2.0 } }; -#endif - -#ifndef NO___ACR -/* Compare mantissa of two multiple precision numbers regardless of the sign - and exponent of the numbers. */ -static int -mcr (const mp_no *x, const mp_no *y, int p) -{ - long i; - long p2 = p; - for (i = 1; i <= p2; i++) - { - if (X[i] == Y[i]) - continue; - else if (X[i] > Y[i]) - return 1; - else - return -1; - } - return 0; -} - -/* Compare the absolute values of two multiple precision numbers. */ -int -__acr (const mp_no *x, const mp_no *y, int p) -{ - long i; - - if (X[0] == 0) - { - if (Y[0] == 0) - i = 0; - else - i = -1; - } - else if (Y[0] == 0) - i = 1; - else - { - if (EX > EY) - i = 1; - else if (EX < EY) - i = -1; - else - i = mcr (x, y, p); - } - - return i; -} -#endif - -#ifndef NO___CPY -/* Copy multiple precision number X into Y. They could be the same - number. */ -void -__cpy (const mp_no *x, mp_no *y, int p) -{ - long i; - - EY = EX; - for (i = 0; i <= p; i++) - Y[i] = X[i]; -} -#endif - -#ifndef NO___MP_DBL -/* Convert a multiple precision number *X into a double precision - number *Y, normalized case (|x| >= 2**(-1022))). X has precision - P, which is positive. */ -static void -norm (const mp_no *x, double *y, int p) -{ -# define R RADIXI - long i; - double c; - mantissa_t a, u, v, z[5]; - if (p < 5) - { - if (p == 1) - c = X[1]; - else if (p == 2) - c = X[1] + R * X[2]; - else if (p == 3) - c = X[1] + R * (X[2] + R * X[3]); - else /* p == 4. */ - c = (X[1] + R * X[2]) + R * R * (X[3] + R * X[4]); - } - else - { - for (a = 1, z[1] = X[1]; z[1] < TWO23; ) - { - a *= 2; - z[1] *= 2; - } - - for (i = 2; i < 5; i++) - { - mantissa_store_t d, r; - d = X[i] * (mantissa_store_t) a; - DIV_RADIX (d, r); - z[i] = r; - z[i - 1] += d; - } - - u = ALIGN_DOWN_TO (z[3], TWO19); - v = z[3] - u; - - if (v == TWO18) - { - if (z[4] == 0) - { - for (i = 5; i <= p; i++) - { - if (X[i] == 0) - continue; - else - { - z[3] += 1; - break; - } - } - } - else - z[3] += 1; - } - - c = (z[1] + R * (z[2] + R * z[3])) / a; - } - - c *= X[0]; - - for (i = 1; i < EX; i++) - c *= RADIX; - for (i = 1; i > EX; i--) - c *= RADIXI; - - *y = c; -# undef R -} - -/* Convert a multiple precision number *X into a double precision - number *Y, Denormal case (|x| < 2**(-1022))). */ -static void -denorm (const mp_no *x, double *y, int p) -{ - long i, k; - long p2 = p; - double c; - mantissa_t u, z[5]; - -# define R RADIXI - if (EX < -44 || (EX == -44 && X[1] < TWO5)) - { - *y = 0; - return; - } - - if (p2 == 1) - { - if (EX == -42) - { - z[1] = X[1] + TWO10; - z[2] = 0; - z[3] = 0; - k = 3; - } - else if (EX == -43) - { - z[1] = TWO10; - z[2] = X[1]; - z[3] = 0; - k = 2; - } - else - { - z[1] = TWO10; - z[2] = 0; - z[3] = X[1]; - k = 1; - } - } - else if (p2 == 2) - { - if (EX == -42) - { - z[1] = X[1] + TWO10; - z[2] = X[2]; - z[3] = 0; - k = 3; - } - else if (EX == -43) - { - z[1] = TWO10; - z[2] = X[1]; - z[3] = X[2]; - k = 2; - } - else - { - z[1] = TWO10; - z[2] = 0; - z[3] = X[1]; - k = 1; - } - } - else - { - if (EX == -42) - { - z[1] = X[1] + TWO10; - z[2] = X[2]; - k = 3; - } - else if (EX == -43) - { - z[1] = TWO10; - z[2] = X[1]; - k = 2; - } - else - { - z[1] = TWO10; - z[2] = 0; - k = 1; - } - z[3] = X[k]; - } - - u = ALIGN_DOWN_TO (z[3], TWO5); - - if (u == z[3]) - { - for (i = k + 1; i <= p2; i++) - { - if (X[i] == 0) - continue; - else - { - z[3] += 1; - break; - } - } - } - - c = X[0] * ((z[1] + R * (z[2] + R * z[3])) - TWO10); - - *y = c * TWOM1032; -# undef R -} - -/* Convert multiple precision number *X into double precision number *Y. The - result is correctly rounded to the nearest/even. */ -void -__mp_dbl (const mp_no *x, double *y, int p) -{ - if (X[0] == 0) - { - *y = 0; - return; - } - - if (__glibc_likely (EX > -42 || (EX == -42 && X[1] >= TWO10))) - norm (x, y, p); - else - denorm (x, y, p); -} -#endif - -/* Get the multiple precision equivalent of X into *Y. If the precision is too - small, the result is truncated. */ -void -SECTION -__dbl_mp (double x, mp_no *y, int p) -{ - long i, n; - long p2 = p; - - /* Sign. */ - if (x == 0) - { - Y[0] = 0; - return; - } - else if (x > 0) - Y[0] = 1; - else - { - Y[0] = -1; - x = -x; - } - - /* Exponent. */ - for (EY = 1; x >= RADIX; EY += 1) - x *= RADIXI; - for (; x < 1; EY -= 1) - x *= RADIX; - - /* Digits. */ - n = MIN (p2, 4); - for (i = 1; i <= n; i++) - { - INTEGER_OF (x, Y[i]); - x *= RADIX; - } - for (; i <= p2; i++) - Y[i] = 0; -} - -/* Add magnitudes of *X and *Y assuming that abs (*X) >= abs (*Y) > 0. The - sign of the sum *Z is not changed. X and Y may overlap but not X and Z or - Y and Z. No guard digit is used. The result equals the exact sum, - truncated. */ -static void -SECTION -add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p) -{ - long i, j, k; - long p2 = p; - mantissa_t zk; - - EZ = EX; - - i = p2; - j = p2 + EY - EX; - k = p2 + 1; - - if (__glibc_unlikely (j < 1)) - { - __cpy (x, z, p); - return; - } - - zk = 0; - - for (; j > 0; i--, j--) - { - zk += X[i] + Y[j]; - if (zk >= RADIX) - { - Z[k--] = zk - RADIX; - zk = 1; - } - else - { - Z[k--] = zk; - zk = 0; - } - } - - for (; i > 0; i--) - { - zk += X[i]; - if (zk >= RADIX) - { - Z[k--] = zk - RADIX; - zk = 1; - } - else - { - Z[k--] = zk; - zk = 0; - } - } - - if (zk == 0) - { - for (i = 1; i <= p2; i++) - Z[i] = Z[i + 1]; - } - else - { - Z[1] = zk; - EZ += 1; - } -} - -/* Subtract the magnitudes of *X and *Y assuming that abs (*x) > abs (*y) > 0. - The sign of the difference *Z is not changed. X and Y may overlap but not X - and Z or Y and Z. One guard digit is used. The error is less than one - ULP. */ -static void -SECTION -sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p) -{ - long i, j, k; - long p2 = p; - mantissa_t zk; - - EZ = EX; - i = p2; - j = p2 + EY - EX; - k = p2; - - /* Y is too small compared to X, copy X over to the result. */ - if (__glibc_unlikely (j < 1)) - { - __cpy (x, z, p); - return; - } - - /* The relevant least significant digit in Y is non-zero, so we factor it in - to enhance accuracy. */ - if (j < p2 && Y[j + 1] > 0) - { - Z[k + 1] = RADIX - Y[j + 1]; - zk = -1; - } - else - zk = Z[k + 1] = 0; - - /* Subtract and borrow. */ - for (; j > 0; i--, j--) - { - zk += (X[i] - Y[j]); - if (zk < 0) - { - Z[k--] = zk + RADIX; - zk = -1; - } - else - { - Z[k--] = zk; - zk = 0; - } - } - - /* We're done with digits from Y, so it's just digits in X. */ - for (; i > 0; i--) - { - zk += X[i]; - if (zk < 0) - { - Z[k--] = zk + RADIX; - zk = -1; - } - else - { - Z[k--] = zk; - zk = 0; - } - } - - /* Normalize. */ - for (i = 1; Z[i] == 0; i++) - ; - EZ = EZ - i + 1; - for (k = 1; i <= p2 + 1; ) - Z[k++] = Z[i++]; - for (; k <= p2; ) - Z[k++] = 0; -} - -/* Add *X and *Y and store the result in *Z. X and Y may overlap, but not X - and Z or Y and Z. One guard digit is used. The error is less than one - ULP. */ -void -SECTION -__add (const mp_no *x, const mp_no *y, mp_no *z, int p) -{ - int n; - - 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]; - } - 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]; - } - else if (n == -1) - { - sub_magnitudes (y, x, z, p); - Z[0] = Y[0]; - } - else - Z[0] = 0; - } -} - -/* Subtract *Y from *X and return the result in *Z. X and Y may overlap but - not X and Z or Y and Z. One guard digit is used. The error is less than - one ULP. */ -void -SECTION -__sub (const mp_no *x, const mp_no *y, mp_no *z, int p) -{ - int n; - - 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]; - } - 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]; - } - else if (n == -1) - { - sub_magnitudes (y, x, z, p); - Z[0] = -Y[0]; - } - else - Z[0] = 0; - } -} - -#ifndef NO__MUL -/* Multiply *X and *Y and store result in *Z. X and Y may overlap but not X - and Z or Y and Z. For P in [1, 2, 3], the exact result is truncated to P - digits. In case P > 3 the error is bounded by 1.001 ULP. */ -void -SECTION -__mul (const mp_no *x, const mp_no *y, mp_no *z, int p) -{ - long i, j, k, ip, ip2; - long p2 = p; - mantissa_store_t zk; - const mp_no *a; - mantissa_store_t *diag; - - /* Is z=0? */ - if (__glibc_unlikely (X[0] * Y[0] == 0)) - { - Z[0] = 0; - return; - } - - /* We need not iterate through all X's and Y's since it's pointless to - multiply zeroes. Here, both are zero... */ - for (ip2 = p2; ip2 > 0; ip2--) - if (X[ip2] != 0 || Y[ip2] != 0) - break; - - a = X[ip2] != 0 ? y : x; - - /* ... and here, at least one of them is still zero. */ - for (ip = ip2; ip > 0; ip--) - if (a->d[ip] != 0) - break; - - /* The product looks like this for p = 3 (as an example): - - - a1 a2 a3 - x b1 b2 b3 - ----------------------------- - a1*b3 a2*b3 a3*b3 - a1*b2 a2*b2 a3*b2 - a1*b1 a2*b1 a3*b1 - - So our K needs to ideally be P*2, but we're limiting ourselves to P + 3 - for P >= 3. We compute the above digits in two parts; the last P-1 - digits and then the first P digits. The last P-1 digits are a sum of - products of the input digits from P to P-k where K is 0 for the least - significant digit and increases as we go towards the left. The product - term is of the form X[k]*X[P-k] as can be seen in the above example. - - The first P digits are also a sum of products with the same product term, - except that the sum is from 1 to k. This is also evident from the above - example. - - Another thing that becomes evident is that only the most significant - ip+ip2 digits of the result are non-zero, where ip and ip2 are the - 'internal precision' of the input numbers, i.e. digits after ip and ip2 - are all 0. */ - - k = (__glibc_unlikely (p2 < 3)) ? p2 + p2 : p2 + 3; - - while (k > ip + ip2 + 1) - Z[k--] = 0; - - zk = 0; - - /* Precompute sums of diagonal elements so that we can directly use them - later. See the next comment to know we why need them. */ - diag = alloca (k * sizeof (mantissa_store_t)); - mantissa_store_t d = 0; - for (i = 1; i <= ip; i++) - { - d += X[i] * (mantissa_store_t) Y[i]; - diag[i] = d; - } - while (i < k) - diag[i++] = d; - - while (k > p2) - { - long lim = k / 2; - - if (k % 2 == 0) - /* We want to add this only once, but since we subtract it in the sum - of products above, we add twice. */ - zk += 2 * X[lim] * (mantissa_store_t) Y[lim]; - - for (i = k - p2, j = p2; i < j; i++, j--) - zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]); - - zk -= diag[k - 1]; - - DIV_RADIX (zk, Z[k]); - k--; - } - - /* The real deal. Mantissa digit Z[k] is the sum of all X[i] * Y[j] where i - goes from 1 -> k - 1 and j goes the same range in reverse. To reduce the - number of multiplications, we halve the range and if k is an even number, - add the diagonal element X[k/2]Y[k/2]. Through the half range, we compute - X[i] * Y[j] as (X[i] + X[j]) * (Y[i] + Y[j]) - X[i] * Y[i] - X[j] * Y[j]. - - This reduction tells us that we're summing two things, the first term - through the half range and the negative of the sum of the product of all - terms of X and Y in the full range. i.e. - - SUM(X[i] * Y[i]) for k terms. This is precalculated above for each k in - a single loop so that it completes in O(n) time and can hence be directly - used in the loop below. */ - while (k > 1) - { - long lim = k / 2; - - if (k % 2 == 0) - /* We want to add this only once, but since we subtract it in the sum - of products above, we add twice. */ - zk += 2 * X[lim] * (mantissa_store_t) Y[lim]; - - for (i = 1, j = k - 1; i < j; i++, j--) - zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]); - - zk -= diag[k - 1]; - - DIV_RADIX (zk, Z[k]); - k--; - } - Z[k] = zk; - - /* Get the exponent sum into an intermediate variable. This is a subtle - optimization, where given enough registers, all operations on the exponent - happen in registers and the result is written out only once into EZ. */ - int e = EX + EY; - - /* Is there a carry beyond the most significant digit? */ - if (__glibc_unlikely (Z[1] == 0)) - { - for (i = 1; i <= p2; i++) - Z[i] = Z[i + 1]; - e--; - } - - EZ = e; - Z[0] = X[0] * Y[0]; -} -#endif - -#ifndef NO__SQR -/* Square *X and store result in *Y. X and Y may not overlap. For P in - [1, 2, 3], the exact result is truncated to P digits. In case P > 3 the - error is bounded by 1.001 ULP. This is a faster special case of - multiplication. */ -void -SECTION -__sqr (const mp_no *x, mp_no *y, int p) -{ - long i, j, k, ip; - mantissa_store_t yk; - - /* Is z=0? */ - if (__glibc_unlikely (X[0] == 0)) - { - Y[0] = 0; - return; - } - - /* We need not iterate through all X's since it's pointless to - multiply zeroes. */ - for (ip = p; ip > 0; ip--) - if (X[ip] != 0) - break; - - k = (__glibc_unlikely (p < 3)) ? p + p : p + 3; - - while (k > 2 * ip + 1) - Y[k--] = 0; - - yk = 0; - - while (k > p) - { - mantissa_store_t yk2 = 0; - long lim = k / 2; - - if (k % 2 == 0) - yk += X[lim] * (mantissa_store_t) X[lim]; - - /* In __mul, this loop (and the one within the next while loop) run - between a range to calculate the mantissa as follows: - - Z[k] = X[k] * Y[n] + X[k+1] * Y[n-1] ... + X[n-1] * Y[k+1] - + X[n] * Y[k] - - For X == Y, we can get away with summing halfway and doubling the - result. For cases where the range size is even, the mid-point needs - to be added separately (above). */ - for (i = k - p, j = p; i < j; i++, j--) - yk2 += X[i] * (mantissa_store_t) X[j]; - - yk += 2 * yk2; - - DIV_RADIX (yk, Y[k]); - k--; - } - - while (k > 1) - { - mantissa_store_t yk2 = 0; - long lim = k / 2; - - if (k % 2 == 0) - yk += X[lim] * (mantissa_store_t) X[lim]; - - /* Likewise for this loop. */ - for (i = 1, j = k - 1; i < j; i++, j--) - yk2 += X[i] * (mantissa_store_t) X[j]; - - yk += 2 * yk2; - - DIV_RADIX (yk, Y[k]); - k--; - } - Y[k] = yk; - - /* Squares are always positive. */ - Y[0] = 1; - - /* Get the exponent sum into an intermediate variable. This is a subtle - optimization, where given enough registers, all operations on the exponent - happen in registers and the result is written out only once into EZ. */ - int e = EX * 2; - - /* Is there a carry beyond the most significant digit? */ - if (__glibc_unlikely (Y[1] == 0)) - { - for (i = 1; i <= p; i++) - Y[i] = Y[i + 1]; - e--; - } - - EY = e; -} -#endif - -/* Invert *X and store in *Y. Relative error bound: - - For P = 2: 1.001 * R ^ (1 - P) - - For P = 3: 1.063 * R ^ (1 - P) - - For P > 3: 2.001 * R ^ (1 - P) - - *X = 0 is not permissible. */ -static void -SECTION -__inv (const mp_no *x, mp_no *y, int p) -{ - long i; - double t; - mp_no z, w; - static const int np1[] = - { 0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 - }; - - __cpy (x, &z, p); - z.e = 0; - __mp_dbl (&z, &t, p); - t = 1 / t; - - /* t == 0 will never happen at this point, since 1/t can only be 0 if t is - infinity, but before the division t == mantissa of x (exponent is 0). We - can instruct the compiler to ignore this case. */ - if (t == 0) - __builtin_unreachable (); - - __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); - } -} - -/* Divide *X by *Y and store result in *Z. X and Y may overlap but not X and Z - or Y and Z. Relative error bound: - - For P = 2: 2.001 * R ^ (1 - P) - - For P = 3: 2.063 * R ^ (1 - P) - - For P > 3: 3.001 * R ^ (1 - P) - - *X = 0 is not permissible. */ -void -SECTION -__dvd (const mp_no *x, const mp_no *y, mp_no *z, int p) -{ - mp_no w; - - if (X[0] == 0) - Z[0] = 0; - else - { - __inv (y, &w, p); - __mul (x, &w, z, p); - } -} diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h deleted file mode 100644 index c28630e148..0000000000 --- a/sysdeps/ieee754/dbl-64/mpa.h +++ /dev/null @@ -1,123 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * Written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ - -/************************************************************************/ -/* MODULE_NAME: mpa.h */ -/* */ -/* FUNCTIONS: */ -/* mcr */ -/* acr */ -/* cpy */ -/* mp_dbl */ -/* dbl_mp */ -/* add */ -/* sub */ -/* mul */ -/* dvd */ -/* */ -/* Arithmetic functions for multiple precision numbers. */ -/* Common types and definition */ -/************************************************************************/ - -#include <mpa-arch.h> - -/* The mp_no structure holds the details of a multi-precision floating point - number. - - - The radix of the number (R) is 2 ^ 24. - - - E: The exponent of the number. - - - D[0]: The sign (-1, 1) or 0 if the value is 0. In the latter case, the - values of the remaining members of the structure are ignored. - - - D[1] - D[p]: The mantissa of the number where: - - 0 <= D[i] < R and - P is the precision of the number and 1 <= p <= 32 - - D[p+1] ... D[39] have no significance. - - - The value of the number is: - - D[1] * R ^ (E - 1) + D[2] * R ^ (E - 2) ... D[p] * R ^ (E - p) - - */ -typedef struct -{ - int e; - mantissa_t d[40]; -} mp_no; - -typedef union -{ - int i[2]; - double d; -} number; - -extern const mp_no __mpone; -extern const mp_no __mptwo; - -#define X x->d -#define Y y->d -#define Z z->d -#define EX x->e -#define EY y->e -#define EZ z->e - -#ifndef RADIXI -# define RADIXI 0x1.0p-24 /* 2^-24 */ -#endif - -#ifndef TWO52 -# define TWO52 0x1.0p52 /* 2^52 */ -#endif - -#define TWO5 TWOPOW (5) /* 2^5 */ -#define TWO8 TWOPOW (8) /* 2^52 */ -#define TWO10 TWOPOW (10) /* 2^10 */ -#define TWO18 TWOPOW (18) /* 2^18 */ -#define TWO19 TWOPOW (19) /* 2^19 */ -#define TWO23 TWOPOW (23) /* 2^23 */ - -#define HALFRAD TWO23 - -#define TWO57 0x1.0p57 /* 2^57 */ -#define TWO71 0x1.0p71 /* 2^71 */ -#define TWOM1032 0x1.0p-1032 /* 2^-1032 */ -#define TWOM1022 0x1.0p-1022 /* 2^-1022 */ - -#define HALF 0x1.0p-1 /* 1/2 */ -#define MHALF -0x1.0p-1 /* -1/2 */ - -int __acr (const mp_no *, const mp_no *, int); -void __cpy (const mp_no *, 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 __sqr (const mp_no *, mp_no *, int); -void __dvd (const mp_no *, const mp_no *, mp_no *, int); - -extern void __mpatan (mp_no *, mp_no *, int); -extern void __mpatan2 (mp_no *, mp_no *, mp_no *, int); -extern void __mpsqrt (mp_no *, mp_no *, int); -extern void __c32 (mp_no *, mp_no *, mp_no *, int); -extern int __mpranred (double, mp_no *, int); diff --git a/sysdeps/ieee754/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c deleted file mode 100644 index 77f84ddbf2..0000000000 --- a/sysdeps/ieee754/dbl-64/mpatan.c +++ /dev/null @@ -1,116 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ -/******************************************************************/ -/* */ -/* MODULE_NAME:mpatan.c */ -/* */ -/* FUNCTIONS:mpatan */ -/* */ -/* FILES NEEDED: mpa.h endian.h mpatan.h */ -/* mpa.c */ -/* */ -/* Multi-Precision Atan function subroutine, for precision p >= 4.*/ -/* The relative error of the result is bounded by 34.32*r**(1-p), */ -/* where r=2**24. */ -/******************************************************************/ - -#include "endian.h" -#include "mpa.h" -#include <math.h> - -#ifndef SECTION -# define SECTION -#endif - -#include "mpatan.h" - -void -SECTION -__mpatan (mp_no *x, mp_no *y, int p) -{ - int i, m, n; - double dx; - mp_no 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 - } - }; - - mp_no mps, mpsm, mpt, mpt1, mpt2, mpt3; - - /* Choose m and initiate mptwoim1. */ - if (EX > 0) - m = 7; - else if (EX < 0) - m = 0; - else - { - __mp_dbl (x, &dx, p); - dx = fabs (dx); - for (m = 6; m > 0; m--) - { - if (dx > __atan_xm[m].d) - break; - } - } - mptwoim1.e = 1; - mptwoim1.d[0] = 1; - - /* Reduce x m times. */ - __sqr (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); - } - __mpsqrt (&mpsm, &mps, p); - mps.d[0] = X[0]; - } - - /* Evaluate a truncated power series for Atan(s). */ - 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] -= 2; - __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); - - /* Compute Atan(x). */ - mptwoim1.d[1] = 1 << m; - __mul (&mptwoim1, &mpt, y, p); -} diff --git a/sysdeps/ieee754/dbl-64/mpatan.h b/sysdeps/ieee754/dbl-64/mpatan.h deleted file mode 100644 index 5f866a77f7..0000000000 --- a/sysdeps/ieee754/dbl-64/mpatan.h +++ /dev/null @@ -1,145 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * Written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ - -/******************************************************************/ -/* */ -/* MODULE_NAME:mpatan.h */ -/* */ -/* common data and variables prototype and definition */ -/******************************************************************/ - -#ifndef MPATAN_H -#define MPATAN_H - -extern const number __atan_xm[8] attribute_hidden; -extern const number __atan_twonm1[33] attribute_hidden; -extern const number __atan_twom[8] attribute_hidden; -extern const int __atan_np[33] attribute_hidden; - - -#ifndef AVOID_MPATAN_H -#ifdef BIG_ENDI - const number - __atan_xm[8] = { /* x[m] */ -/**/ {{0x00000000, 0x00000000} }, /* 0.0 */ -/**/ {{0x3f8930be, 0x00000000} }, /* 0.0123 */ -/**/ {{0x3f991687, 0x00000000} }, /* 0.0245 */ -/**/ {{0x3fa923a2, 0x00000000} }, /* 0.0491 */ -/**/ {{0x3fb930be, 0x00000000} }, /* 0.0984 */ -/**/ {{0x3fc95810, 0x00000000} }, /* 0.198 */ -/**/ {{0x3fda7ef9, 0x00000000} }, /* 0.414 */ -/**/ {{0x3ff00000, 0x00000000} }, /* 1.0 */ - }; - const number - __atan_twonm1[33] = { /* 2n-1 */ -/**/ {{0x00000000, 0x00000000} }, /* 0 */ -/**/ {{0x00000000, 0x00000000} }, /* 0 */ -/**/ {{0x00000000, 0x00000000} }, /* 0 */ -/**/ {{0x00000000, 0x00000000} }, /* 0 */ -/**/ {{0x40260000, 0x00000000} }, /* 11 */ -/**/ {{0x402e0000, 0x00000000} }, /* 15 */ -/**/ {{0x40330000, 0x00000000} }, /* 19 */ -/**/ {{0x40350000, 0x00000000} }, /* 21 */ -/**/ {{0x40390000, 0x00000000} }, /* 25 */ -/**/ {{0x403d0000, 0x00000000} }, /* 29 */ -/**/ {{0x40408000, 0x00000000} }, /* 33 */ -/**/ {{0x40428000, 0x00000000} }, /* 37 */ -/**/ {{0x40448000, 0x00000000} }, /* 41 */ -/**/ {{0x40468000, 0x00000000} }, /* 45 */ -/**/ {{0x40488000, 0x00000000} }, /* 49 */ -/**/ {{0x404a8000, 0x00000000} }, /* 53 */ -/**/ {{0x404b8000, 0x00000000} }, /* 55 */ -/**/ {{0x404d8000, 0x00000000} }, /* 59 */ -/**/ {{0x404f8000, 0x00000000} }, /* 63 */ -/**/ {{0x4050c000, 0x00000000} }, /* 67 */ -/**/ {{0x4051c000, 0x00000000} }, /* 71 */ -/**/ {{0x4052c000, 0x00000000} }, /* 75 */ -/**/ {{0x4053c000, 0x00000000} }, /* 79 */ -/**/ {{0x4054c000, 0x00000000} }, /* 83 */ -/**/ {{0x40554000, 0x00000000} }, /* 85 */ -/**/ {{0x40564000, 0x00000000} }, /* 89 */ -/**/ {{0x40574000, 0x00000000} }, /* 93 */ -/**/ {{0x40584000, 0x00000000} }, /* 97 */ -/**/ {{0x40594000, 0x00000000} }, /* 101 */ -/**/ {{0x405a4000, 0x00000000} }, /* 105 */ -/**/ {{0x405b4000, 0x00000000} }, /* 109 */ -/**/ {{0x405c4000, 0x00000000} }, /* 113 */ -/**/ {{0x405d4000, 0x00000000} }, /* 117 */ - }; - -#else -#ifdef LITTLE_ENDI - - const number - __atan_xm[8] = { /* x[m] */ -/**/ {{0x00000000, 0x00000000} }, /* 0.0 */ -/**/ {{0x00000000, 0x3f8930be} }, /* 0.0123 */ -/**/ {{0x00000000, 0x3f991687} }, /* 0.0245 */ -/**/ {{0x00000000, 0x3fa923a2} }, /* 0.0491 */ -/**/ {{0x00000000, 0x3fb930be} }, /* 0.0984 */ -/**/ {{0x00000000, 0x3fc95810} }, /* 0.198 */ -/**/ {{0x00000000, 0x3fda7ef9} }, /* 0.414 */ -/**/ {{0x00000000, 0x3ff00000} }, /* 1.0 */ - }; - const number -__atan_twonm1[33] = { /* 2n-1 */ -/**/ {{0x00000000, 0x00000000} }, /* 0 */ -/**/ {{0x00000000, 0x00000000} }, /* 0 */ -/**/ {{0x00000000, 0x00000000} }, /* 0 */ -/**/ {{0x00000000, 0x00000000} }, /* 0 */ -/**/ {{0x00000000, 0x40260000} }, /* 11 */ -/**/ {{0x00000000, 0x402e0000} }, /* 15 */ -/**/ {{0x00000000, 0x40330000} }, /* 19 */ -/**/ {{0x00000000, 0x40350000} }, /* 21 */ -/**/ {{0x00000000, 0x40390000} }, /* 25 */ -/**/ {{0x00000000, 0x403d0000} }, /* 29 */ -/**/ {{0x00000000, 0x40408000} }, /* 33 */ -/**/ {{0x00000000, 0x40428000} }, /* 37 */ -/**/ {{0x00000000, 0x40448000} }, /* 41 */ -/**/ {{0x00000000, 0x40468000} }, /* 45 */ -/**/ {{0x00000000, 0x40488000} }, /* 49 */ -/**/ {{0x00000000, 0x404a8000} }, /* 53 */ -/**/ {{0x00000000, 0x404b8000} }, /* 55 */ -/**/ {{0x00000000, 0x404d8000} }, /* 59 */ -/**/ {{0x00000000, 0x404f8000} }, /* 63 */ -/**/ {{0x00000000, 0x4050c000} }, /* 67 */ -/**/ {{0x00000000, 0x4051c000} }, /* 71 */ -/**/ {{0x00000000, 0x4052c000} }, /* 75 */ -/**/ {{0x00000000, 0x4053c000} }, /* 79 */ -/**/ {{0x00000000, 0x4054c000} }, /* 83 */ -/**/ {{0x00000000, 0x40554000} }, /* 85 */ -/**/ {{0x00000000, 0x40564000} }, /* 89 */ -/**/ {{0x00000000, 0x40574000} }, /* 93 */ -/**/ {{0x00000000, 0x40584000} }, /* 97 */ -/**/ {{0x00000000, 0x40594000} }, /* 101 */ -/**/ {{0x00000000, 0x405a4000} }, /* 105 */ -/**/ {{0x00000000, 0x405b4000} }, /* 109 */ -/**/ {{0x00000000, 0x405c4000} }, /* 113 */ -/**/ {{0x00000000, 0x405d4000} }, /* 117 */ - }; - -#endif -#endif - - const int - __atan_np[33] = { 0, 0, 0, 0, 6, 8,10,11,13,15,17,19,21,23,25,27,28, - 30,32,34,36,38,40,42,43,45,47,49,51,53,55,57,59}; - -#endif -#endif diff --git a/sysdeps/ieee754/dbl-64/mpatan2.c b/sysdeps/ieee754/dbl-64/mpatan2.c deleted file mode 100644 index 68348145af..0000000000 --- a/sysdeps/ieee754/dbl-64/mpatan2.c +++ /dev/null @@ -1,67 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ -/******************************************************************/ -/* MODULE_NAME: mpatan2.c */ -/* */ -/* FUNCTIONS:mpatan2 */ -/* */ -/* FILES NEEDED: mpa.h */ -/* mpa.c mpatan.c mpsqrt.c */ -/* */ -/* Multi-Precision Atan2(y,x) function subroutine, */ -/* for precision p >= 4. */ -/* y=0 is not permitted if x<=0. No error messages are given. */ -/* The relative error of the result is bounded by 44.84*r**(1-p) */ -/* if x <= 0, y != 0 and by 37.33*r**(1-p) if x>0. here r=2**24. */ -/* */ -/******************************************************************/ - -#include "mpa.h" - -#ifndef SECTION -# define SECTION -#endif - -/* Multi-Precision Atan2 (y, x) function subroutine, for p >= 4. - y = 0 is not permitted if x <= 0. No error messages are given. */ -void -SECTION -__mpatan2 (mp_no *y, mp_no *x, mp_no *z, int p) -{ - mp_no mpt1, mpt2, mpt3; - - if (X[0] <= 0) - { - __dvd (x, y, &mpt1, p); - __mul (&mpt1, &mpt1, &mpt2, p); - 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); - } - else - { - __dvd (y, x, &mpt1, p); - __mpatan (&mpt1, z, p); - } -} diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c deleted file mode 100644 index 2c32490c1f..0000000000 --- a/sysdeps/ieee754/dbl-64/mpsqrt.c +++ /dev/null @@ -1,111 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ -/****************************************************************************/ -/* MODULE_NAME:mpsqrt.c */ -/* */ -/* FUNCTION:mpsqrt */ -/* fastiroot */ -/* */ -/* FILES NEEDED:endian.h mpa.h mpsqrt.h */ -/* mpa.c */ -/* Multi-Precision square root function subroutine for precision p >= 4. */ -/* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */ -/* */ -/****************************************************************************/ -#include "endian.h" -#include "mpa.h" - -#ifndef SECTION -# define SECTION -#endif - -#include "mpsqrt.h" - -/****************************************************************************/ -/* Multi-Precision square root function subroutine for precision p >= 4. */ -/* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */ -/* Routine receives two pointers to Multi Precision numbers: */ -/* x (left argument) and y (next argument). Routine also receives precision */ -/* p as integer. Routine computes sqrt(*x) and stores result in *y */ -/****************************************************************************/ - -static double fastiroot (double); - -void -SECTION -__mpsqrt (mp_no *x, mp_no *y, int p) -{ - int i, m, ey; - double dx, dy; - static const mp_no mphalf = {0, {1.0, HALFRAD}}; - static const mp_no mp3halfs = {1, {1.0, 1.0, HALFRAD}}; - mp_no mpxn, mpz, mpu, mpt1, mpt2; - - 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 = __mpsqrt_mp[p]; - for (i = 0; i < m; i++) - { - __sqr (&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; -} - -/***********************************************************/ -/* Compute a double precision approximation for 1/sqrt(x) */ -/* with the relative error bounded by 2**-51. */ -/***********************************************************/ -static double -SECTION -fastiroot (double x) -{ - union - { - int i[2]; - double d; - } p, q; - double y, z, t; - int n; - static const double c0 = 0.99674, c1 = -0.53380; - static const double c2 = 0.45472, c3 = -0.21553; - - p.d = x; - p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF) | 0x3FE00000; - q.d = x; - y = p.d; - z = y - 1.0; - n = (q.i[HIGH_HALF] - p.i[HIGH_HALF]) >> 1; - z = ((c3 * z + c2) * z + c1) * z + c0; /* 2**-7 */ - z = z * (1.5 - 0.5 * y * z * z); /* 2**-14 */ - p.d = z * (1.5 - 0.5 * y * z * z); /* 2**-28 */ - p.i[HIGH_HALF] -= n; - t = x * p.d; - return p.d * (1.5 - 0.5 * p.d * t); -} diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.h b/sysdeps/ieee754/dbl-64/mpsqrt.h deleted file mode 100644 index d66fc99395..0000000000 --- a/sysdeps/ieee754/dbl-64/mpsqrt.h +++ /dev/null @@ -1,38 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * Written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ - -/******************************************************************/ -/* */ -/* MODULE_NAME:mpatan.h */ -/* */ -/* common data and variables prototype and definition */ -/******************************************************************/ - -#ifndef MPSQRT_H -#define MPSQRT_H - -extern const int __mpsqrt_mp[33] attribute_hidden; - - -#ifndef AVOID_MPSQRT_H - const int __mpsqrt_mp[33] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4, - 4,4,4,4,4,4,4,4,4}; -#endif - -#endif diff --git a/sysdeps/ieee754/dbl-64/mptan.c b/sysdeps/ieee754/dbl-64/mptan.c deleted file mode 100644 index 6cc3ee58a0..0000000000 --- a/sysdeps/ieee754/dbl-64/mptan.c +++ /dev/null @@ -1,63 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ -/**********************************************************************/ -/* MODULE_NAME:mptan.c */ -/* */ -/* FUNCTION: mptan */ -/* */ -/* FILES NEEDED: endian.h mpa.h */ -/* mpa.c sincos32.c branred.c */ -/* */ -/* Multi-Precision tan() function subroutine, for p=32. It is based */ -/* on the routines mpranred() and c32(). mpranred() performs range */ -/* reduction of a double number x into a multiple precision number */ -/* y, such that y=x-n*pi/2, abs(y)<pi/4, n=0,+-1,+-2,.... c32() */ -/* computes both sin(y), cos(y). tan(x) is either sin(y)/cos(y) */ -/* or -cos(y)/sin(y). The precision of the result is of about 559 */ -/* significant bits. */ -/* */ -/**********************************************************************/ -#include "endian.h" -#include "mpa.h" - -#ifndef SECTION -# define SECTION -#endif - -void -SECTION -__mptan (double x, mp_no *mpy, int p) -{ - int n; - mp_no mpw, mpc, mps; - - /* Negative or positive result. */ - n = __mpranred (x, &mpw, p) & 0x00000001; - /* Computing sin(x) and cos(x). */ - __c32 (&mpw, &mpc, &mps, p); - /* Second or fourth quarter of unit circle. */ - if (n) - { - __dvd (&mpc, &mps, mpy, p); - 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 deleted file mode 100644 index 44a313ad76..0000000000 --- a/sysdeps/ieee754/dbl-64/sincos32.c +++ /dev/null @@ -1,307 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ -/****************************************************************/ -/* MODULE_NAME: sincos32.c */ -/* */ -/* FUNCTIONS: ss32 */ -/* cc32 */ -/* c32 */ -/* sin32 */ -/* cos32 */ -/* mpsin */ -/* mpcos */ -/* mpranred */ -/* mpsin1 */ -/* mpcos1 */ -/* */ -/* FILES NEEDED: endian.h mpa.h sincos32.h */ -/* mpa.c */ -/* */ -/* Multi Precision sin() and cos() function with p=32 for sin()*/ -/* cos() arcsin() and arccos() routines */ -/* In addition mpranred() routine performs range reduction of */ -/* a double number x into multi precision number y, */ -/* such that y=x-n*pi/2, abs(y)<pi/4, n=0,+-1,+-2,.... */ -/****************************************************************/ -#include "endian.h" -#include "mpa.h" -#include "sincos32.h" -#include <math.h> -#include <math_private.h> -#include <stap-probe.h> - -#ifndef SECTION -# define SECTION -#endif - -/* Compute Multi-Precision sin() function for given p. Receive Multi Precision - number x and result stored at y. */ -static void -SECTION -ss32 (mp_no *x, mp_no *y, int p) -{ - int i; - double a; - mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}}; - for (i = 1; i <= p; i++) - mpk.d[i] = 0; - - __sqr (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 (x, &sum, y, p); -} - -/* Compute Multi-Precision cos() function for given p. Receive Multi Precision - number x and result stored at y. */ -static void -SECTION -cc32 (mp_no *x, mp_no *y, int p) -{ - int i; - double a; - mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}}; - for (i = 1; i <= p; i++) - mpk.d[i] = 0; - - __sqr (x, &x2, p); - mpk.d[1] = 27.0; - __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 (&x2, &sum, y, p); -} - -/* Compute both sin(x), cos(x) as Multi precision numbers. */ -void -SECTION -__c32 (mp_no *x, mp_no *y, mp_no *z, int p) -{ - mp_no u, t, t1, t2, c, s; - int i; - __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 (&__mptwo, &c, &t1, p); - __mul (&t1, &c, &t2, p); - __add (&t2, &t2, &c, p); - } - __sub (&__mpone, &c, y, p); - __cpy (&s, z, p); -} - -/* Compute sin() of double-length number (X + DX) as Multi Precision number and - return result as double. If REDUCE_RANGE is true, X is assumed to be the - original input and DX is ignored. */ -double -SECTION -__mpsin (double x, double dx, bool reduce_range) -{ - double y; - mp_no a, b, c, s; - int n; - int p = 32; - - if (reduce_range) - { - n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */ - __c32 (&a, &c, &s, p); - } - else - { - n = -1; - __dbl_mp (x, &b, p); - __dbl_mp (dx, &c, p); - __add (&b, &c, &a, p); - if (x > 0.8) - { - __sub (&hp, &a, &b, p); - __c32 (&b, &s, &c, p); - } - else - __c32 (&a, &c, &s, p); /* b = sin(x+dx) */ - } - - /* Convert result based on which quarter of unit circle y is in. */ - switch (n) - { - case 1: - __mp_dbl (&c, &y, p); - break; - - case 3: - __mp_dbl (&c, &y, p); - y = -y; - break; - - case 2: - __mp_dbl (&s, &y, p); - y = -y; - break; - - /* Quadrant not set, so the result must be sin (X + DX), which is also in - S. */ - case 0: - default: - __mp_dbl (&s, &y, p); - } - LIBC_PROBE (slowsin, 3, &x, &dx, &y); - return y; -} - -/* Compute cos() of double-length number (X + DX) as Multi Precision number and - return result as double. If REDUCE_RANGE is true, X is assumed to be the - original input and DX is ignored. */ -double -SECTION -__mpcos (double x, double dx, bool reduce_range) -{ - double y; - mp_no a, b, c, s; - int n; - int p = 32; - - if (reduce_range) - { - n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */ - __c32 (&a, &c, &s, p); - } - else - { - n = -1; - __dbl_mp (x, &b, p); - __dbl_mp (dx, &c, p); - __add (&b, &c, &a, p); - if (x > 0.8) - { - __sub (&hp, &a, &b, p); - __c32 (&b, &s, &c, p); - } - else - __c32 (&a, &c, &s, p); /* a = cos(x+dx) */ - } - - /* Convert result based on which quarter of unit circle y is in. */ - switch (n) - { - case 1: - __mp_dbl (&s, &y, p); - y = -y; - break; - - case 3: - __mp_dbl (&s, &y, p); - break; - - case 2: - __mp_dbl (&c, &y, p); - y = -y; - break; - - /* Quadrant not set, so the result must be cos (X + DX), which is also - stored in C. */ - case 0: - default: - __mp_dbl (&c, &y, p); - } - LIBC_PROBE (slowcos, 3, &x, &dx, &y); - return y; -} - -/* Perform range reduction of a double number x into multi precision number y, - such that y = x - n * pi / 2, abs (y) < pi / 4, n = 0, +-1, +-2, ... - Return int which indicates in which quarter of circle x is. */ -int -SECTION -__mpranred (double x, mp_no *y, int p) -{ - number v; - double t, xn; - int i, k, n; - mp_no a, b, c; - - if (fabs (x) < 2.8e14) - { - t = (x * hpinv.d + toint.d); - 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); - return n; - } - else - { - /* If x is very big more precision required. */ - __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); - 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] >= HALFRAD) - { - t += 1.0; - __sub (&c, &__mpone, &b, p); - __mul (&b, &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); - } -} diff --git a/sysdeps/ieee754/dbl-64/sincos32.h b/sysdeps/ieee754/dbl-64/sincos32.h deleted file mode 100644 index e112329c19..0000000000 --- a/sysdeps/ieee754/dbl-64/sincos32.h +++ /dev/null @@ -1,81 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * Written by International Business Machines Corp. - * Copyright (C) 2001-2021 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 - * the Free Software Foundation; either version 2.1 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 Lesser 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, see <https://www.gnu.org/licenses/>. - */ - -/******************************************************************/ -/* */ -/* MODULE_NAME:sincos32.h */ -/* */ -/* common data and variables prototype and definition */ -/******************************************************************/ - -#ifndef SINCOS32_H -#define SINCOS32_H - -#ifdef BIG_ENDI -static const number -/**/ hpinv = {{0x3FE45F30, 0x6DC9C883}}, /* 0.63661977236758138 */ -/**/ toint = {{0x43380000, 0x00000000}}; /* 6755399441055744 */ - -#else -#ifdef LITTLE_ENDI -static const number -/**/ hpinv = {{0x6DC9C883, 0x3FE45F30}}, /* 0.63661977236758138 */ -/**/ toint = {{0x00000000, 0x43380000}}; /* 6755399441055744 */ - -#endif -#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}}; - -static const double toverp[75] = { - 10680707.0, 7228996.0, 1387004.0, 2578385.0, 16069853.0, - 12639074.0, 9804092.0, 4427841.0, 16666979.0, 11263675.0, - 12935607.0, 2387514.0, 4345298.0, 14681673.0, 3074569.0, - 13734428.0, 16653803.0, 1880361.0, 10960616.0, 8533493.0, - 3062596.0, 8710556.0, 7349940.0, 6258241.0, 3772886.0, - 3769171.0, 3798172.0, 8675211.0, 12450088.0, 3874808.0, - 9961438.0, 366607.0, 15675153.0, 9132554.0, 7151469.0, - 3571407.0, 2607881.0, 12013382.0, 4155038.0, 6285869.0, - 7677882.0, 13102053.0, 15825725.0, 473591.0, 9065106.0, - 15363067.0, 6271263.0, 9264392.0, 5636912.0, 4652155.0, - 7056368.0, 13614112.0, 10155062.0, 1944035.0, 9527646.0, - 15080200.0, 6658437.0, 6231200.0, 6832269.0, 16767104.0, - 5075751.0, 3212806.0, 1398474.0, 7579849.0, 6349435.0, - 12618859.0, 4703257.0, 12806093.0, 14477321.0, 2786137.0, - 12875403.0, 9837734.0, 14528324.0, 13719321.0, 343717.0 }; - -#endif |