Previous commit was missing deleted files in sysdeps/ieee754/dbl-64.
Finally remove all mpa related files, headers, declarations, probes, unused
tables and update makefiles.
Reviewed-By: Paul Zimmermann <Paul.Zimmermann@inria.fr>
+++ /dev/null
-/*
- * 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] */
-}
+++ /dev/null
-
-/*
- * 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
+++ /dev/null
-/*
- * 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];
- }
-}
+++ /dev/null
-
-/*
- * 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
+++ /dev/null
-/* 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))
+++ /dev/null
-/*
- * 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);
- }
-}
+++ /dev/null
-/*
- * 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);
+++ /dev/null
-/*
- * 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);
-}
+++ /dev/null
-/*
- * 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
+++ /dev/null
-/*
- * 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);
- }
-}
+++ /dev/null
-/*
- * 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);
-}
+++ /dev/null
-/*
- * 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
+++ /dev/null
-/*
- * 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);
-}
+++ /dev/null
-/*
- * 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);
- }
-}
+++ /dev/null
-/*
- * 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
git://sourceware.org / glibc.git / commitdiff