diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/sincos32.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/sincos32.c | 307 |
1 files changed, 0 insertions, 307 deletions
diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c deleted file mode 100644 index 44a313ad76..0000000000 --- a/sysdeps/ieee754/dbl-64/sincos32.c +++ /dev/null @@ -1,307 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2021 Free Software Foundation, Inc. - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU Lesser General Public License as published by - * the Free Software Foundation; either version 2.1 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with this program; if not, see <https://www.gnu.org/licenses/>. - */ -/****************************************************************/ -/* MODULE_NAME: sincos32.c */ -/* */ -/* FUNCTIONS: ss32 */ -/* cc32 */ -/* c32 */ -/* sin32 */ -/* cos32 */ -/* mpsin */ -/* mpcos */ -/* mpranred */ -/* mpsin1 */ -/* mpcos1 */ -/* */ -/* FILES NEEDED: endian.h mpa.h sincos32.h */ -/* mpa.c */ -/* */ -/* Multi Precision sin() and cos() function with p=32 for sin()*/ -/* cos() arcsin() and arccos() routines */ -/* In addition mpranred() routine performs range reduction of */ -/* a double number x into multi precision number y, */ -/* such that y=x-n*pi/2, abs(y)<pi/4, n=0,+-1,+-2,.... */ -/****************************************************************/ -#include "endian.h" -#include "mpa.h" -#include "sincos32.h" -#include <math.h> -#include <math_private.h> -#include <stap-probe.h> - -#ifndef SECTION -# define SECTION -#endif - -/* Compute Multi-Precision sin() function for given p. Receive Multi Precision - number x and result stored at y. */ -static void -SECTION -ss32 (mp_no *x, mp_no *y, int p) -{ - int i; - double a; - mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}}; - for (i = 1; i <= p; i++) - mpk.d[i] = 0; - - __sqr (x, &x2, p); - __cpy (&oofac27, &gor, p); - __cpy (&gor, &sum, p); - for (a = 27.0; a > 1.0; a -= 2.0) - { - mpk.d[1] = a * (a - 1.0); - __mul (&gor, &mpk, &mpt1, p); - __cpy (&mpt1, &gor, p); - __mul (&x2, &sum, &mpt1, p); - __sub (&gor, &mpt1, &sum, p); - } - __mul (x, &sum, y, p); -} - -/* Compute Multi-Precision cos() function for given p. Receive Multi Precision - number x and result stored at y. */ -static void -SECTION -cc32 (mp_no *x, mp_no *y, int p) -{ - int i; - double a; - mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}}; - for (i = 1; i <= p; i++) - mpk.d[i] = 0; - - __sqr (x, &x2, p); - mpk.d[1] = 27.0; - __mul (&oofac27, &mpk, &gor, p); - __cpy (&gor, &sum, p); - for (a = 26.0; a > 2.0; a -= 2.0) - { - mpk.d[1] = a * (a - 1.0); - __mul (&gor, &mpk, &mpt1, p); - __cpy (&mpt1, &gor, p); - __mul (&x2, &sum, &mpt1, p); - __sub (&gor, &mpt1, &sum, p); - } - __mul (&x2, &sum, y, p); -} - -/* Compute both sin(x), cos(x) as Multi precision numbers. */ -void -SECTION -__c32 (mp_no *x, mp_no *y, mp_no *z, int p) -{ - mp_no u, t, t1, t2, c, s; - int i; - __cpy (x, &u, p); - u.e = u.e - 1; - cc32 (&u, &c, p); - ss32 (&u, &s, p); - for (i = 0; i < 24; i++) - { - __mul (&c, &s, &t, p); - __sub (&s, &t, &t1, p); - __add (&t1, &t1, &s, p); - __sub (&__mptwo, &c, &t1, p); - __mul (&t1, &c, &t2, p); - __add (&t2, &t2, &c, p); - } - __sub (&__mpone, &c, y, p); - __cpy (&s, z, p); -} - -/* Compute sin() of double-length number (X + DX) as Multi Precision number and - return result as double. If REDUCE_RANGE is true, X is assumed to be the - original input and DX is ignored. */ -double -SECTION -__mpsin (double x, double dx, bool reduce_range) -{ - double y; - mp_no a, b, c, s; - int n; - int p = 32; - - if (reduce_range) - { - n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */ - __c32 (&a, &c, &s, p); - } - else - { - n = -1; - __dbl_mp (x, &b, p); - __dbl_mp (dx, &c, p); - __add (&b, &c, &a, p); - if (x > 0.8) - { - __sub (&hp, &a, &b, p); - __c32 (&b, &s, &c, p); - } - else - __c32 (&a, &c, &s, p); /* b = sin(x+dx) */ - } - - /* Convert result based on which quarter of unit circle y is in. */ - switch (n) - { - case 1: - __mp_dbl (&c, &y, p); - break; - - case 3: - __mp_dbl (&c, &y, p); - y = -y; - break; - - case 2: - __mp_dbl (&s, &y, p); - y = -y; - break; - - /* Quadrant not set, so the result must be sin (X + DX), which is also in - S. */ - case 0: - default: - __mp_dbl (&s, &y, p); - } - LIBC_PROBE (slowsin, 3, &x, &dx, &y); - return y; -} - -/* Compute cos() of double-length number (X + DX) as Multi Precision number and - return result as double. If REDUCE_RANGE is true, X is assumed to be the - original input and DX is ignored. */ -double -SECTION -__mpcos (double x, double dx, bool reduce_range) -{ - double y; - mp_no a, b, c, s; - int n; - int p = 32; - - if (reduce_range) - { - n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */ - __c32 (&a, &c, &s, p); - } - else - { - n = -1; - __dbl_mp (x, &b, p); - __dbl_mp (dx, &c, p); - __add (&b, &c, &a, p); - if (x > 0.8) - { - __sub (&hp, &a, &b, p); - __c32 (&b, &s, &c, p); - } - else - __c32 (&a, &c, &s, p); /* a = cos(x+dx) */ - } - - /* Convert result based on which quarter of unit circle y is in. */ - switch (n) - { - case 1: - __mp_dbl (&s, &y, p); - y = -y; - break; - - case 3: - __mp_dbl (&s, &y, p); - break; - - case 2: - __mp_dbl (&c, &y, p); - y = -y; - break; - - /* Quadrant not set, so the result must be cos (X + DX), which is also - stored in C. */ - case 0: - default: - __mp_dbl (&c, &y, p); - } - LIBC_PROBE (slowcos, 3, &x, &dx, &y); - return y; -} - -/* Perform range reduction of a double number x into multi precision number y, - such that y = x - n * pi / 2, abs (y) < pi / 4, n = 0, +-1, +-2, ... - Return int which indicates in which quarter of circle x is. */ -int -SECTION -__mpranred (double x, mp_no *y, int p) -{ - number v; - double t, xn; - int i, k, n; - mp_no a, b, c; - - if (fabs (x) < 2.8e14) - { - t = (x * hpinv.d + toint.d); - xn = t - toint.d; - v.d = t; - n = v.i[LOW_HALF] & 3; - __dbl_mp (xn, &a, p); - __mul (&a, &hp, &b, p); - __dbl_mp (x, &c, p); - __sub (&c, &b, y, p); - return n; - } - else - { - /* If x is very big more precision required. */ - __dbl_mp (x, &a, p); - a.d[0] = 1.0; - k = a.e - 5; - if (k < 0) - k = 0; - b.e = -k; - b.d[0] = 1.0; - for (i = 0; i < p; i++) - b.d[i + 1] = toverp[i + k]; - __mul (&a, &b, &c, p); - t = c.d[c.e]; - for (i = 1; i <= p - c.e; i++) - c.d[i] = c.d[i + c.e]; - for (i = p + 1 - c.e; i <= p; i++) - c.d[i] = 0; - c.e = 0; - if (c.d[1] >= HALFRAD) - { - t += 1.0; - __sub (&c, &__mpone, &b, p); - __mul (&b, &hp, y, p); - } - else - __mul (&c, &hp, y, p); - n = (int) t; - if (x < 0) - { - y->d[0] = -y->d[0]; - n = -n; - } - return (n & 3); - } -} |