AuroraMiddleware/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2024-11-10 07:10:06 +00:00
Code Issues Packages Projects Releases Wiki Activity

math: Remove mpa files (part 2) [BZ #15267]

Browse Source 
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>
This commit is contained in:
Wilco Dijkstra 2021-03-11 15:36:14 +00:00
parent 47ad14d789
commit 92cfc9ad82
15 changed files with 0 additions and 2452 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

81
sysdeps/ieee754/dbl-64/doasin.c
View File

@ -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] */
}

63
sysdeps/ieee754/dbl-64/doasin.h
View File

@ -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

217
sysdeps/ieee754/dbl-64/dosincos.c
View File

@ -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];
}
}

80
sysdeps/ieee754/dbl-64/dosincos.h
View File

@ -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

47
sysdeps/ieee754/dbl-64/mpa-arch.h
View File

@ -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))

913
sysdeps/ieee754/dbl-64/mpa.c
View File

@ -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);
}
}

123
sysdeps/ieee754/dbl-64/mpa.h
View File

@ -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);

116
sysdeps/ieee754/dbl-64/mpatan.c
View File

@ -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);
}

145
sysdeps/ieee754/dbl-64/mpatan.h
View File

@ -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

67
sysdeps/ieee754/dbl-64/mpatan2.c
View File

@ -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);
}
}

111
sysdeps/ieee754/dbl-64/mpsqrt.c
View File

@ -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);
}

38
sysdeps/ieee754/dbl-64/mpsqrt.h
View File

@ -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

63
sysdeps/ieee754/dbl-64/mptan.c
View File

@ -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);
}

307
sysdeps/ieee754/dbl-64/sincos32.c
View File

@ -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);
}
}

81
sysdeps/ieee754/dbl-64/sincos32.h
View File

@ -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