/* * IBM Accurate Mathematical Library * written by International Business Machines Corp. * Copyright (C) 2001-2022 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:usncs.c */ /* */ /* FUNCTIONS: usin */ /* ucos */ /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */ /* branred.c sincos.tbl */ /* */ /* An ultimate sin and cos routine. Given an IEEE double machine number x */ /* it computes sin(x) or cos(x) with ~0.55 ULP. */ /* Assumption: Machine arithmetic operations are performed in */ /* round to nearest mode of IEEE 754 standard. */ /* */ /****************************************************************************/ #include <errno.h> #include <float.h> #include "endian.h" #include "mydefs.h" #include "usncs.h" #include <math.h> #define attribute_hidden #if !defined(__always_inline) #define __always_inline #endif /* Helper macros to compute sin of the input values. */ #define POLYNOMIAL2(xx) … #define POLYNOMIAL(xx) … /* The computed polynomial is a variation of the Taylor series expansion for sin(x): x - x^3/3! + x^5/5! - x^7/7! + x^9/9! - dx*x^2/2 + dx The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so on. The result is returned to LHS. */ #define TAYLOR_SIN(xx, x, dx) … #define SINCOS_TABLE_LOOKUP(u, sn, ssn, cs, ccs) … #ifndef SECTION #define SECTION #endif extern const union { … } __sincostab attribute_hidden; static const double sn3 = …, sn5 = …, cs2 = …, cs4 = …, cs6 = …; int __branred (double x, double *a, double *aa); /* Given a number partitioned into X and DX, this function computes the cosine of the number by combining the sin and cos of X (as computed by a variation of the Taylor series) with the values looked up from the sin/cos table to get the result. */ static __always_inline double do_cos (double x, double dx) { … } /* Given a number partitioned into X and DX, this function computes the sine of the number by combining the sin and cos of X (as computed by a variation of the Taylor series) with the values looked up from the sin/cos table to get the result. */ static __always_inline double do_sin (double x, double dx) { … } /* Reduce range of x to within PI/2 with abs (x) < 105414350. The high part is written to *a, the low part to *da. Range reduction is accurate to 136 bits so that when x is large and *a very close to zero, all 53 bits of *a are correct. */ static __always_inline int4 reduce_sincos (double x, double *a, double *da) { … } /* Compute sin or cos (A + DA) for the given quadrant N. */ static __always_inline double do_sincos (double a, double da, int4 n) { … } /*******************************************************************/ /* An ultimate sin routine. Given an IEEE double machine number x */ /* it computes the rounded value of sin(x). */ /*******************************************************************/ #ifndef IN_SINCOS double SECTION glibc_sin (double x) { … } /*******************************************************************/ /* An ultimate cos routine. Given an IEEE double machine number x */ /* it computes the rounded value of cos(x). */ /*******************************************************************/ double SECTION glibc_cos (double x) { … } #endif