From 82edd901fee21bf934a1b53bcf44de2bd08ce6e0 Mon Sep 17 00:00:00 2001 From: Ryan Houdek Date: Sun, 30 Mar 2025 10:14:57 -0700 Subject: [PATCH] FEXCore/Common: Adds cephes math library Only the few transcendental functions that FEX needs. Disabled when building on x86-64. --- External/cephes/CMakeLists.txt | 29 ++ External/cephes/LICENSE | 118 ++++++ External/cephes/include/cephes_128bit.h | 10 + External/cephes/src/Impl.cpp | 36 ++ External/cephes/src/atanll.c | 222 +++++++++++ External/cephes/src/constll.c | 32 ++ External/cephes/src/exp2ll.c | 125 +++++++ External/cephes/src/floorll.c | 479 ++++++++++++++++++++++++ External/cephes/src/log2ll.c | 206 ++++++++++ External/cephes/src/mconf.h | 176 +++++++++ External/cephes/src/mtherr.c | 102 +++++ External/cephes/src/polevll.c | 97 +++++ External/cephes/src/sinll.c | 270 +++++++++++++ External/cephes/src/tanll.c | 214 +++++++++++ 14 files changed, 2116 insertions(+) create mode 100644 External/cephes/CMakeLists.txt create mode 100644 External/cephes/LICENSE create mode 100644 External/cephes/include/cephes_128bit.h create mode 100644 External/cephes/src/Impl.cpp create mode 100644 External/cephes/src/atanll.c create mode 100644 External/cephes/src/constll.c create mode 100644 External/cephes/src/exp2ll.c create mode 100644 External/cephes/src/floorll.c create mode 100644 External/cephes/src/log2ll.c create mode 100644 External/cephes/src/mconf.h create mode 100644 External/cephes/src/mtherr.c create mode 100644 External/cephes/src/polevll.c create mode 100644 External/cephes/src/sinll.c create mode 100644 External/cephes/src/tanll.c diff --git a/External/cephes/CMakeLists.txt b/External/cephes/CMakeLists.txt new file mode 100644 index 000000000..1024f169d --- /dev/null +++ b/External/cephes/CMakeLists.txt @@ -0,0 +1,29 @@ +set (SRCS + src/Impl.cpp + ) + +if (_M_X86_64 OR MINGW_BUILD) + set (LOW_PRECISION TRUE) +endif() + +if (NOT LOW_PRECISION) + list(APPEND SRCS + src/atanll.c + src/constll.c + src/exp2ll.c + src/floorll.c + src/log2ll.c + src/mtherr.c + src/polevll.c + src/sinll.c + src/tanll.c) +endif() + +add_library(cephes_128bit STATIC ${SRCS}) +target_include_directories(cephes_128bit PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include/) +if (NOT LOW_PRECISION) + target_compile_options(cephes_128bit PRIVATE -fno-builtin) + target_compile_definitions(cephes_128bit PRIVATE -DLOW_PRECISION=0) +else() + target_compile_definitions(cephes_128bit PRIVATE -DLOW_PRECISION=1) +endif() diff --git a/External/cephes/LICENSE b/External/cephes/LICENSE new file mode 100644 index 000000000..3cfe13ae3 --- /dev/null +++ b/External/cephes/LICENSE @@ -0,0 +1,118 @@ +The cephes math library is BSD licensed. +The source can be accessed from https://www.netlib.org/cephes/ + +Original license from https://www.netlib.org/cephes/readme : +> Some software in this archive may be from the book _Methods and +> Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster +> International, 1989) or from the Cephes Mathematical Library, a +> commercial product. In either event, it is copyrighted by the author. +> What you see here may be used freely but it comes with no support or +> guarantee. +> +> The two known misprints in the book are repaired here in the +> source listings for the gamma function and the incomplete beta +> integral. +> +> +> Stephen L. Moshier +> moshier@na-net.ornl.gov + +The author was e-mailed and they allowed it to be relicensed under BSD. +Resources: +https://bugs.gentoo.org/687276 +https://lists.debian.org/debian-legal/2004/12/msg00295.html +https://github.com/deepmind/torch-cephes/blob/master/LICENSE.txt +https://github.com/nearform/node-cephes/blob/master/LICENSE + +E-mail snippit from torch-cephes source: + +Return-Path: +X-Original-To: julien@cornebise.com +Delivered-To: julien@cornebise.com +Received: from atl4mhob11.myregisteredsite.com (atl4mhob11.myregisteredsite.com [209.17.115.49]) + by cornebise.com (Postfix) with ESMTP id D47B139FC0 + for ; Fri, 25 Oct 2013 16:32:40 +0200 (CEST) +Received: from mailpod1.hostingplatform.com ([10.30.71.116]) + by atl4mhob11.myregisteredsite.com (8.14.4/8.14.4) with ESMTP id r9PEWcwQ003543 + for ; Fri, 25 Oct 2013 10:32:38 -0400 +Received: (qmail 11948 invoked by uid 0); 25 Oct 2013 12:36:20 -0000 +X-TCPREMOTEIP: 76.24.25.74 +X-Authenticated-UID: steve@moshier.net +Received: from unknown (HELO d510.local) (steve@moshier.net@76.24.25.74) + by 0 with ESMTPA; 25 Oct 2013 12:36:20 -0000 +Date: Fri, 25 Oct 2013 08:36:19 -0400 (EDT) +From: Stephen Moshier +X-X-Sender: steve@d510 +To: Julien Cornebise +Subject: Re: Cephes: permission to wrap+distribute for Lua +In-Reply-To: <52653AD3.1010004@cornebise.com> +Message-ID: +References: <52653AD3.1010004@cornebise.com> +User-Agent: Alpine 2.02 (DEB 1266 2009-07-14) +MIME-Version: 1.0 +Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed + + +Julien, thank you for writing. +BSD license is fine, modification is OK. +There are more build scripts available in the web site distributions than +there are on the Netlib. I think there is an update to Planck's radiation +function that I haven't sent to Netlib yet. But Netlib is a more stable +site, so it is better to cite that as a reference. + + +On Mon, 21 Oct 2013, Julien Cornebise wrote: + +> -----BEGIN PGP SIGNED MESSAGE----- +> Hash: SHA1 +> +> Dear Mr Moshier +> +> I am a researcher in mathematics and machine learning in London, and +> am writing about your awesome Cephes library, whom I found at the +> heart of Scipy. +> +> It is so useful that, with your permission, I would like to wrap it +> for Lua and Torch (a machine learning overlay to Lua, specialized in +> neural nets, see http://www.torch.ch). I would like to distribute it +> as a package for Torch, including your source code along the wrapping +> code. +> This wouldbe a public package, distributed under BSD License. I have +> put a first draft on github: +> https://github.com/jucor/torch-cephes +> +> Hence my three questions, please: +> +> 1/ How would you like to be acknowledged, beyond the comments that are +> already in your code? Do you have any standard header/disclaimer that +> I could add to the documentation? +> +> 2/ At the moment, your code is left untouched. However, if I ever need +> to modify bits of the code, what are the conditions/restrictions? +> Nothing huge -- I definitely do not want to mess with it: I was +> planning to use the natural completion of some functions on the +> completed real line (e.g. CDF returing 1 when called with "infinity", +> or quantiles returning -Infinity when called with 0), either natively +> if supported, or by setting a specific flag via mtherr(). +> +> 3/ I am currently using the source from Netlib. Do you recommend using +> the source from your website instead ? +> +> Thank you very much for your attention, +> and, more importantly, for the time and effort your poured into Cephes. +> +> Best regards, +> +> Julien Cornebise, Ph.D. +> London, UK +> http://www.cornebise.com/julien +> -----BEGIN PGP SIGNATURE----- +> Version: GnuPG v1.4.14 (Darwin) +> Comment: GPGTools - http://gpgtools.org +> Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/ +> +> iEYEARECAAYFAlJlOtEACgkQKYR3gC0rw/gIpQCfZKu6+iDh9ghhm6QfsLXnldKN +> BuIAn2zZHu1c/IrRAevhjM7N7xGg0LHO +> =WeP5 +> -----END PGP SIGNATURE----- + diff --git a/External/cephes/include/cephes_128bit.h b/External/cephes/include/cephes_128bit.h new file mode 100644 index 000000000..cc666a173 --- /dev/null +++ b/External/cephes/include/cephes_128bit.h @@ -0,0 +1,10 @@ +#pragma once + +namespace FEXCore::cephes_128bit { + long double atan2l(long double y, long double x); + long double cosl(long double x); + long double exp2l(long double x); + long double log2l(long double x); + long double sinl(long double x); + long double tanl(long double x); +} diff --git a/External/cephes/src/Impl.cpp b/External/cephes/src/Impl.cpp new file mode 100644 index 000000000..c5fa76011 --- /dev/null +++ b/External/cephes/src/Impl.cpp @@ -0,0 +1,36 @@ +#include "cephes_128bit.h" + +#if !LOW_PRECISION +extern "C" { +// cephes_128bit functions +long double atan2l(long double y, long double x); +long double cosl(long double x); +long double exp2l(long double x); +long double log2l(long double x); +long double sinl(long double x); +long double tanl(long double x); +} +#else +#include +#endif + +namespace FEXCore::cephes_128bit { + long double atan2l(long double y, long double x) { + return ::atan2l(y, x); + } + long double cosl(long double x) { + return ::cosl(x); + } + long double exp2l(long double x) { + return ::exp2l(x); + } + long double log2l(long double x) { + return ::log2l(x); + } + long double sinl(long double x) { + return ::sinl(x); + } + long double tanl(long double x) { + return ::tanl(x); + } +} diff --git a/External/cephes/src/atanll.c b/External/cephes/src/atanll.c new file mode 100644 index 000000000..15f889ab4 --- /dev/null +++ b/External/cephes/src/atanll.c @@ -0,0 +1,222 @@ +/* atanl.c + * + * Inverse circular tangent, 128-bit long double precision + * (arctangent) + * + * + * + * SYNOPSIS: + * + * long double x, y, atanl(); + * + * y = atanl( x ); + * + * + * + * DESCRIPTION: + * + * Returns radian angle between -pi/2 and +pi/2 whose tangent + * is x. + * + * Range reduction is from four intervals into the interval + * from zero to tan( pi/8 ). The approximant uses a rational + * function of degree 3/4 of the form x + x**3 P(x)/Q(x). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10, 10 100,000 2.6e-34 6.5e-35 + * + */ + /* atan2l() + * + * Quadrant correct inverse circular tangent, + * long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, z, atan2l(); + * + * z = atan2l( y, x ); + * + * + * + * DESCRIPTION: + * + * Returns radian angle whose tangent is y/x. + * Define compile time symbol ANSIC = 1 for ANSI standard, + * range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range + * 0 to 2PI, args (x,y). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10, 10 100,000 3.2e-34 5.9e-35 + * See atan.c. + * + */ + +/* atan.c */ + + +/* +Cephes Math Library Release 2.2: December, 1990 +Copyright 1984, 1990 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + + +#include "mconf.h" + +/* arctan(x) = x + x^3 P(x^2) + * Theoretical peak relative error = 3.0e-36 + * relative peak error spread = 6.6e-8 + */ +static long double P[9] = { +-6.635810778635296712545011270011752799963E-4L, +-8.768423468036849091777415076702113400070E-1L, +-2.548067867495502632615671450650071218995E1L, +-2.497759878476618348858065206895055957104E2L, +-1.148164399808514330375280133523543970854E3L, +-2.792272753241044941703278827346430350236E3L, +-3.696264445691821235400930243493001671932E3L, +-2.514829758941713674909996882101723647996E3L, +-6.880597774405940432145577545328795037141E2L +}; +static long double Q[8] = { +/* 1.000000000000000000000000000000000000000E0L, */ + 3.566239794444800849656497338030115886153E1L, + 4.308348370818927353321556740027020068897E2L, + 2.494680540950601626662048893678584497900E3L, + 7.928572347062145288093560392463784743935E3L, + 1.458510242529987155225086911411015961174E4L, + 1.547394317752562611786521896296215170819E4L, + 8.782996876218210302516194604424986107121E3L, + 2.064179332321782129643673263598686441900E3L +}; + +/* tan( 3*pi/8 ) */ +static long double T3P8 = 2.414213562373095048801688724209698078569672L; + +/* tan( pi/8 ) */ +static long double TP8 = 0.414213562373095048801688724209698078569672L; + + + + +long double atanl(x) +long double x; +{ +extern long double PIO2L, PIO4L; +long double y, z; +long double polevll(), p1evll(); +short sign; + +/* make argument positive and save the sign */ +sign = 1; +if( x < 0.0L ) + { + sign = -1; + x = -x; + } + +/* range reduction */ +if( x > T3P8 ) + { + y = PIO2L; + x = -( 1.0L/x ); + } + +else if( x > TP8 ) + { + y = PIO4L; + x = (x-1.0L)/(x+1.0L); + } +else + y = 0.0L; + +/* rational form in x**2 */ +z = x * x; +y = y + ( polevll( z, P, 8 ) / p1evll( z, Q, 8 ) ) * z * x + x; + +if( sign < 0 ) + y = -y; + +return(y); +} + +/* atan2 */ + + +extern long double PIL, PIO2L; + +#if ANSIC +long double atan2l( y, x ) +#else +long double atan2l( x, y ) +#endif +long double x, y; +{ +long double z, w; +short code; +long double atanl(); + + +code = 0; +w = 0.0L; + +if( x < 0.0L ) + code = 2; +if( y < 0.0L ) + code |= 1; + +if( x == 0.0L ) + { + if( code & 1 ) + { +#if ANSIC + return( -PIO2L ); +#else + return( 3.0L*PIO2L ); +#endif + } + if( y == 0.0L ) + return( 0.0L ); + return( PIO2L ); + } + +if( y == 0.0L ) + { + if( code & 2 ) + return( PIL ); + return( 0.0L ); + } + + +switch( code ) + { +#if ANSIC + case 0: + case 1: w = 0.0L; break; + case 2: w = PIL; break; + case 3: w = -PIL; break; +#else + case 0: w = 0.0L; break; + case 1: w = 2.0L * PIL; break; + case 2: + case 3: w = PIL; break; +#endif + } + +z = atanl( y/x ); + +return( w + z ); +} diff --git a/External/cephes/src/constll.c b/External/cephes/src/constll.c new file mode 100644 index 000000000..4a89a4151 --- /dev/null +++ b/External/cephes/src/constll.c @@ -0,0 +1,32 @@ +/* (1 - 2^-113) 2^16384 */ +long double MAXNUML = 1.189731495357231765085759326628007016196469e4932L; + +/* 2^-113 */ +long double MACHEPL = 9.629649721936179265279889712924636592690508e-35L; + +/* (1 + 2^-112) 2^-16382 */ +long double UFTHRESHL = 3.362103143112093506262677817321753250115591e-4932L; + +/* 2^-16494 */ +long double MINNUML = 6.475175119438025110924438958227646552499569e-4966L; + +/* ln(MAXNUM) */ +long double MAXLOGL = 1.1356523406294143949491931077970764891253E4L; + +/* ln(MINNUM) */ +long double MINLOGL = -1.143276959615573793352782661133116431383730e4L; + +/* ln(UFTHRESH) */ +/* long double MINLOGL = -1.135513711193302405887309661372784853802025e4L; */ + +long double PIL = 3.141592653589793238462643383279502884197169L; + +long double PIO2L = 1.570796326794896619231321691639751442098585L; + +long double PIO4L = 0.7853981633974483096156608458198757210492923L; + +long double LOGE2L = 0.6931471805599453094172321214581765680755001L; + +long double LOG2EL = 1.442695040888963407359924681001892137426646L; + +long double INFINITYL = 1.0L / 0.0L; diff --git a/External/cephes/src/exp2ll.c b/External/cephes/src/exp2ll.c new file mode 100644 index 000000000..76001bb16 --- /dev/null +++ b/External/cephes/src/exp2ll.c @@ -0,0 +1,125 @@ +/* exp2l.c + * + * Base 2 exponential function, 128-bit long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, exp2l(); + * + * y = exp2l( x ); + * + * + * + * DESCRIPTION: + * + * Returns 2 raised to the x power. + * + * Range reduction is accomplished by separating the argument + * into an integer k and fraction f such that + * x k f + * 2 = 2 2. + * + * A Pade' form + * + * 1 + 2x P(x**2) / (Q(x**2) - x P(x**2) ) + * + * approximates 2**x in the basic range [-0.5, 0.5]. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE +-16300 100,000 2.0e-34 4.8e-35 + * + * + * See exp.c for comments on error amplification. + * + * + * ERROR MESSAGES: + * + * message condition value returned + * exp2l underflow x < -16382 0.0 + * exp2l overflow x >= 16384 MAXNUM + * + */ + + +/* +Cephes Math Library Release 2.2: January, 1991 +Copyright 1984, 1991 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + + + +#include "mconf.h" +static char fname[] = {"exp2l"}; + +/* Pade' coefficients for 2^x - 1 + Theoretical peak relative error = 1.4e-40, + relative peak error spread = 6.8e-14 + */ +static long double P[5] = { + 1.587171580015525194694938306936721666031E2L, + 6.185032670011643762127954396427045467506E5L, + 5.677513871931844661829755443994214173883E8L, + 1.530625323728429161131811299626419117557E11L, + 9.079594442980146270952372234833529694788E12L +}; +static long double Q[5] = { +/* 1.000000000000000000000000000000000000000E0L, */ + 1.236602014442099053716561665053645270207E4L, + 2.186249607051644894762167991800811827835E7L, + 1.092141473886177435056423606755843616331E10L, + 1.490560994263653042761789432690793026977E12L, + 2.619817175234089411411070339065679229869E13L +}; + + +#define MAXL2 16384.0L +#define MINL2 -16382.0L + + +extern long double MAXNUML; + +long double exp2l(x) +long double x; +{ +long double px, xx; +int n; +long double polevll(), p1evll(), floorl(), ldexpl(); + +if( x >= MAXL2) + { + mtherr( fname, OVERFLOW ); + return( MAXNUML ); + } + +if( x < MINL2 ) + { + mtherr( fname, UNDERFLOW ); + return(0.0L); + } + +xx = x; /* save x */ +/* separate into integer and fractional parts */ +px = floorl(x+0.5L); +n = px; +x = x - px; + +/* rational approximation + * exp2(x) = 1.0 + 2xP(xx)/(Q(xx) - P(xx)) + * where xx = x**2 + */ +xx = x * x; +px = x * polevll( xx, P, 4 ); +x = px / ( p1evll( xx, Q, 5 ) - px ); +x = 1.0L + ldexpl( x, 1 ); + +/* scale by power of 2 */ +x = ldexpl( x, n ); +return(x); +} diff --git a/External/cephes/src/floorll.c b/External/cephes/src/floorll.c new file mode 100644 index 000000000..42d79e24f --- /dev/null +++ b/External/cephes/src/floorll.c @@ -0,0 +1,479 @@ +/* ceill() + * floorl() + * frexpl() + * ldexpl() + * fabsl() + * signbitl() + * isnanl() + * isfinitel() + * + * Floating point numeric utilities + * + * + * + * SYNOPSIS: + * + * long double x, y; + * long double ceill(), floorl(), frexpl(), ldexpl(), fabsl(); + * int signbitl(), isnanl(), isfinitel(); + * int expnt, n; + * + * y = floorl(x); + * y = ceill(x); + * y = frexpl( x, &expnt ); + * y = ldexpl( x, n ); + * y = fabsl( x ); + * + * + * + * DESCRIPTION: + * + * All four routines return a long double precision floating point + * result. + * + * floorl() returns the largest integer less than or equal to x. + * It truncates toward minus infinity. + * + * ceill() returns the smallest integer greater than or equal + * to x. It truncates toward plus infinity. + * + * frexpl() extracts the exponent from x. It returns an integer + * power of two to expnt and the significand between 0.5 and 1 + * to y. Thus x = y * 2**expn. + * + * ldexpl() multiplies x by 2**n. + * + * fabsl() returns the absolute value of its argument. + * + * signbitl(x) returns 1 if the sign bit of x is 1, else 0. + * + * These functions are part of the standard C run time library + * for some but not all C compilers. The ones supplied are + * written in C for IEEE arithmetic. They should + * be used only if your compiler library does not already have + * them. + * + * The IEEE versions assume that denormal numbers are implemented + * in the arithmetic. Some modifications will be required if + * the arithmetic has abrupt rather than gradual underflow. + */ + + +/* +Cephes Math Library Release 2.2: July, 1992 +Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + + +#include "mconf.h" +#define DENORMAL 1 + +#ifdef UNK +char *unkmsg = "ceill(), floorl(), frexpl(), ldexpl() must be rewritten!\n"; +#undef UNK +#define MIEEE 1 +#define EXPOFS 0 +#endif + +#ifdef IBMPC +#define NBITS 113 +#define EXPOFS 7 +#endif + +#ifdef MIEEE +#define NBITS 113 +#define EXPOFS 0 +#endif + +extern long double MAXNUML; + + +long double fabsl(x) +long double x; +{ +if( x < 0 ) + return( -x ); +else + return( x ); +} + + + +long double ceill(x) +long double x; +{ +long double y; +long double floorl(); + +#ifdef UNK +mtherr( "ceill", DOMAIN ); +return(0.0L); +#endif + +y = floorl(x); +if( y < x ) + y += 1.0L; +return(y); +} + + + + +/* Bit clearing masks: */ + +static unsigned short bmask[] = { +0xffff, +0xfffe, +0xfffc, +0xfff8, +0xfff0, +0xffe0, +0xffc0, +0xff80, +0xff00, +0xfe00, +0xfc00, +0xf800, +0xf000, +0xe000, +0xc000, +0x8000, +0x0000, +}; + + + +long double floorl(x) +long double x; +{ +union + { + long double y; + unsigned short sh[8]; + } u; +int e, j; + +#ifdef UNK +mtherr( "floor", DOMAIN ); +return(0.0L); +#endif + +u.y = x; +/* find the exponent (power of 2) */ +e = (u.sh[EXPOFS] & 0x7fff) - 0x3fff; + +if( e < 0 ) + { + if( u.y < 0 ) + return( -1.0L ); + else + return( 0.0L ); + } + +#ifdef IBMPC +j = 0; +#endif + +#ifdef MIEEE +j = 7; +#endif + +e = (NBITS - 1) - e; +/* clean out 16 bits at a time */ +while( e >= 16 ) + { +#ifdef IBMPC + u.sh[j++] = 0; +#endif + +#ifdef MIEEE + u.sh[j--] = 0; +#endif + e -= 16; + } + +/* clear the remaining bits */ +if( e > 0 ) + u.sh[j] &= bmask[e]; + +if( (x < 0.0L) && (u.y != x) ) + u.y -= 1.0L; + +return(u.y); +} + + + +long double frexpl( x, pw2 ) +long double x; +int *pw2; +{ +union + { + long double y; + unsigned short sh[8]; + } u; +int i, k; + +u.y = x; + +#ifdef UNK +mtherr( "frexp", DOMAIN ); +return(0.0L); +#endif + +/* find the exponent (power of 2) */ +i = u.sh[EXPOFS] & 0x7fff; + +if( i == 0 ) + { + if( u.y == 0.0L ) + { + *pw2 = 0; + return(0.0L); + } +/* Number is denormal or zero */ +#if DENORMAL +/* Handle denormal number. */ +do + { + u.y *= 2.0L; + i -= 1; + k = u.sh[EXPOFS] & 0x7fff; + } +while( (k == 0) && (i > -115) ); +i = i + k; +#else + *pw2 = 0; + return(0.0L); +#endif /* DENORMAL */ + } + +*pw2 = i - 0x3ffe; +u.sh[EXPOFS] = 0x3ffe; +return( u.y ); +} + + + + + + +long double ldexpl( x, pw2 ) +long double x; +int pw2; +{ +union + { + long double y; + unsigned short sh[8]; + } u; +long e; + +#ifdef UNK +mtherr( "ldexp", DOMAIN ); +return(0.0L); +#endif + +u.y = x; +while( (e = (u.sh[EXPOFS] & 0x7fffL)) == 0 ) + { +#if DENORMAL + if( u.y == 0.0L ) + { + return( 0.0L ); + } +/* Input is denormal. */ + if( pw2 > 0 ) + { + u.y *= 2.0L; + pw2 -= 1; + } + if( pw2 < 0 ) + { + if( pw2 < -113 ) + return(0.0L); + u.y *= 0.5L; + pw2 += 1; + } + if( pw2 == 0 ) + return(u.y); +#else + return( 0.0L ); +#endif + } + +e = e + pw2; + +/* Handle overflow */ +if( e > 0x7ffeL ) + { + e = u.sh[EXPOFS]; + u.y = 0.0L; + u.sh[EXPOFS] = e | 0x7fff; + return( u.y ); + } +u.sh[EXPOFS] &= 0x8000; +/* Handle denormalized results */ +if( e < 1 ) + { +#if DENORMAL + if( e < -113 ) + return(0.0L); + u.sh[EXPOFS] |= 1; + while( e < 1 ) + { + u.y *= 0.5L; + e += 1; + } + e = 0; +#else + return(0.0L); +#endif + } + +u.sh[EXPOFS] |= e & 0x7fff; +return(u.y); +} + +/* Return 1 if x is a number that is Not a Number, else return 0. */ + +int isnanl(x) +long double x; +{ +#ifdef NANS +union + { + long double d; + unsigned short s[8]; + unsigned int i[4]; + } u; + +u.d = x; + +if( sizeof(int) == 4 ) + { +#ifdef IBMPC + if( ((u.s[7] & 0x7fff) == 0x7fff) + && ((u.i[3] & 0x7fff) | u.i[2] | u.i[1] | u.i[0])) + return 1; +#endif +#ifdef MIEEE + if( ((u.i[0] & 0x7fff0000) == 0x7fff0000) + && ((u.i[0] & 0x7fff) | u.i[1] | u.i[2] | u.i[3])) + return 1; +#endif + return(0); + } +else + { /* size int not 4 */ +#ifdef IBMPC + if( (u.s[7] & 0x7fff) == 0x7fff) + { + if((u.s[6] & 0x7fff) | u.s[5] | u.s[4] | u.s[3] | u.s[2] | u.s[1] | u.s[0]) + return(1); + } +#endif +#ifdef MIEEE + if( (u.s[0] & 0x7fff) == 0x7fff) + { + if((u.s[1] & 0x7fff) | (u.s[2] & 0x7fff) | u.s[3] | u.s[4] | u.s[5] | u.s[6] | u.s[7]) + return(1); + } +#endif + return(0); + } /* size int not 4 */ + +#else +/* No NANS. */ +return(0); +#endif +} + + +/* Return 1 if x is not infinite and is not a NaN. */ + +int isfinitel(x) +long double x; +{ +#ifdef INFINITIES +union + { + long double d; + unsigned short s[8]; + unsigned int i[4]; + } u; + +u.d = x; + +if( sizeof(int) == 4 ) + { +#ifdef IBMPC + if( (u.s[7] & 0x7fff) != 0x7fff) + return 1; +#endif +#ifdef MIEEE + if( (u.i[0] & 0x7fff0000) != 0x7fff0000) + return 1; +#endif + return(0); + } +else + { +#ifdef IBMPC + if( (u.s[7] & 0x7fff) != 0x7fff) + return 1; +#endif +#ifdef MIEEE + if( (u.s[0] & 0x7fff) != 0x7fff) + return 1; +#endif + return(0); + } +#else +/* No INFINITY. */ +return(1); +#endif +} + + +/* Return 1 if the sign bit of x is 1, else 0. */ + +int signbitl(x) +long double x; +{ +union + { + long double d; + short s[8]; + int i[4]; + } u; + +u.d = x; + +if( sizeof(int) == 4 ) + { +#ifdef IBMPC + return( u.s[7] < 0 ); +#endif +#ifdef DEC +error no such DEC format +#endif +#ifdef MIEEE + return( u.i[0] < 0 ); +#endif + } +else + { +#ifdef IBMPC + return( u.s[7] < 0 ); +#endif +#ifdef DEC +error no such DEC format +#endif +#ifdef MIEEE + return( u.s[0] < 0 ); +#endif + } +} diff --git a/External/cephes/src/log2ll.c b/External/cephes/src/log2ll.c new file mode 100644 index 000000000..2a0c006bb --- /dev/null +++ b/External/cephes/src/log2ll.c @@ -0,0 +1,206 @@ +/* log2l.c + * + * Base 2 logarithm, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, log2l(); + * + * y = log2l( x ); + * + * + * + * DESCRIPTION: + * + * Returns the base 2 logarithm of x. + * + * The argument is separated into its exponent and fractional + * parts. If the exponent is between -1 and +1, the (natural) + * logarithm of the fraction is approximated by + * + * log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x). + * + * Otherwise, setting z = 2(x-1)/x+1), + * + * log(x) = z + z**3 P(z)/Q(z). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0.5, 2.0 100,000 1.3e-34 4.5e-35 + * IEEE exp(+-10000) 100,000 9.6e-35 4.0e-35 + * + * In the tests over the interval exp(+-10000), the logarithms + * of the random arguments were uniformly distributed over + * [-10000, +10000]. + * + * ERROR MESSAGES: + * + * log singularity: x = 0; returns MINLOG + * log domain: x < 0; returns MINLOG + */ + +/* +Cephes Math Library Release 2.2: January, 1991 +Copyright 1984, 1991 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +#include "mconf.h" +static char fname[] = {"log2l"}; + +/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x) + * 1/sqrt(2) <= x < sqrt(2) + * Theoretical peak relative error = 5.3e-37, + * relative peak error spread = 2.3e-14 + */ +static long double P[13] = { + 1.538612243596254322971797716843006400388E-6L, + 4.998469661968096229986658302195402690910E-1L, + 2.321125933898420063925789532045674660756E1L, + 4.114517881637811823002128927449878962058E2L, + 3.824952356185897735160588078446136783779E3L, + 2.128857716871515081352991964243375186031E4L, + 7.594356839258970405033155585486712125861E4L, + 1.797628303815655343403735250238293741397E5L, + 2.854829159639697837788887080758954924001E5L, + 3.007007295140399532324943111654767187848E5L, + 2.014652742082537582487669938141683759923E5L, + 7.771154681358524243729929227226708890930E4L, + 1.313572404063446165910279910527789794488E4L +}; +static long double Q[12] = { +/* 1.000000000000000000000000000000000000000E0L, */ + 4.839208193348159620282142911143429644326E1L, + 9.104928120962988414618126155557301584078E2L, + 9.147150349299596453976674231612674085381E3L, + 5.605842085972455027590989944010492125825E4L, + 2.248234257620569139969141618556349415120E5L, + 6.132189329546557743179177159925690841200E5L, + 1.158019977462989115839826904108208787040E6L, + 1.514882452993549494932585972882995548426E6L, + 1.347518538384329112529391120390701166528E6L, + 7.777690340007566932935753241556479363645E5L, + 2.626900195321832660448791748036714883242E5L, + 3.940717212190338497730839731583397586124E4L +}; + +/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), + * where z = 2(x-1)/(x+1) + * 1/sqrt(2) <= x < sqrt(2) + * Theoretical peak relative error = 1.1e-35, + * relative peak error spread 1.1e-9 + */ +static long double R[6] = { +-8.828896441624934385266096344596648080902E-1L, + 8.057002716646055371965756206836056074715E1L, +-2.024301798136027039250415126250455056397E3L, + 2.048819892795278657810231591630928516206E4L, +-8.977257995689735303686582344659576526998E4L, + 1.418134209872192732479751274970992665513E5L +}; +static long double S[6] = { +/* 1.000000000000000000000000000000000000000E0L, */ +-1.186359407982897997337150403816839480438E2L, + 3.998526750980007367835804959888064681098E3L, +-5.748542087379434595104154610899551484314E4L, + 4.001557694070773974936904547424676279307E5L, +-1.332535117259762928288745111081235577029E6L, + 1.701761051846631278975701529965589676574E6L +}; +/* log2(e) - 1 */ +#define LOG2EA 4.4269504088896340735992468100189213742664595E-1L + +#define SQRTH 7.071067811865475244008443621048490392848359E-1L + +extern long double MINLOGL; +long double frexpl(), ldexpl(), polevll(), p1evll(); + + +long double log2l(x) +long double x; +{ +VOLATILE long double z; +long double y; +int e; + +/* Test for domain */ +if( x <= 0.0L ) + { + if( x == 0.0L ) + mtherr( fname, SING ); + else + mtherr( fname, DOMAIN ); + return( -16384.0L ); + } + +/* separate mantissa from exponent */ + +/* Note, frexp is used so that denormal numbers + * will be handled properly. + */ +x = frexpl( x, &e ); + + +/* logarithm using log(x) = z + z**3 P(z)/Q(z), + * where z = 2(x-1)/x+1) + */ +if( (e > 2) || (e < -2) ) +{ +if( x < SQRTH ) + { /* 2( 2x-1 )/( 2x+1 ) */ + e -= 1; + z = x - 0.5L; + y = 0.5L * z + 0.5L; + } +else + { /* 2 (x-1)/(x+1) */ + z = x - 0.5L; + z -= 0.5L; + y = 0.5L * x + 0.5L; + } +x = z / y; +z = x*x; +y = x * ( z * polevll( z, R, 5 ) / p1evll( z, S, 6 ) ); +goto done; +} + + +/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ + +if( x < SQRTH ) + { + e -= 1; + x = ldexpl( x, 1 ) - 1.0L; /* 2x - 1 */ + } +else + { + x = x - 1.0L; + } +z = x*x; +y = x * ( z * polevll( x, P, 12 ) / p1evll( x, Q, 12 ) ); +y = y - ldexpl( z, -1 ); /* -0.5x^2 + ... */ + +done: + +/* Multiply log of fraction by log2(e) + * and base 2 exponent by 1 + * + * ***CAUTION*** + * + * This sequence of operations is critical and it may + * be horribly defeated by some compiler optimizers. + */ +z = y * LOG2EA; +z += x * LOG2EA; +z += y; +z += x; +z += e; +return( z ); +} + diff --git a/External/cephes/src/mconf.h b/External/cephes/src/mconf.h new file mode 100644 index 000000000..dc215fcea --- /dev/null +++ b/External/cephes/src/mconf.h @@ -0,0 +1,176 @@ +/* mconf.h + * + * Common include file for math routines + * + * + * + * SYNOPSIS: + * + * #include "mconf.h" + * + * + * + * DESCRIPTION: + * + * This file contains definitions for error codes that are + * passed to the common error handling routine mtherr() + * (which see). + * + * The file also includes a conditional assembly definition + * for the type of computer arithmetic (IEEE, DEC, Motorola + * IEEE, or UNKnown). + * + * For Digital Equipment PDP-11 and VAX computers, certain + * IBM systems, and others that use numbers with a 56-bit + * significand, the symbol DEC should be defined. In this + * mode, most floating point constants are given as arrays + * of octal integers to eliminate decimal to binary conversion + * errors that might be introduced by the compiler. + * + * For little-endian computers, such as IBM PC, that follow the + * IEEE Standard for Binary Floating Point Arithmetic (ANSI/IEEE + * Std 754-1985), the symbol IBMPC should be defined. These + * numbers have 53-bit significands. In this mode, constants + * are provided as arrays of hexadecimal 16 bit integers. + * + * Big-endian IEEE format is denoted MIEEE. On some RISC + * systems such as Sun SPARC, double precision constants + * must be stored on 8-byte address boundaries. Since integer + * arrays may be aligned differently, the MIEEE configuration + * may fail on such machines. + * + * To accommodate other types of computer arithmetic, all + * constants are also provided in a normal decimal radix + * which one can hope are correctly converted to a suitable + * format by the available C language compiler. To invoke + * this mode, define the symbol UNK. + * + * An important difference among these modes is a predefined + * set of machine arithmetic constants for each. The numbers + * MACHEP (the machine roundoff error), MAXNUM (largest number + * represented), and several other parameters are preset by + * the configuration symbol. Check the file const.c to + * ensure that these values are correct for your computer. + * + * Configurations NANS, INFINITIES, MINUSZERO, and DENORMAL + * may fail on many systems. Verify that they are supposed + * to work on your computer. + */ + +/* +Cephes Math Library Release 2.3: June, 1995 +Copyright 1984, 1987, 1989, 1995 by Stephen L. Moshier +*/ + + +/* Constant definitions for math error conditions + */ + +#define DOMAIN 1 /* argument domain error */ +#define SING 2 /* argument singularity */ +#define OVERFLOW 3 /* overflow range error */ +#define UNDERFLOW 4 /* underflow range error */ +#define TLOSS 5 /* total loss of precision */ +#define PLOSS 6 /* partial loss of precision */ + +#define EDOM 33 +#define ERANGE 34 + +/* Complex numeral. */ +typedef struct + { + double r; + double i; + } cmplx; + +typedef struct + { + float r; + float i; + } cmplxf; + +/* Long double complex numeral. */ + +typedef struct + { + long double r; + long double i; + } cmplxl; + + +/* Type of computer arithmetic */ + +/* PDP-11, Pro350, VAX: + */ +/* #define DEC 1 */ + +/* Intel IEEE, low order words come first: + */ +#define IBMPC 1 + +/* Motorola IEEE, high order words come first + * (Sun 680x0 workstation): + */ +/* #define MIEEE 1 */ + +/* UNKnown arithmetic, invokes coefficients given in + * normal decimal format. Beware of range boundary + * problems (MACHEP, MAXLOG, etc. in const.c) and + * roundoff problems in pow.c: + * (Sun SPARCstation) + */ +/* #define UNK 1 */ + +/* If you define UNK, then be sure to set BIGENDIAN properly. */ +/* #define BIGENDIAN 1 */ + +/* Define this `volatile' if your compiler thinks + * that floating point arithmetic obeys the associative + * and distributive laws. It will defeat some optimizations + * (but probably not enough of them). + * + * #define VOLATILE volatile + */ +#define VOLATILE + +/* For 12-byte long doubles on an i386, pad a 16-bit short 0 + * to the end of real constants initialized by integer arrays. + * + * #define XPD 0, + * + * Otherwise, the type is 10 bytes long and XPD should be + * defined blank (e.g., Microsoft C). + * + * #define XPD + */ +#define XPD 0, + +/* Define to support tiny denormal numbers, else undefine. */ +#define DENORMAL 1 + +/* Define to ask for infinity support, else undefine. */ +#define INFINITIES 1 + +/* Define to ask for support of numbers that are Not-a-Number, + else undefine. This may automatically define INFINITIES in some files. */ +#define NANS 1 + +/* Define to distinguish between -0.0 and +0.0. */ +#define MINUSZERO 1 + +/* Define 1 for ANSI C atan2() function + and ANSI prototypes for float arguments. + See atan.c and clog.c. */ +#define ANSIC 1 + +/* Get ANSI function prototypes, if you want them. */ +#ifdef __STDC__ +#define ANSIPROT +/* #include "protos.h" */ +int mtherr (char *, int); +#else +int mtherr(); +#endif + +/* Variable for error reporting. See mtherr.c. */ +extern int merror; diff --git a/External/cephes/src/mtherr.c b/External/cephes/src/mtherr.c new file mode 100644 index 000000000..0d94cd14d --- /dev/null +++ b/External/cephes/src/mtherr.c @@ -0,0 +1,102 @@ +/* mtherr.c + * + * Library common error handling routine + * + * + * + * SYNOPSIS: + * + * char *fctnam; + * int code; + * int mtherr(); + * + * mtherr( fctnam, code ); + * + * + * + * DESCRIPTION: + * + * This routine may be called to report one of the following + * error conditions (in the include file mconf.h). + * + * Mnemonic Value Significance + * + * DOMAIN 1 argument domain error + * SING 2 function singularity + * OVERFLOW 3 overflow range error + * UNDERFLOW 4 underflow range error + * TLOSS 5 total loss of precision + * PLOSS 6 partial loss of precision + * EDOM 33 Unix domain error code + * ERANGE 34 Unix range error code + * + * The default version of the file prints the function name, + * passed to it by the pointer fctnam, followed by the + * error condition. The display is directed to the standard + * output device. The routine then returns to the calling + * program. Users may wish to modify the program to abort by + * calling exit() under severe error conditions such as domain + * errors. + * + * Since all error conditions pass control to this function, + * the display may be easily changed, eliminated, or directed + * to an error logging device. + * + * SEE ALSO: + * + * mconf.h + * + */ + +/* +Cephes Math Library Release 2.0: April, 1987 +Copyright 1984, 1987 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +#include +#include "mconf.h" + +int merror = 0; + +/* Notice: the order of appearance of the following + * messages is bound to the error codes defined + * in mconf.h. + */ +static char *ermsg[7] = { +"unknown", /* error code 0 */ +"domain", /* error code 1 */ +"singularity", /* et seq. */ +"overflow", +"underflow", +"total loss of precision", +"partial loss of precision" +}; + + +int mtherr( name, code ) +char *name; +int code; +{ + +/* Display string passed by calling program, + * which is supposed to be the name of the + * function in which the error occurred: + */ +printf( "\n%s ", name ); + +/* Set global error message word */ +merror = code; + +/* Display error message defined + * by the code argument. + */ +if( (code <= 0) || (code >= 7) ) + code = 0; +printf( "%s error\n", ermsg[code] ); + +/* Return to calling + * program + */ +return( 0 ); +} diff --git a/External/cephes/src/polevll.c b/External/cephes/src/polevll.c new file mode 100644 index 000000000..58588c0c8 --- /dev/null +++ b/External/cephes/src/polevll.c @@ -0,0 +1,97 @@ +/* polevll.c + * p1evll.c + * + * Evaluate polynomial + * + * + * + * SYNOPSIS: + * + * int N; + * long double x, y, coef[N+1], polevl[]; + * + * y = polevll( x, coef, N ); + * + * + * + * DESCRIPTION: + * + * Evaluates polynomial of degree N: + * + * 2 N + * y = C + C x + C x +...+ C x + * 0 1 2 N + * + * Coefficients are stored in reverse order: + * + * coef[0] = C , ..., coef[N] = C . + * N 0 + * + * The function p1evll() assumes that coef[N] = 1.0 and is + * omitted from the array. Its calling arguments are + * otherwise the same as polevll(). + * + * + * SPEED: + * + * In the interest of speed, there are no checks for out + * of bounds arithmetic. This routine is used by most of + * the functions in the library. Depending on available + * equipment features, the user may wish to rewrite the + * program in microcode or assembly language. + * + */ + + +/* +Cephes Math Library Release 2.2: July, 1992 +Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ +#include "mconf.h" + + +/* Polynomial evaluator: + * P[0] x^n + P[1] x^(n-1) + ... + P[n] + */ +long double polevll( x, PP, n ) +long double x; +void *PP; +int n; +{ +register long double y; +long double *P; + +P = (long double *) PP; +y = *P++; +do + { + y = y * x + *P++; + } +while( --n ); +return(y); +} + + + +/* Polynomial evaluator: + * x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n] + */ +long double p1evll( x, PP, n ) +long double x; +void *PP; +int n; +{ +register long double y; +long double *P; + +P = (long double *) PP; +n -= 1; +y = x + *P++; +do + { + y = y * x + *P++; + } +while( --n ); +return( y ); +} diff --git a/External/cephes/src/sinll.c b/External/cephes/src/sinll.c new file mode 100644 index 000000000..e61bbc428 --- /dev/null +++ b/External/cephes/src/sinll.c @@ -0,0 +1,270 @@ +/* sinl.c + * + * Circular sine, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, sinl(); + * + * y = sinl( x ); + * + * + * + * DESCRIPTION: + * + * Range reduction is into intervals of pi/4. The reduction + * error is nearly eliminated by contriving an extended precision + * modular arithmetic. + * + * Two polynomial approximating functions are employed. + * Between 0 and pi/4 the sine is approximated by the Cody + * and Waite polynomial form + * x + x^3 P(x^2) . + * Between pi/4 and pi/2 the cosine is represented as + * 1 - .5 x^2 + x^4 Q(x^2) . + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE +-3.6e16 100,000 2.0e-34 5.3e-35 + * + * ERROR MESSAGES: + * + * message condition value returned + * sin total loss x > 2^55 0.0 + * + */ + /* cosl.c + * + * Circular cosine, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, cosl(); + * + * y = cosl( x ); + * + * + * + * DESCRIPTION: + * + * Range reduction is into intervals of pi/4. The reduction + * error is nearly eliminated by contriving an extended precision + * modular arithmetic. + * + * Two polynomial approximating functions are employed. + * Between 0 and pi/4 the cosine is approximated by + * 1 - .5 x^2 + x^4 Q(x^2) . + * Between pi/4 and pi/2 the sine is represented by the Cody + * and Waite polynomial form + * x + x^3 P(x^2) . + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE +-3.6e16 100,000 2.0e-34 5.2e-35 + * + * ERROR MESSAGES: + * + * message condition value returned + * cos total loss x > 2^55 0.0 + */ + +/* sin.c */ + +/* +Cephes Math Library Release 2.2: December, 1990 +Copyright 1985, 1990 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +#include "mconf.h" + +/* sin(x) = x + x^3 P(x^2) + * Theoretical peak relative error = 5.6e-39 + * relative peak error spread = 1.7e-9 + */ + +static long double sincof[12] = { + 6.410290407010279602425714995528976754871E-26L, +-3.868105354403065333804959405965295962871E-23L, + 1.957294039628045847156851410307133941611E-20L, +-8.220635246181818130416407184286068307901E-18L, + 2.811457254345322887443598804951004537784E-15L, +-7.647163731819815869711749952353081768709E-13L, + 1.605904383682161459812515654720205050216E-10L, +-2.505210838544171877505034150892770940116E-8L, + 2.755731922398589065255731765498970284004E-6L, +-1.984126984126984126984126984045294307281E-4L, + 8.333333333333333333333333333333119885283E-3L, +-1.666666666666666666666666666666666647199E-1L +}; +/* cos(x) = 1 - .5 x^2 + x^2 (x^2 P(x^2)) + * Theoretical peak relative error = 2.1e-37, + * relative peak error spread = 1.4e-8 + */ +static long double coscof[11] = { + 1.601961934248327059668321782499768648351E-24L, +-8.896621117922334603659240022184527001401E-22L, + 4.110317451243694098169570731967589555498E-19L, +-1.561920696747074515985647487260202922160E-16L, + 4.779477332386900932514186378501779328195E-14L, +-1.147074559772972328629102981460088437917E-11L, + 2.087675698786809897637922200570559726116E-9L, +-2.755731922398589065255365968070684102298E-7L, + 2.480158730158730158730158440896461945271E-5L, +-1.388888888888888888888888888765724370132E-3L, + 4.166666666666666666666666666666459301466E-2L +}; +/* +static long double DP1 = 7.853981554508209228515625E-1L; +static long double DP2 = 7.94662735614792836713604629039764404296875E-9L; +static long double DP3 = 3.0616169978683829430651648306875026455243736148E-17L; +static long double lossth = 5.49755813888e11L; +*/ +static long double DP1 = + 7.853981633974483067550664827649598009884357452392578125E-1L; +static long double DP2 = + 2.8605943630549158983813312792950660807511260829685741796657E-18L; +static long double DP3 = + 2.1679525325309452561992610065108379921905808E-35L; + +static long double lossth = 3.6028797018963968E16L; /* 2^55 */ + +extern long double PIO4L; + + +long double sinl(x) +long double x; +{ +long double y, z, zz; +int j, sign; +long double polevll(), floorl(), ldexpl(); + +/* make argument positive but save the sign */ +sign = 1; +if( x < 0 ) + { + x = -x; + sign = -1; + } + +if( x > lossth ) + { + mtherr( "sinl", TLOSS ); + return(0.0L); + } + +y = floorl( x/PIO4L ); /* integer part of x/PIO4 */ + +/* strip high bits of integer part to prevent integer overflow */ +z = ldexpl( y, -4 ); +z = floorl(z); /* integer part of y/8 */ +z = y - ldexpl( z, 4 ); /* y - 16 * (y/16) */ + +j = z; /* convert to integer for tests on the phase angle */ +/* map zeros to origin */ +if( j & 1 ) + { + j += 1; + y += 1.0L; + } +j = j & 07; /* octant modulo 360 degrees */ +/* reflect in x axis */ +if( j > 3) + { + sign = -sign; + j -= 4; + } + +/* Extended precision modular arithmetic */ +z = ((x - y * DP1) - y * DP2) - y * DP3; + +zz = z * z; +if( (j==1) || (j==2) ) + { + y = 1.0L - ldexpl(zz,-1) + zz * zz * polevll( zz, coscof, 10 ); + } +else + { + y = z + z * (zz * polevll( zz, sincof, 11 )); + } + +if(sign < 0) + y = -y; + +return(y); +} + + + + + +long double cosl(x) +long double x; +{ +long double y, z, zz; +long i; +int j, sign; +long double polevll(), floorl(), ldexpl(); + + +/* make argument positive */ +sign = 1; +if( x < 0 ) + x = -x; + +if( x > lossth ) + { + mtherr( "cosl", TLOSS ); + return(0.0L); + } + +y = floorl( x/PIO4L ); +z = ldexpl( y, -4 ); +z = floorl(z); /* integer part of y/8 */ +z = y - ldexpl( z, 4 ); /* y - 16 * (y/16) */ + +/* integer and fractional part modulo one octant */ +i = z; +if( i & 1 ) /* map zeros to origin */ + { + i += 1; + y += 1.0L; + } +j = i & 07; +if( j > 3) + { + j -=4; + sign = -sign; + } + +if( j > 1 ) + sign = -sign; + +/* Extended precision modular arithmetic */ +z = ((x - y * DP1) - y * DP2) - y * DP3; + +zz = z * z; +if( (j==1) || (j==2) ) + { + y = z + z * (zz * polevll( zz, sincof, 11 )); + } +else + { + y = 1.0L - ldexpl(zz,-1) + zz * zz * polevll( zz, coscof, 10 ); + } + +if(sign < 0) + y = -y; + +return(y); +} diff --git a/External/cephes/src/tanll.c b/External/cephes/src/tanll.c new file mode 100644 index 000000000..cf79ca913 --- /dev/null +++ b/External/cephes/src/tanll.c @@ -0,0 +1,214 @@ +/* tanl.c + * + * Circular tangent, 128-bit long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, tanl(); + * + * y = tanl( x ); + * + * + * + * DESCRIPTION: + * + * Returns the circular tangent of the radian argument x. + * + * Range reduction is modulo pi/4. A rational function + * x + x**3 P(x**2)/Q(x**2) + * is employed in the basic interval [0, pi/4]. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE +-3.6e16 100,000 3.0e-34 7.2e-35 + * + * ERROR MESSAGES: + * + * message condition value returned + * tan total loss x > 2^55 0.0 + * + */ + /* cotl.c + * + * Circular cotangent, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, cotl(); + * + * y = cotl( x ); + * + * + * + * DESCRIPTION: + * + * Returns the circular cotangent of the radian argument x. + * + * Range reduction is modulo pi/4. A rational function + * x + x**3 P(x**2)/Q(x**2) + * is employed in the basic interval [0, pi/4]. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE +-3.6e16 100,000 2.9e-34 7.2e-35 + * + * + * ERROR MESSAGES: + * + * message condition value returned + * cot total loss x > 2^55 0.0 + * cot singularity x = 0 MAXNUM + * + */ + +/* +Cephes Math Library Release 2.2: December, 1990 +Copyright 1984, 1990 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +#include "mconf.h" + +/* tan(x) = x + x^3 P(x^2) + * 0 <= |x| <= pi/4 + * Theoretical peak relative error = 4.3e-38 + * relative peak error spread = 6.1e-11 + */ +static long double P[6] = { +-9.889929415807650724957118893791829849557E-1L, + 1.272297782199996882828849455156962260810E3L, +-4.249691853501233575668486667664718192660E5L, + 5.160188250214037865511600561074819366815E7L, +-2.307030822693734879744223131873392503321E9L, + 2.883414728874239697964612246732416606301E10L +}; +static long double Q[6] = { +/* 1.000000000000000000000000000000000000000E0L, */ +-1.317243702830553658702531997959756728291E3L, + 4.529422062441341616231663543669583527923E5L, +-5.733709132766856723608447733926138506824E7L, + 2.758476078803232151774723646710890525496E9L, +-4.152206921457208101480801635640958361612E10L, + 8.650244186622719093893836740197250197602E10L +}; + +static long double DP1 = + 7.853981633974483067550664827649598009884357452392578125E-1L; +static long double DP2 = + 2.8605943630549158983813312792950660807511260829685741796657E-18L; +static long double DP3 = + 2.1679525325309452561992610065108379921905808E-35L; + +static long double lossth = 3.6028797018963968E16L; /* 2^55 */ + +extern long double PIO4L; +extern long double MAXNUML; +static long double tancotl(); + +long double tanl(x) +long double x; +{ + +return( tancotl(x,0) ); +} + + +long double cotl(x) +long double x; +{ + +if( x == 0.0L ) + { + mtherr( "cotl", SING ); + return( MAXNUML ); + } +return( tancotl(x,1) ); +} + + +static long double tancotl( xx, cotflg ) +long double xx; +int cotflg; +{ +long double x, y, z, zz; +int j, sign; +long double polevll(), p1evll(), floorl(), ldexpl(); + +/* make argument positive but save the sign */ +if( xx < 0.0L ) + { + x = -xx; + sign = -1; + } +else + { + x = xx; + sign = 1; + } + +if( x > lossth ) + { + if( cotflg ) + mtherr( "cotl", TLOSS ); + else + mtherr( "tanl", TLOSS ); + return(0.0L); + } + +/* compute x mod PIO4 */ +y = floorl( x/PIO4L ); + +/* strip high bits of integer part */ +z = ldexpl( y, -4 ); +z = floorl(z); /* integer part of y/16 */ +z = y - ldexpl( z, 4 ); /* y - 16 * (y/16) */ + +/* integer and fractional part modulo one octant */ +j = z; + +/* map zeros and singularities to origin */ +if( j & 1 ) + { + j += 1; + y += 1.0L; + } + +z = ((x - y * DP1) - y * DP2) - y * DP3; + +zz = z * z; + +if( zz > 1.0e-20L ) + y = z + z * (zz * polevll( zz, P, 5 )/p1evll(zz, Q, 6)); +else + y = z; + +if( j & 2 ) + { + if( cotflg ) + y = -y; + else + y = -1.0L/y; + } +else + { + if( cotflg ) + y = 1.0L/y; + } + +if( sign < 0 ) + y = -y; + +return( y ); +}