diff options
Diffstat (limited to 'REORG.TODO/sysdeps/ieee754/dbl-64/sincos32.c')
| -rw-r--r-- | REORG.TODO/sysdeps/ieee754/dbl-64/sincos32.c | 369 |
1 files changed, 369 insertions, 0 deletions
diff --git a/REORG.TODO/sysdeps/ieee754/dbl-64/sincos32.c b/REORG.TODO/sysdeps/ieee754/dbl-64/sincos32.c new file mode 100644 index 0000000000..9cd8e2f97f --- /dev/null +++ b/REORG.TODO/sysdeps/ieee754/dbl-64/sincos32.c @@ -0,0 +1,369 @@ +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001-2017 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 <http://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); +} + +/* Receive double x and two double results of sin(x) and return result which is + more accurate, computing sin(x) with multi precision routine c32. */ +double +SECTION +__sin32 (double x, double res, double res1) +{ + int p; + mp_no a, b, c; + p = 32; + __dbl_mp (res, &a, p); + __dbl_mp (0.5 * (res1 - res), &b, p); + __add (&a, &b, &c, p); + if (x > 0.8) + { + __sub (&hp, &c, &a, p); + __c32 (&a, &b, &c, p); + } + else + __c32 (&c, &a, &b, p); /* b=sin(0.5*(res+res1)) */ + __dbl_mp (x, &c, p); /* c = x */ + __sub (&b, &c, &a, p); + /* if a > 0 return min (res, res1), otherwise return max (res, res1). */ + if ((a.d[0] > 0 && res >= res1) || (a.d[0] <= 0 && res <= res1)) + res = res1; + LIBC_PROBE (slowasin, 2, &res, &x); + return res; +} + +/* Receive double x and two double results of cos(x) and return result which is + more accurate, computing cos(x) with multi precision routine c32. */ +double +SECTION +__cos32 (double x, double res, double res1) +{ + int p; + mp_no a, b, c; + p = 32; + __dbl_mp (res, &a, p); + __dbl_mp (0.5 * (res1 - res), &b, p); + __add (&a, &b, &c, p); + if (x > 2.4) + { + __sub (&pi, &c, &a, p); + __c32 (&a, &b, &c, p); + b.d[0] = -b.d[0]; + } + else if (x > 0.8) + { + __sub (&hp, &c, &a, p); + __c32 (&a, &c, &b, p); + } + else + __c32 (&c, &b, &a, p); /* b=cos(0.5*(res+res1)) */ + __dbl_mp (x, &c, p); /* c = x */ + __sub (&b, &c, &a, p); + /* if a > 0 return max (res, res1), otherwise return min (res, res1). */ + if ((a.d[0] > 0 && res <= res1) || (a.d[0] <= 0 && res >= res1)) + res = res1; + LIBC_PROBE (slowacos, 2, &res, &x); + return res; +} + +/* 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); + } +} |