mirror of
https://github.com/FEX-Emu/FEX.git
synced 2026-10-06 10:00:16 +02:00
FEXCore/Common: Adds cephes math library
Only the few transcendental functions that FEX needs. Disabled when building on x86-64.
This commit is contained in:
14 files changed
+2116
No files matched your search
Vendored
+29
@@ -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()
|
||||
Vendored
+118
@@ -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: <steve@moshier.net>
|
||||
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 <julien@cornebise.com>; 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 <julien@cornebise.com>; 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 <steve@moshier.net>
|
||||
X-X-Sender: steve@d510
|
||||
To: Julien Cornebise <julien@cornebise.com>
|
||||
Subject: Re: Cephes: permission to wrap+distribute for Lua
|
||||
In-Reply-To: <52653AD3.1010004@cornebise.com>
|
||||
Message-ID: <alpine.DEB.2.02.1310250827040.17646@d510>
|
||||
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-----
|
||||
|
||||
+10
@@ -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);
|
||||
}
|
||||
Vendored
+36
@@ -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 <cmath>
|
||||
#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);
|
||||
}
|
||||
}
|
||||
Vendored
+222
@@ -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 );
|
||||
}
|
||||
Vendored
+32
@@ -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;
|
||||
Vendored
+125
@@ -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);
|
||||
}
|
||||
Vendored
+479
@@ -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
|
||||
}
|
||||
}
|
||||
Vendored
+206
@@ -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 );
|
||||
}
|
||||
|
||||
Vendored
+176
@@ -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;
|
||||
Vendored
+102
@@ -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 <stdio.h>
|
||||
#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 );
|
||||
}
|
||||
Vendored
+97
@@ -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 );
|
||||
}
|
||||
Vendored
+270
@@ -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);
|
||||
}
|
||||
Vendored
+214
@@ -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 );
|
||||
}
|
||||
Reference in new issue
Block a user