Mercurial > hg > octave-nkf > gnulib-hg
changeset 4239:b33b8603715c
Transcendental functions for 'long double', from Paolo Bonzini.
author | Bruno Haible <bruno@clisp.org> |
---|---|
date | Tue, 18 Feb 2003 17:05:23 +0000 |
parents | b0b208dbaf73 |
children | 6918c3fd6d40 |
files | ChangeLog lib/ChangeLog lib/acosl.c lib/asinl.c lib/atanl.c lib/ceill.c lib/cosl.c lib/expl.c lib/floorl.c lib/frexpl.c lib/ldexpl.c lib/logl.c lib/mathl.h lib/sincosl.c lib/sinl.c lib/sqrtl.c lib/tanl.c lib/trigl.c lib/trigl.h m4/ChangeLog m4/mathl.m4 modules/mathl |
diffstat | 22 files changed, 3430 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,7 @@ +2003-02-18 Paolo Bonzino <bonzini@gnu.org> + + * modules/mathl: New file. + 2003-02-17 Bruno Haible <bruno@clisp.org> * modules/mkdtemp: New module.
--- a/lib/ChangeLog +++ b/lib/ChangeLog @@ -1,3 +1,23 @@ +2003-02-18 Paolo Bonzino <bonzini@gnu.org> + + * mathl.h: New file. + * acosl.c: New file. + * asinl.c: New file. + * atanl.c: New file. + * ceill.c: New file. + * cosl.c: New file. + * expl.c: New file. + * floorl.c: New file. + * frexpl.c: New file. + * ldexpl.c: New file. + * logl.c: New file. + * sincosl.c: New file. + * sinl.c: New file. + * sqrtl.c: New file. + * tanl.c: New file. + * trigl.c: New file. + * trigl.h: New file. + 2003-02-17 Bruno Haible <bruno@clisp.org> * mkdtemp.h: New file, from GNU gettext.
new file mode 100644 --- /dev/null +++ b/lib/acosl.c @@ -0,0 +1,224 @@ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "mathl.h" + +/* + Long double expansions contributed by + Stephen L. Moshier <moshier@na-net.ornl.gov> +*/ + +/* asin(x) + * Method : + * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... + * we approximate asin(x) on [0,0.5] by + * asin(x) = x + x*x^2*R(x^2) + * Between .5 and .625 the approximation is + * asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x) + * For x in [0.625,1] + * asin(x) = pi/2-2*asin(sqrt((1-x)/2)) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + * + */ + + +#include <math.h> + +static const long double + one = 1.0L, + huge = 1.0e+4932L, + pi = 3.1415926535897932384626433832795028841972L, + pio2_hi = 1.5707963267948966192313216916397514420986L, + pio2_lo = 4.3359050650618905123985220130216759843812E-35L, + pio4_hi = 7.8539816339744830961566084581987569936977E-1L, + + /* coefficient for R(x^2) */ + + /* asin(x) = x + x^3 pS(x^2) / qS(x^2) + 0 <= x <= 0.5 + peak relative error 1.9e-35 */ + pS0 = -8.358099012470680544198472400254596543711E2L, + pS1 = 3.674973957689619490312782828051860366493E3L, + pS2 = -6.730729094812979665807581609853656623219E3L, + pS3 = 6.643843795209060298375552684423454077633E3L, + pS4 = -3.817341990928606692235481812252049415993E3L, + pS5 = 1.284635388402653715636722822195716476156E3L, + pS6 = -2.410736125231549204856567737329112037867E2L, + pS7 = 2.219191969382402856557594215833622156220E1L, + pS8 = -7.249056260830627156600112195061001036533E-1L, + pS9 = 1.055923570937755300061509030361395604448E-3L, + + qS0 = -5.014859407482408326519083440151745519205E3L, + qS1 = 2.430653047950480068881028451580393430537E4L, + qS2 = -4.997904737193653607449250593976069726962E4L, + qS3 = 5.675712336110456923807959930107347511086E4L, + qS4 = -3.881523118339661268482937768522572588022E4L, + qS5 = 1.634202194895541569749717032234510811216E4L, + qS6 = -4.151452662440709301601820849901296953752E3L, + qS7 = 5.956050864057192019085175976175695342168E2L, + qS8 = -4.175375777334867025769346564600396877176E1L, + /* 1.000000000000000000000000000000000000000E0 */ + + /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x) + -0.0625 <= x <= 0.0625 + peak relative error 3.3e-35 */ + rS0 = -5.619049346208901520945464704848780243887E0L, + rS1 = 4.460504162777731472539175700169871920352E1L, + rS2 = -1.317669505315409261479577040530751477488E2L, + rS3 = 1.626532582423661989632442410808596009227E2L, + rS4 = -3.144806644195158614904369445440583873264E1L, + rS5 = -9.806674443470740708765165604769099559553E1L, + rS6 = 5.708468492052010816555762842394927806920E1L, + rS7 = 1.396540499232262112248553357962639431922E1L, + rS8 = -1.126243289311910363001762058295832610344E1L, + rS9 = -4.956179821329901954211277873774472383512E-1L, + rS10 = 3.313227657082367169241333738391762525780E-1L, + + sS0 = -4.645814742084009935700221277307007679325E0L, + sS1 = 3.879074822457694323970438316317961918430E1L, + sS2 = -1.221986588013474694623973554726201001066E2L, + sS3 = 1.658821150347718105012079876756201905822E2L, + sS4 = -4.804379630977558197953176474426239748977E1L, + sS5 = -1.004296417397316948114344573811562952793E2L, + sS6 = 7.530281592861320234941101403870010111138E1L, + sS7 = 1.270735595411673647119592092304357226607E1L, + sS8 = -1.815144839646376500705105967064792930282E1L, + sS9 = -7.821597334910963922204235247786840828217E-2L, + /* 1.000000000000000000000000000000000000000E0 */ + + asinr5625 = 5.9740641664535021430381036628424864397707E-1L; + + +long double +acosl (long double x) +{ + long double t, p, q; + + if (x < 0.0L) + { + t = pi - acosl(-x); + if (huge + x > one) /* return with inexact */ + return t; + } + + if (x >= 1.0L) /* |x|>= 1 */ + { + if (x == 1.0L) + return 0.0L; /* return zero */ + + return (x - x) / (x - x); /* asin(|x|>1) is NaN */ + } + + else if (x < 0.5L) /* |x| < 0.5 */ + { + if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* |x| < 2**-57 */ + /* acos(0)=+-pi/2 with inexact */ + return x * pio2_hi + x * pio2_lo; + + t = x * x; + p = (((((((((pS9 * t + + pS8) * t + + pS7) * t + + pS6) * t + + pS5) * t + + pS4) * t + + pS3) * t + + pS2) * t + + pS1) * t + + pS0) * t; + + q = (((((((( t + + qS8) * t + + qS7) * t + + qS6) * t + + qS5) * t + + qS4) * t + + qS3) * t + + qS2) * t + + qS1) * t + + qS0; + + return pio2_hi - (x + x * (p / q) - pio2_lo); + } + + else if (x < 0.625) /* 0.625 */ + { + t = x - 0.5625; + p = ((((((((((rS10 * t + + rS9) * t + + rS8) * t + + rS7) * t + + rS6) * t + + rS5) * t + + rS4) * t + + rS3) * t + + rS2) * t + + rS1) * t + + rS0) * t; + + q = ((((((((( t + + sS9) * t + + sS8) * t + + sS7) * t + + sS6) * t + + sS5) * t + + sS4) * t + + sS3) * t + + sS2) * t + + sS1) * t + + sS0; + + return (pio2_hi - asinr5625) - (p / q - pio2_lo); + } + else + return 2 * asinl(sqrtl((1-x)/2)); +} + +#if 0 +main() +{ + printf ("%.18Lg %.18Lg\n", + acosl(1.0L), + 1.5707963267948966192313216916397514420984L - + 1.5707963267948966192313216916397514420984L); + printf ("%.18Lg %.18Lg\n", + acosl(0.7071067811865475244008443621048490392848L), + 1.5707963267948966192313216916397514420984L - + 0.7853981633974483096156608458198757210492L); + printf ("%.18Lg %.18Lg\n", + acosl(0.5L), + 1.5707963267948966192313216916397514420984L - + 0.5235987755982988730771072305465838140328L); + printf ("%.18Lg %.18Lg\n", + acosl(0.3090169943749474241022934171828190588600L), + 1.5707963267948966192313216916397514420984L - + 0.3141592653589793238462643383279502884196L); + printf ("%.18Lg %.18Lg\n", + acosl(-1.0L), + 1.5707963267948966192313216916397514420984L - + -1.5707963267948966192313216916397514420984L); + printf ("%.18Lg %.18Lg\n", + acosl(-0.7071067811865475244008443621048490392848L), + 1.5707963267948966192313216916397514420984L - + -0.7853981633974483096156608458198757210492L); + printf ("%.18Lg %.18Lg\n", + acosl(-0.5L), + 1.5707963267948966192313216916397514420984L - + -0.5235987755982988730771072305465838140328L); + printf ("%.18Lg %.18Lg\n", + acosl(-0.3090169943749474241022934171828190588600L), + 1.5707963267948966192313216916397514420984L - + -0.3141592653589793238462643383279502884196L); +} +#endif
new file mode 100644 --- /dev/null +++ b/lib/asinl.c @@ -0,0 +1,218 @@ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "mathl.h" + +/* + Long double expansions contributed by + Stephen L. Moshier <moshier@na-net.ornl.gov> +*/ + +/* asin(x) + * Method : + * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... + * we approximate asin(x) on [0,0.5] by + * asin(x) = x + x*x^2*R(x^2) + * Between .5 and .625 the approximation is + * asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x) + * For x in [0.625,1] + * asin(x) = pi/2-2*asin(sqrt((1-x)/2)) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + * + */ + + +#include <math.h> + +static const long double + one = 1.0L, + huge = 1.0e+4932L, + pio2_hi = 1.5707963267948966192313216916397514420986L, + pio2_lo = 4.3359050650618905123985220130216759843812E-35L, + pio4_hi = 7.8539816339744830961566084581987569936977E-1L, + + /* coefficient for R(x^2) */ + + /* asin(x) = x + x^3 pS(x^2) / qS(x^2) + 0 <= x <= 0.5 + peak relative error 1.9e-35 */ + pS0 = -8.358099012470680544198472400254596543711E2L, + pS1 = 3.674973957689619490312782828051860366493E3L, + pS2 = -6.730729094812979665807581609853656623219E3L, + pS3 = 6.643843795209060298375552684423454077633E3L, + pS4 = -3.817341990928606692235481812252049415993E3L, + pS5 = 1.284635388402653715636722822195716476156E3L, + pS6 = -2.410736125231549204856567737329112037867E2L, + pS7 = 2.219191969382402856557594215833622156220E1L, + pS8 = -7.249056260830627156600112195061001036533E-1L, + pS9 = 1.055923570937755300061509030361395604448E-3L, + + qS0 = -5.014859407482408326519083440151745519205E3L, + qS1 = 2.430653047950480068881028451580393430537E4L, + qS2 = -4.997904737193653607449250593976069726962E4L, + qS3 = 5.675712336110456923807959930107347511086E4L, + qS4 = -3.881523118339661268482937768522572588022E4L, + qS5 = 1.634202194895541569749717032234510811216E4L, + qS6 = -4.151452662440709301601820849901296953752E3L, + qS7 = 5.956050864057192019085175976175695342168E2L, + qS8 = -4.175375777334867025769346564600396877176E1L, + /* 1.000000000000000000000000000000000000000E0 */ + + /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x) + -0.0625 <= x <= 0.0625 + peak relative error 3.3e-35 */ + rS0 = -5.619049346208901520945464704848780243887E0L, + rS1 = 4.460504162777731472539175700169871920352E1L, + rS2 = -1.317669505315409261479577040530751477488E2L, + rS3 = 1.626532582423661989632442410808596009227E2L, + rS4 = -3.144806644195158614904369445440583873264E1L, + rS5 = -9.806674443470740708765165604769099559553E1L, + rS6 = 5.708468492052010816555762842394927806920E1L, + rS7 = 1.396540499232262112248553357962639431922E1L, + rS8 = -1.126243289311910363001762058295832610344E1L, + rS9 = -4.956179821329901954211277873774472383512E-1L, + rS10 = 3.313227657082367169241333738391762525780E-1L, + + sS0 = -4.645814742084009935700221277307007679325E0L, + sS1 = 3.879074822457694323970438316317961918430E1L, + sS2 = -1.221986588013474694623973554726201001066E2L, + sS3 = 1.658821150347718105012079876756201905822E2L, + sS4 = -4.804379630977558197953176474426239748977E1L, + sS5 = -1.004296417397316948114344573811562952793E2L, + sS6 = 7.530281592861320234941101403870010111138E1L, + sS7 = 1.270735595411673647119592092304357226607E1L, + sS8 = -1.815144839646376500705105967064792930282E1L, + sS9 = -7.821597334910963922204235247786840828217E-2L, + /* 1.000000000000000000000000000000000000000E0 */ + + asinr5625 = 5.9740641664535021430381036628424864397707E-1L; + + +long double +asinl (long double x) +{ + long double y, t, p, q; + int sign; + + sign = 1; + y = x; + if (x < 0.0L) + { + sign = -1; + y = -x; + } + + if (y >= 1.0L) /* |x|>= 1 */ + { + if (y == 1.0L) + /* asin(1)=+-pi/2 with inexact */ + return x * pio2_hi + x * pio2_lo; + + return (x - x) / (x - x); /* asin(|x|>1) is NaN */ + } + else if (y < 0.5L) /* |x| < 0.5 */ + { + if (y < 0.000000000000000006938893903907228377647697925567626953125L) /* |x| < 2**-57 */ + if (huge + y > one) + return y; /* return x with inexact if x!=0 */ + + t = x * x; + p = (((((((((pS9 * t + + pS8) * t + + pS7) * t + + pS6) * t + + pS5) * t + + pS4) * t + + pS3) * t + + pS2) * t + + pS1) * t + + pS0) * t; + + q = (((((((( t + + qS8) * t + + qS7) * t + + qS6) * t + + qS5) * t + + qS4) * t + + qS3) * t + + qS2) * t + + qS1) * t + + qS0; + + return x + x * (p / q); + } + + else if (y < 0.625) /* 0.625 */ + { + t = y - 0.5625; + p = ((((((((((rS10 * t + + rS9) * t + + rS8) * t + + rS7) * t + + rS6) * t + + rS5) * t + + rS4) * t + + rS3) * t + + rS2) * t + + rS1) * t + + rS0) * t; + + q = ((((((((( t + + sS9) * t + + sS8) * t + + sS7) * t + + sS6) * t + + sS5) * t + + sS4) * t + + sS3) * t + + sS2) * t + + sS1) * t + + sS0; + t = asinr5625 + p / q; + } + else + t = pio2_hi + pio2_lo - 2 * asinl(sqrtl((1-y)/2)); + + return t * sign; +} + +#if 0 +main() +{ + printf ("%.18Lg %.18Lg\n", + asinl(1.0L), + 1.5707963267948966192313216916397514420984L); + printf ("%.18Lg %.18Lg\n", + asinl(0.7071067811865475244008443621048490392848L), + 0.7853981633974483096156608458198757210492L); + printf ("%.18Lg %.18Lg\n", + asinl(0.5L), + 0.5235987755982988730771072305465838140328L); + printf ("%.18Lg %.18Lg\n", + asinl(0.3090169943749474241022934171828190588600L), + 0.3141592653589793238462643383279502884196L); + printf ("%.18Lg %.18Lg\n", + asinl(-1.0L), + -1.5707963267948966192313216916397514420984L); + printf ("%.18Lg %.18Lg\n", + asinl(-0.7071067811865475244008443621048490392848L), + -0.7853981633974483096156608458198757210492L); + printf ("%.18Lg %.18Lg\n", + asinl(-0.5L), + -0.5235987755982988730771072305465838140328L); + printf ("%.18Lg %.18Lg\n", + asinl(-0.3090169943749474241022934171828190588600L), + -0.3141592653589793238462643383279502884196L); +} +#endif
new file mode 100644 --- /dev/null +++ b/lib/atanl.c @@ -0,0 +1,200 @@ +/* s_atanl.c + * + * Inverse circular tangent for 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. + * + * The function uses a rational approximation of the form + * t + t^3 P(t^2)/Q(t^2), optimized for |t| < 0.09375. + * + * The argument is reduced using the identity + * arctan x - arctan u = arctan ((x-u)/(1 + ux)) + * and an 83-entry lookup table for arctan u, with u = 0, 1/8, ..., 10.25. + * Use of the table improves the execution speed of the routine. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -19, 19 4e5 1.7e-34 5.4e-35 + * + * + * WARNING: + * + * This program uses integer operations on bit fields of floating-point + * numbers. It does not work with data structures other than the + * structure assumed. + * + */ + +/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov> */ + +#include "mathl.h" + +#include <math.h> + +/* arctan(k/8), k = 0, ..., 82 */ +static const long double atantbl[84] = { + 0.0000000000000000000000000000000000000000E0L, + 1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125) */ + 2.4497866312686415417208248121127581091414E-1L, + 3.5877067027057222039592006392646049977698E-1L, + 4.6364760900080611621425623146121440202854E-1L, + 5.5859931534356243597150821640166127034645E-1L, + 6.4350110879328438680280922871732263804151E-1L, + 7.1882999962162450541701415152590465395142E-1L, + 7.8539816339744830961566084581987572104929E-1L, + 8.4415398611317100251784414827164750652594E-1L, + 8.9605538457134395617480071802993782702458E-1L, + 9.4200004037946366473793717053459358607166E-1L, + 9.8279372324732906798571061101466601449688E-1L, + 1.0191413442663497346383429170230636487744E0L, + 1.0516502125483736674598673120862998296302E0L, + 1.0808390005411683108871567292171998202703E0L, + 1.1071487177940905030170654601785370400700E0L, + 1.1309537439791604464709335155363278047493E0L, + 1.1525719972156675180401498626127513797495E0L, + 1.1722738811284763866005949441337046149712E0L, + 1.1902899496825317329277337748293183376012E0L, + 1.2068173702852525303955115800565576303133E0L, + 1.2220253232109896370417417439225704908830E0L, + 1.2360594894780819419094519711090786987027E0L, + 1.2490457723982544258299170772810901230778E0L, + 1.2610933822524404193139408812473357720101E0L, + 1.2722973952087173412961937498224804940684E0L, + 1.2827408797442707473628852511364955306249E0L, + 1.2924966677897852679030914214070816845853E0L, + 1.3016288340091961438047858503666855921414E0L, + 1.3101939350475556342564376891719053122733E0L, + 1.3182420510168370498593302023271362531155E0L, + 1.3258176636680324650592392104284756311844E0L, + 1.3329603993374458675538498697331558093700E0L, + 1.3397056595989995393283037525895557411039E0L, + 1.3460851583802539310489409282517796256512E0L, + 1.3521273809209546571891479413898128509842E0L, + 1.3578579772154994751124898859640585287459E0L, + 1.3633001003596939542892985278250991189943E0L, + 1.3684746984165928776366381936948529556191E0L, + 1.3734007669450158608612719264449611486510E0L, + 1.3780955681325110444536609641291551522494E0L, + 1.3825748214901258580599674177685685125566E0L, + 1.3868528702577214543289381097042486034883E0L, + 1.3909428270024183486427686943836432060856E0L, + 1.3948567013423687823948122092044222644895E0L, + 1.3986055122719575950126700816114282335732E0L, + 1.4021993871854670105330304794336492676944E0L, + 1.4056476493802697809521934019958079881002E0L, + 1.4089588955564736949699075250792569287156E0L, + 1.4121410646084952153676136718584891599630E0L, + 1.4152014988178669079462550975833894394929E0L, + 1.4181469983996314594038603039700989523716E0L, + 1.4209838702219992566633046424614466661176E0L, + 1.4237179714064941189018190466107297503086E0L, + 1.4263547484202526397918060597281265695725E0L, + 1.4288992721907326964184700745371983590908E0L, + 1.4313562697035588982240194668401779312122E0L, + 1.4337301524847089866404719096698873648610E0L, + 1.4360250423171655234964275337155008780675E0L, + 1.4382447944982225979614042479354815855386E0L, + 1.4403930189057632173997301031392126865694E0L, + 1.4424730991091018200252920599377292525125E0L, + 1.4444882097316563655148453598508037025938E0L, + 1.4464413322481351841999668424758804165254E0L, + 1.4483352693775551917970437843145232637695E0L, + 1.4501726582147939000905940595923466567576E0L, + 1.4519559822271314199339700039142990228105E0L, + 1.4536875822280323362423034480994649820285E0L, + 1.4553696664279718992423082296859928222270E0L, + 1.4570043196511885530074841089245667532358E0L, + 1.4585935117976422128825857356750737658039E0L, + 1.4601391056210009726721818194296893361233E0L, + 1.4616428638860188872060496086383008594310E0L, + 1.4631064559620759326975975316301202111560E0L, + 1.4645314639038178118428450961503371619177E0L, + 1.4659193880646627234129855241049975398470E0L, + 1.4672716522843522691530527207287398276197E0L, + 1.4685896086876430842559640450619880951144E0L, + 1.4698745421276027686510391411132998919794E0L, + 1.4711276743037345918528755717617308518553E0L, + 1.4723501675822635384916444186631899205983E0L, + 1.4735431285433308455179928682541563973416E0L, /* arctan(10.25) */ + 1.5707963267948966192313216916397514420986E0L /* pi/2 */ +}; + + +/* arctan t = t + t^3 p(t^2) / q(t^2) + |t| <= 0.09375 + peak relative error 5.3e-37 */ + +static const long double + p0 = -4.283708356338736809269381409828726405572E1L, + p1 = -8.636132499244548540964557273544599863825E1L, + p2 = -5.713554848244551350855604111031839613216E1L, + p3 = -1.371405711877433266573835355036413750118E1L, + p4 = -8.638214309119210906997318946650189640184E-1L, + q0 = 1.285112506901621042780814422948906537959E2L, + q1 = 3.361907253914337187957855834229672347089E2L, + q2 = 3.180448303864130128268191635189365331680E2L, + q3 = 1.307244136980865800160844625025280344686E2L, + q4 = 2.173623741810414221251136181221172551416E1L; + /* q5 = 1.000000000000000000000000000000000000000E0 */ + + +long double +atanl (long double x) +{ + int k, sign; + long double t, u, p, q; + + sign = (x < 0) ? -1 : 1; + + /* Check for zero or NaN. */ + if (x != x || x == 0.0) + return x + x; + + /* Check for infinity. */ + if (x + x == x) + return sign * atantbl[83]; + + x *= sign; + + if (x >= 10.25) /* 10.25 */ + { + k = 83; + t = -1.0/x; + } + else + { + /* Index of nearest table element. + Roundoff to integer is asymmetrical to avoid cancellation when t < 0 + (cf. fdlibm). */ + k = 8.0 * x + 0.25; + u = 0.125 * k; + /* Small arctan argument. */ + t = (x - u) / (1.0 + x * u); + } + + /* Arctan of small argument t. */ + u = t * t; + p = ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0; + q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0; + u = t * u * p / q + t; + + /* arctan x = arctan u + arctan t */ + return sign * (atantbl[k] + u); +}
new file mode 100644 --- /dev/null +++ b/lib/ceill.c @@ -0,0 +1,45 @@ +/* Emulation for ceill. + Contributed by Paolo Bonzini + + Copyright 2002, 2003 Free Software Foundation, Inc. + + This file is part of gnulib. + + gnulib is free software; you can redistribute it and/or modify it + under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation; either version 2.1, or (at your option) + any later version. + + gnulib is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public + License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with gnulib; see the file COPYING.LIB. If not, write to the Free + Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, + USA. + */ + +#include <float.h> + +#include "mathl.h" + +/* To compute the integer part of X, sum a big enough + integer so that the precision of the floating point + number is exactly 1. */ + +long double +ceill(long double x) +{ + long double y; + if (x < 0.0L) + y = -(1.0L / LDBL_EPSILON - x - 1.0 / LDBL_EPSILON); + else + y = 1.0L / LDBL_EPSILON + x - 1.0 / LDBL_EPSILON; + + if (y < x) + return y + 1.0L; + else + return y; +}
new file mode 100644 --- /dev/null +++ b/lib/cosl.c @@ -0,0 +1,97 @@ +/* s_cosl.c -- long double version of s_sin.c. + * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* sinl(x) + * Return sine function of x. + * + * kernel function: + * __kernel_sinl ... sine function on [-pi/4,pi/4] + * __kernel_cosl ... cose function on [-pi/4,pi/4] + * __ieee754_rem_pio2l ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ + +#include <math.h> + +#include "mathl.h" + +#include "trigl.h" +#ifdef HAVE_SINL +#include "trigl.c" +#include "sincosl.c" +#endif + +long double cosl(long double x) +{ + long double y[2],z=0.0L; + int n; + + /* |x| ~< pi/4 */ + if(x >= -0.7853981633974483096156608458198757210492 && + x <= 0.7853981633974483096156608458198757210492) + return kernel_cosl(x, z); + + /* sinl(Inf or NaN) is NaN, sinl(0) is 0 */ + else if ((x + x == x && x != 0.0) || x != x) + return x-x; /* NaN */ + + /* argument reduction needed */ + else { + n = ieee754_rem_pio2l(x,y); + switch(n&3) { + case 0: return kernel_cosl(y[0],y[1]); + case 1: return -kernel_sinl(y[0],y[1],1); + case 2: return -kernel_cosl(y[0],y[1]); + default: + return kernel_sinl(y[0],y[1],1); + } + } +} + +#if 0 +int +main () +{ + printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492)); + printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *29)); + printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *2)); + printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *30)); + printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *4)); + printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *32)); + printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *2/3)); + printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *4/3)); +} +#endif
new file mode 100644 --- /dev/null +++ b/lib/expl.c @@ -0,0 +1,136 @@ +/* Emulation for expl. + Contributed by Paolo Bonzini + + Copyright 2002, 2003 Free Software Foundation, Inc. + + This file is part of gnulib. + + gnulib is free software; you can redistribute it and/or modify it + under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation; either version 2.1, or (at your option) + any later version. + + gnulib is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public + License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with gnulib; see the file COPYING.LIB. If not, write to the Free + Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, + USA. + */ + +#include <float.h> +#include <math.h> + +#include "mathl.h" + +static const long double C[] = { +/* Chebyshev polynom coeficients for (exp(x)-1)/x */ +#define P1 C[0] +#define P2 C[1] +#define P3 C[2] +#define P4 C[3] +#define P5 C[4] +#define P6 C[5] + 0.5L, + 1.66666666666666666666666666666666683E-01L, + 4.16666666666666666666654902320001674E-02L, + 8.33333333333333333333314659767198461E-03L, + 1.38888888889899438565058018857254025E-03L, + 1.98412698413981650382436541785404286E-04L, + +/* Smallest integer x for which e^x overflows. */ +#define himark C[6] + 11356.523406294143949491931077970765L, + +/* Largest integer x for which e^x underflows. */ +#define lomark C[7] +-11433.4627433362978788372438434526231L, + +/* very small number */ +#define TINY C[8] + 1.0e-4900L, + +/* 2^16383 */ +#define TWO16383 C[9] + 5.94865747678615882542879663314003565E+4931L}; + +long double +expl (long double x) +{ + /* Check for usual case. */ + if (x < himark && x > lomark) + { + int exponent; + long double t, x22; + int k = 1; + long double result = 1.0; + + /* Compute an integer power of e with a granularity of 0.125. */ + exponent = (int) floorl (x * 8.0L); + x -= exponent / 8.0L; + + if (x > 0.0625) + { + exponent++; + x -= 0.125L; + } + + if (exponent < 0) + { + t = 0.8824969025845954028648921432290507362220L; /* e^-0.25 */ + exponent = -exponent; + } + else + t = 1.1331484530668263168290072278117938725655L; /* e^0.25 */ + + while (exponent) + { + if (exponent & k) + { + result *= t; + exponent ^= k; + } + t *= t; + k <<= 1; + } + + /* Approximate (e^x - 1)/x, using a seventh-degree polynomial, + with maximum error in [-2^-16-2^-53,2^-16+2^-53] + less than 4.8e-39. */ + x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6))))); + + return result + result * x22; + } + /* Exceptional cases: */ + else if (x < himark) + { + if (x + x == x) + /* e^-inf == 0, with no error. */ + return 0; + else + /* Underflow */ + return TINY * TINY; + } + else + /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ + return TWO16383*x; +} + +#if 0 +int +main () +{ + printf ("%.16Lg\n", expl(1.0L)); + printf ("%.16Lg\n", expl(-1.0L)); + printf ("%.16Lg\n", expl(2.0L)); + printf ("%.16Lg\n", expl(4.0L)); + printf ("%.16Lg\n", expl(-2.0L)); + printf ("%.16Lg\n", expl(-4.0L)); + printf ("%.16Lg\n", expl(0.0625L)); + printf ("%.16Lg\n", expl(0.3L)); + printf ("%.16Lg\n", expl(0.6L)); +} +#endif
new file mode 100644 --- /dev/null +++ b/lib/floorl.c @@ -0,0 +1,45 @@ +/* Emulation for floorl. + Contributed by Paolo Bonzini + + Copyright 2002, 2003 Free Software Foundation, Inc. + + This file is part of gnulib. + + gnulib is free software; you can redistribute it and/or modify it + under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation; either version 2.1, or (at your option) + any later version. + + gnulib is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public + License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with gnulib; see the file COPYING.LIB. If not, write to the Free + Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, + USA. + */ + +#include <float.h> + +#include "mathl.h" + +/* To compute the integer part of X, sum a big enough + integer so that the precision of the floating point + number is exactly 1. */ + +long double +floorl(long double x) +{ + long double y; + if (x < 0.0L) + y = -(1.0L / LDBL_EPSILON - x - 1.0 / LDBL_EPSILON); + else + y = 1.0L / LDBL_EPSILON + x - 1.0 / LDBL_EPSILON; + + if (y > x) + return y - 1.0L; + else + return y; +}
new file mode 100644 --- /dev/null +++ b/lib/frexpl.c @@ -0,0 +1,103 @@ +/* Emulation for frexpl. + Contributed by Paolo Bonzini + + Copyright 2002, 2003 Free Software Foundation, Inc. + + This file is part of gnulib. + + gnulib is free software; you can redistribute it and/or modify it + under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation; either version 2.1, or (at your option) + any later version. + + gnulib is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public + License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with gnulib; see the file COPYING.LIB. If not, write to the Free + Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, + USA. + */ + +#include <float.h> +#include <math.h> + +#include "mathl.h" + +/* Binary search. Quite inefficient but portable. */ +long double +frexpl(long double x, int *exp) +{ + long double exponents[20], *next; + int exponent, bit; + + /* Check for zero, nan and infinity. */ + if (x != x || x + x == x ) + { + *exp = 0; + return x; + } + + if (x < 0) + return -frexpl(-x, exp); + + exponent = 0; + if (x > 1.0) + { + for (next = exponents, exponents[0] = 2.0L, bit = 1; + *next <= x + x; + bit <<= 1, next[1] = next[0] * next[0], next++); + + for (; next >= exponents; bit >>= 1, next--) + if (x + x >= *next) + { + x /= *next; + exponent |= bit; + } + + } + + else if (x < 0.5) + { + for (next = exponents, exponents[0] = 0.5L, bit = 1; + *next > x; + bit <<= 1, next[1] = next[0] * next[0], next++); + + for (; next >= exponents; bit >>= 1, next--) + if (x < *next) + { + x /= *next; + exponent |= bit; + } + + exponent = -exponent; + } + + *exp = exponent; + return x; +} + +#if 0 +int +main() +{ + long double x; + int y; + x = frexpl(0.0L / 0.0L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(1.0L / 0.0L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(-1.0L / 0.0L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(0.5L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(0.75L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(3.6L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(17.8L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(8.0L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(0.3L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(0.49L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(0.049L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(0.0245L, &y); printf ("%.6Lg %d\n", x, y); + x = frexpl(0.0625L, &y); printf ("%.6Lg %d\n", x, y); +} +#endif +
new file mode 100644 --- /dev/null +++ b/lib/ldexpl.c @@ -0,0 +1,63 @@ +/* Emulation for ldexpl. + Contributed by Paolo Bonzini + + Copyright 2002, 2003 Free Software Foundation, Inc. + + This file is part of gnulib. + + gnulib is free software; you can redistribute it and/or modify it + under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation; either version 2.1, or (at your option) + any later version. + + gnulib is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public + License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with gnulib; see the file COPYING.LIB. If not, write to the Free + Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, + USA. + */ + +#include <float.h> +#include <math.h> + +#include "mathl.h" + +long double +ldexpl(long double x, int exp) +{ + long double factor; + int bit; + + /* Check for zero, nan and infinity. */ + if (x != x || x + x == x ) + return x; + + if (exp < 0) + { + exp = -exp; + factor = 0.5L; + } + else + factor = 2.0L; + + for (bit = 1; bit <= exp; bit <<= 1, factor *= factor) + if (exp & bit) + x *= factor; + + return x; +} + +#if 0 +int +main() +{ + long double x; + int y; + for (y = 0; y < 29; y++) + printf ("%5d %.16Lg %.16Lg\n", y, ldexpl(0.8L, y), ldexpl(0.8L, -y) * ldexpl(0.8L, y)); +} +#endif
new file mode 100644 --- /dev/null +++ b/lib/logl.c @@ -0,0 +1,238 @@ +/* logll.c + * + * Natural logarithm for 128-bit long double precision. + * + * + * + * SYNOPSIS: + * + * long double x, y, logl(); + * + * y = logl( x ); + * + * + * + * DESCRIPTION: + * + * Returns the base e (2.718...) logarithm of x. + * + * The argument is separated into its exponent and fractional + * parts. Use of a lookup table increases the speed of the routine. + * The program uses logarithms tabulated at intervals of 1/128 to + * cover the domain from approximately 0.7 to 1.4. + * + * On the interval [-1/128, +1/128] the logarithm of 1+x is approximated by + * log(1+x) = x - 0.5 x^2 + x^3 P(x) . + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0.875, 1.125 100000 1.2e-34 4.1e-35 + * IEEE 0.125, 8 100000 1.2e-34 4.1e-35 + * + * + * WARNING: + * + * This program uses integer operations on bit fields of floating-point + * numbers. It does not work with data structures other than the + * structure assumed. + * + */ + +#include <math.h> + +#include "mathl.h" + +/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov> */ + +/* log(1+x) = x - .5 x^2 + x^3 l(x) + -.0078125 <= x <= +.0078125 + peak relative error 1.2e-37 */ +static const long double +l3 = 3.333333333333333333333333333333336096926E-1L, +l4 = -2.499999999999999999999999999486853077002E-1L, +l5 = 1.999999999999999999999999998515277861905E-1L, +l6 = -1.666666666666666666666798448356171665678E-1L, +l7 = 1.428571428571428571428808945895490721564E-1L, +l8 = -1.249999999999999987884655626377588149000E-1L, +l9 = 1.111111111111111093947834982832456459186E-1L, +l10 = -1.000000000000532974938900317952530453248E-1L, +l11 = 9.090909090915566247008015301349979892689E-2L, +l12 = -8.333333211818065121250921925397567745734E-2L, +l13 = 7.692307559897661630807048686258659316091E-2L, +l14 = -7.144242754190814657241902218399056829264E-2L, +l15 = 6.668057591071739754844678883223432347481E-2L; + +/* Lookup table of ln(t) - (t-1) + t = 0.5 + (k+26)/128) + k = 0, ..., 91 */ +static const long double logtbl[92] = { +-5.5345593589352099112142921677820359632418E-2L, +-5.2108257402767124761784665198737642086148E-2L, +-4.8991686870576856279407775480686721935120E-2L, +-4.5993270766361228596215288742353061431071E-2L, +-4.3110481649613269682442058976885699556950E-2L, +-4.0340872319076331310838085093194799765520E-2L, +-3.7682072451780927439219005993827431503510E-2L, +-3.5131785416234343803903228503274262719586E-2L, +-3.2687785249045246292687241862699949178831E-2L, +-3.0347913785027239068190798397055267411813E-2L, +-2.8110077931525797884641940838507561326298E-2L, +-2.5972247078357715036426583294246819637618E-2L, +-2.3932450635346084858612873953407168217307E-2L, +-2.1988775689981395152022535153795155900240E-2L, +-2.0139364778244501615441044267387667496733E-2L, +-1.8382413762093794819267536615342902718324E-2L, +-1.6716169807550022358923589720001638093023E-2L, +-1.5138929457710992616226033183958974965355E-2L, +-1.3649036795397472900424896523305726435029E-2L, +-1.2244881690473465543308397998034325468152E-2L, +-1.0924898127200937840689817557742469105693E-2L, +-9.6875626072830301572839422532631079809328E-3L, +-8.5313926245226231463436209313499745894157E-3L, +-7.4549452072765973384933565912143044991706E-3L, +-6.4568155251217050991200599386801665681310E-3L, +-5.5356355563671005131126851708522185605193E-3L, +-4.6900728132525199028885749289712348829878E-3L, +-3.9188291218610470766469347968659624282519E-3L, +-3.2206394539524058873423550293617843896540E-3L, +-2.5942708080877805657374888909297113032132E-3L, +-2.0385211375711716729239156839929281289086E-3L, +-1.5522183228760777967376942769773768850872E-3L, +-1.1342191863606077520036253234446621373191E-3L, +-7.8340854719967065861624024730268350459991E-4L, +-4.9869831458030115699628274852562992756174E-4L, +-2.7902661731604211834685052867305795169688E-4L, +-1.2335696813916860754951146082826952093496E-4L, +-3.0677461025892873184042490943581654591817E-5L, +#define ZERO logtbl[38] + 0.0000000000000000000000000000000000000000E0L, +-3.0359557945051052537099938863236321874198E-5L, +-1.2081346403474584914595395755316412213151E-4L, +-2.7044071846562177120083903771008342059094E-4L, +-4.7834133324631162897179240322783590830326E-4L, +-7.4363569786340080624467487620270965403695E-4L, +-1.0654639687057968333207323853366578860679E-3L, +-1.4429854811877171341298062134712230604279E-3L, +-1.8753781835651574193938679595797367137975E-3L, +-2.3618380914922506054347222273705859653658E-3L, +-2.9015787624124743013946600163375853631299E-3L, +-3.4938307889254087318399313316921940859043E-3L, +-4.1378413103128673800485306215154712148146E-3L, +-4.8328735414488877044289435125365629849599E-3L, +-5.5782063183564351739381962360253116934243E-3L, +-6.3731336597098858051938306767880719015261E-3L, +-7.2169643436165454612058905294782949315193E-3L, +-8.1090214990427641365934846191367315083867E-3L, +-9.0486422112807274112838713105168375482480E-3L, +-1.0035177140880864314674126398350812606841E-2L, +-1.1067990155502102718064936259435676477423E-2L, +-1.2146457974158024928196575103115488672416E-2L, +-1.3269969823361415906628825374158424754308E-2L, +-1.4437927104692837124388550722759686270765E-2L, +-1.5649743073340777659901053944852735064621E-2L, +-1.6904842527181702880599758489058031645317E-2L, +-1.8202661505988007336096407340750378994209E-2L, +-1.9542647000370545390701192438691126552961E-2L, +-2.0924256670080119637427928803038530924742E-2L, +-2.2346958571309108496179613803760727786257E-2L, +-2.3810230892650362330447187267648486279460E-2L, +-2.5313561699385640380910474255652501521033E-2L, +-2.6856448685790244233704909690165496625399E-2L, +-2.8438398935154170008519274953860128449036E-2L, +-3.0058928687233090922411781058956589863039E-2L, +-3.1717563112854831855692484086486099896614E-2L, +-3.3413836095418743219397234253475252001090E-2L, +-3.5147290019036555862676702093393332533702E-2L, +-3.6917475563073933027920505457688955423688E-2L, +-3.8723951502862058660874073462456610731178E-2L, +-4.0566284516358241168330505467000838017425E-2L, +-4.2444048996543693813649967076598766917965E-2L, +-4.4356826869355401653098777649745233339196E-2L, +-4.6304207416957323121106944474331029996141E-2L, +-4.8285787106164123613318093945035804818364E-2L, +-5.0301169421838218987124461766244507342648E-2L, +-5.2349964705088137924875459464622098310997E-2L, +-5.4431789996103111613753440311680967840214E-2L, +-5.6546268881465384189752786409400404404794E-2L, +-5.8693031345788023909329239565012647817664E-2L, +-6.0871713627532018185577188079210189048340E-2L, +-6.3081958078862169742820420185833800925568E-2L, +-6.5323413029406789694910800219643791556918E-2L, +-6.7595732653791419081537811574227049288168E-2L +}; + +/* ln(2) = ln2a + ln2b with extended precision. */ +static const long double + ln2a = 6.93145751953125e-1L, + ln2b = 1.4286068203094172321214581765680755001344E-6L; + +long double +logl(long double x) +{ + long double z, y, w, u, t; + unsigned int m; + int k, e; + + /* Check for IEEE special cases. */ + + /* log(0) = -infinity. */ + if (x == 0.0L) + return -0.5L / 0.0L; + + /* log ( x < 0 ) = NaN */ + if (x < 0.0L) + return (x - x) / (x - x); + + /* log (infinity or NaN) */ + if (x + x == x || x != x) + return x + x; + + /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625 */ + x = frexpl(x, &e); + if (x < 0.703125L) + { + x += x; + e--; + } + + /* On this interval the table is not used due to cancellation error. */ + if ((x <= 1.0078125L) && (x >= 0.9921875L)) + { + k = 64; + t = 1.0L; + z = x - 1.0L; + } + else + { + k = floorl((x - 0.5L) * 128.0L); + t = 0.5L + k / 128.0L; + z = (x - t) / t; + } + + /* Series expansion of log(1+z). */ + w = z * z; + y = ((((((((((((l15 * z + + l14) * z + + l13) * z + + l12) * z + + l11) * z + + l10) * z + + l9) * z + + l8) * z + + l7) * z + + l6) * z + + l5) * z + + l4) * z + + l3) * z * w; + y -= 0.5 * w; + y += e * ln2b; /* Base 2 exponent offset times ln(2). */ + y += z; + y += logtbl[k-26]; /* log(t) - (t-1) */ + y += (t - 1.0L); + y += e * ln2a; + return y; +}
new file mode 100644 --- /dev/null +++ b/lib/mathl.h @@ -0,0 +1,41 @@ +/* Declaration for long double functions. + Contributed by Paolo Bonzini + + Copyright 2002, 2003 Free Software Foundation, Inc. + + This file is part of gnulib. + + gnulib is free software; you can redistribute it and/or modify it + under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation; either version 2.1, or (at your option) + any later version. + + gnulib is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public + License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with gnulib; see the file COPYING.LIB. If not, write to the Free + Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, + USA. + */ + +#ifndef GNULIB_MATHL_H +#define GNULIB_MATHL_H + +extern long double acosl (long double x); +extern long double asinl (long double x); +extern long double atanl (long double x); +extern long double ceill (long double x); +extern long double cosl (long double x); +extern long double expl (long double x); +extern long double floorl (long double x); +extern long double frexpl (long double x, int *exp); +extern long double ldexpl (long double x, int exp); +extern long double logl (long double x); +extern long double sinl (long double x); +extern long double sqrtl (long double x); +extern long double tanl (long double x); + +#endif
new file mode 100644 --- /dev/null +++ b/lib/sincosl.c @@ -0,0 +1,897 @@ +/* Quad-precision floating point trigonometric functions on <-pi/4,pi/4>. + Copyright (C) 1999 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Jakub Jelinek <jj@ultra.linux.cz> + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include <math.h> +#include <float.h> + +#include "mathl.h" + +static const long double sin_c[] = { +#define ONE sin_c[0] + 1.00000000000000000000000000000000000E+00L, /* 3fff0000000000000000000000000000 */ + +/* cos x ~ ONE + x^2 ( SCOS1 + SCOS2 * x^2 + ... + SCOS4 * x^6 + SCOS5 * x^8 ) + x in <0,1/256> */ +#define SCOS1 sin_c[1] +#define SCOS2 sin_c[2] +#define SCOS3 sin_c[3] +#define SCOS4 sin_c[4] +#define SCOS5 sin_c[5] + -5.00000000000000000000000000000000000E-01L, /* bffe0000000000000000000000000000 */ + 4.16666666666666666666666666556146073E-02L, /* 3ffa5555555555555555555555395023 */ + -1.38888888888888888888309442601939728E-03L, /* bff56c16c16c16c16c16a566e42c0375 */ + 2.48015873015862382987049502531095061E-05L, /* 3fefa01a01a019ee02dcf7da2d6d5444 */ + -2.75573112601362126593516899592158083E-07L, /* bfe927e4f5dce637cb0b54908754bde0 */ + +/* sin x ~ ONE * x + x^3 ( SIN1 + SIN2 * x^2 + ... + SIN7 * x^12 + SIN8 * x^14 ) + x in <0,0.1484375> */ +#define SIN1 sin_c[6] +#define SIN2 sin_c[7] +#define SIN3 sin_c[8] +#define SIN4 sin_c[9] +#define SIN5 sin_c[10] +#define SIN6 sin_c[11] +#define SIN7 sin_c[12] +#define SIN8 sin_c[13] + -1.66666666666666666666666666666666538e-01L, /* bffc5555555555555555555555555550 */ + 8.33333333333333333333333333307532934e-03L, /* 3ff811111111111111111111110e7340 */ + -1.98412698412698412698412534478712057e-04L, /* bff2a01a01a01a01a01a019e7a626296 */ + 2.75573192239858906520896496653095890e-06L, /* 3fec71de3a556c7338fa38527474b8f5 */ + -2.50521083854417116999224301266655662e-08L, /* bfe5ae64567f544e16c7de65c2ea551f */ + 1.60590438367608957516841576404938118e-10L, /* 3fde6124613a811480538a9a41957115 */ + -7.64716343504264506714019494041582610e-13L, /* bfd6ae7f3d5aef30c7bc660b060ef365 */ + 2.81068754939739570236322404393398135e-15L, /* 3fce9510115aabf87aceb2022a9a9180 */ + +/* sin x ~ ONE * x + x^3 ( SSIN1 + SSIN2 * x^2 + ... + SSIN4 * x^6 + SSIN5 * x^8 ) + x in <0,1/256> */ +#define SSIN1 sin_c[14] +#define SSIN2 sin_c[15] +#define SSIN3 sin_c[16] +#define SSIN4 sin_c[17] +#define SSIN5 sin_c[18] + -1.66666666666666666666666666666666659E-01L, /* bffc5555555555555555555555555555 */ + 8.33333333333333333333333333146298442E-03L, /* 3ff81111111111111111111110fe195d */ + -1.98412698412698412697726277416810661E-04L, /* bff2a01a01a01a01a019e7121e080d88 */ + 2.75573192239848624174178393552189149E-06L, /* 3fec71de3a556c640c6aaa51aa02ab41 */ + -2.50521016467996193495359189395805639E-08L, /* bfe5ae644ee90c47dc71839de75b2787 */ +}; + +static const long double cos_c[] = { +/* cos x ~ ONE + x^2 ( COS1 + COS2 * x^2 + ... + COS7 * x^12 + COS8 * x^14 ) + x in <0,0.1484375> */ +#define COS1 cos_c[0] +#define COS2 cos_c[1] +#define COS3 cos_c[2] +#define COS4 cos_c[3] +#define COS5 cos_c[4] +#define COS6 cos_c[5] +#define COS7 cos_c[6] +#define COS8 cos_c[7] + -4.99999999999999999999999999999999759E-01L, /* bffdfffffffffffffffffffffffffffb */ + 4.16666666666666666666666666651287795E-02L, /* 3ffa5555555555555555555555516f30 */ + -1.38888888888888888888888742314300284E-03L, /* bff56c16c16c16c16c16c16a463dfd0d */ + 2.48015873015873015867694002851118210E-05L, /* 3fefa01a01a01a01a0195cebe6f3d3a5 */ + -2.75573192239858811636614709689300351E-07L, /* bfe927e4fb7789f5aa8142a22044b51f */ + 2.08767569877762248667431926878073669E-09L, /* 3fe21eed8eff881d1e9262d7adff4373 */ + -1.14707451049343817400420280514614892E-11L, /* bfda9397496922a9601ed3d4ca48944b */ + 4.77810092804389587579843296923533297E-14L, /* 3fd2ae5f8197cbcdcaf7c3fb4523414c */ + +}; + +#define SINCOSL_COS_HI 0 +#define SINCOSL_COS_LO 1 +#define SINCOSL_SIN_HI 2 +#define SINCOSL_SIN_LO 3 +static const long double sincosl_table[]; + +long double +kernel_sinl (long double x, long double y, int iy) +{ + long double h, l, z, sin_l, cos_l_m1; + int index, sign; + + sign = 1; + if (x < 0) + { + x = -x; + sign = -1; + } + + if (x < 0.148375L) /* |x| < 0.1484375 */ + { + /* Argument is small enough to approximate it by a Chebyshev + polynomial of degree 17. */ + if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* |x| < 2^-57 */ + if (!((int) x)) + return x; /* generate inexact */ + + z = x * x; + return x + (x * (z * (SIN1 + z * (SIN2 + z * (SIN3 + z * (SIN4 + + z * (SIN5 + + z * + (SIN6 + + z * + (SIN7 + + z * + SIN8))))))))); + } + else + { + /* So that we don't have to use too large polynomial, we find + l and h such that x = l + h, where fabsl(l) <= 1.0/256 with 83 + possible values for h. We look up cosl(h) and sinl(h) in + pre-computed tables, compute cosl(l) and sinl(l) using a + Chebyshev polynomial of degree 10(11) and compute + sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l). */ + x -= 0.1484375L; + index = (int) (x * 128L + 0.5L); + h = index / 128.0L; + if (iy) + l = y - (h - x); + else + l = x - h; + + z = l * l; + sin_l = + l * (ONE + + z * (SSIN1 + + z * (SSIN2 + z * (SSIN3 + z * (SSIN4 + z * SSIN5))))); + cos_l_m1 = + z * (SCOS1 + z * (SCOS2 + z * (SCOS3 + z * (SCOS4 + z * SCOS5)))); + + index *= 4; + z = + sincosl_table[index + SINCOSL_SIN_HI] + + (sincosl_table[index + SINCOSL_SIN_LO] + + (sincosl_table[index + SINCOSL_SIN_HI] * cos_l_m1) + + (sincosl_table[index + SINCOSL_COS_HI] * sin_l)); + return z * sign; + } +} + +long double +kernel_cosl (long double x, long double y) +{ + long double h, l, z, sin_l, cos_l_m1; + int index; + + if (x < 0) + x = -x; + + if (x < 0.1484375L) /* |x| < 0.1484375 */ + { + /* Argument is small enough to approximate it by a Chebyshev + polynomial of degree 16. */ + if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* |x| < 2^-57 */ + if (!((int) x)) + return ONE; /* generate inexact */ + z = x * x; + return ONE + (z * (COS1 + z * (COS2 + z * (COS3 + z * (COS4 + + z * (COS5 + + z * (COS6 + + z * + (COS7 + + z * + COS8)))))))); + } + else + { + /* So that we don't have to use too large polynomial, we find + l and h such that x = l + h, where fabsl(l) <= 1.0/256 with 83 + possible values for h. We look up cosl(h) and sinl(h) in + pre-computed tables, compute cosl(l) and sinl(l) using a + Chebyshev polynomial of degree 10(11) and compute + sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l). */ + x -= 0.1484375L; + index = (int) (x * 128L + 0.5L); + h = index / 128.0L; + l = y - (h - x); + z = l * l; + sin_l = + l * (ONE + + z * (SSIN1 + + z * (SSIN2 + z * (SSIN3 + z * (SSIN4 + z * SSIN5))))); + cos_l_m1 = + z * (SCOS1 + z * (SCOS2 + z * (SCOS3 + z * (SCOS4 + z * SCOS5)))); + + index *= 4; + z = sincosl_table [index + SINCOSL_COS_HI] + + (sincosl_table [index + SINCOSL_COS_LO] + - (sincosl_table [index + SINCOSL_SIN_HI] * sin_l) + - (sincosl_table [index + SINCOSL_COS_HI] * cos_l_m1)); + return z; + } +} + + +/* For 0.1484375 + n/128.0, n=0..82 this table contains + first 113 bits of cosine, then at least 113 additional + bits and the same for sine. + 0.1484375+82.0/128.0 is the smallest number among above defined numbers + larger than pi/4. + Computed using gmp. + */ + +const long double sincosl_table[] = { + +/* x = 1.48437500000000000000000000000000000e-01L 3ffc3000000000000000000000000000 */ +/* cos(x) = 0.fd2f5320e1b790209b4dda2f98f79caaa7b873aff1014b0fbc5243766d03cb006bc837c4358 */ + 9.89003367927322909016887196069562069e-01L, /* 3ffefa5ea641c36f2041369bb45f31ef */ + 2.15663692029265697782289400027743703e-35L, /* 3f8bcaaa7b873aff1014b0fbc5243767 */ +/* sin(x) = 0.25dc50bc95711d0d9787d108fd438cf5959ee0bfb7a1e36e8b1a112968f356657420e9cc9ea */ + 1.47892995873409608580026675734609314e-01L, /* 3ffc2ee285e4ab88e86cbc3e8847ea1c */ + 9.74950446464233268291647449768590886e-36L, /* 3f8a9eb2b3dc17f6f43c6dd16342252d */ + +/* x = 1.56250000000000000000000000000000000e-01 3ffc4000000000000000000000000000 */ +/* cos(x) = 0.fce1a053e621438b6d60c76e8c45bf0a9dc71aa16f922acc10e95144ec796a249813c9cb649 */ + 9.87817783816471944100503034363211317e-01L, /* 3ffef9c340a7cc428716dac18edd188b */ + 4.74271307836705897892468107620526395e-35L, /* 3f8cf854ee38d50b7c915660874a8a27 */ +/* sin(x) = 0.27d66258bacd96a3eb335b365c87d59438c5142bb56a489e9b8db9d36234ffdebb6bdc22d8e */ + 1.55614992773556041209920643203516258e-01L, /* 3ffc3eb312c5d66cb51f599ad9b2e43f */ + -7.83989563419287980121718050629497270e-36L, /* bf8a4d78e75d7a8952b6ec2c8e48c594 */ + +/* x = 1.64062500000000000000000000000000000e-01 3ffc5000000000000000000000000000 */ +/* cos(x) = 0.fc8ffa01ba6807417e05962b0d9fdf1fddb0cc4c07d22e19e08019bffa50a6c7acdb40307a3 */ + 9.86571908399497588757337407495308409e-01L, /* 3ffef91ff40374d00e82fc0b2c561b40 */ + -2.47327949936985362476252401212720725e-35L, /* bf8c070112799d9fc16e8f30fbff3200 */ +/* sin(x) = 0.29cfd49b8be4f665276cab01cbf0426934906c3dd105473b226e410b1450f62e53ff7c6cce1 */ + 1.63327491736612850846866172454354370e-01L, /* 3ffc4e7ea4dc5f27b3293b65580e5f82 */ + 1.81380344301155485770367902300754350e-36L, /* 3f88349a48361ee882a39d913720858a */ + +/* x = 1.71875000000000000000000000000000000e-01 3ffc6000000000000000000000000000 */ +/* cos(x) = 0.fc3a6170f767ac735d63d99a9d439e1db5e59d3ef153a4265d5855850ed82b536bf361b80e3 */ + 9.85265817718213816204294709759578994e-01L, /* 3ffef874c2e1eecf58e6bac7b3353a87 */ + 2.26568029505818066141517497778527952e-35L, /* 3f8be1db5e59d3ef153a4265d5855851 */ +/* sin(x) = 0.2bc89f9f424de5485de7ce03b2514952b9faf5648c3244d4736feb95dbb9da49f3b58a9253b */ + 1.71030022031395019281347969239834331e-01L, /* 3ffc5e44fcfa126f2a42ef3e701d928a */ + 7.01395875187487608875416030203241317e-36L, /* 3f8a2a573f5eac9186489a8e6dfd72bb */ + +/* x = 1.79687500000000000000000000000000000e-01 3ffc7000000000000000000000000000 */ +/* cos(x) = 0.fbe0d7f7fef11e70aa43b8abf4f6a457cea20c8f3f676b47781f9821bbe9ce04b3c7b981c0b */ + 9.83899591489663972178309351416487245e-01L, /* 3ffef7c1afeffde23ce154877157e9ed */ + 2.73414318948066207810486330723761265e-35L, /* 3f8c22be75106479fb3b5a3bc0fcc10e */ +/* sin(x) = 0.2dc0bb80b49a97ffb34e8dd1f8db9df7af47ed2dcf58b12c8e7827e048cae929da02c04ecac */ + 1.78722113535153659375356241864180724e-01L, /* 3ffc6e05dc05a4d4bffd9a746e8fc6dd */ + -1.52906926517265103202547561260594148e-36L, /* bf8804285c09691853a769b8c3ec0fdc */ + +/* x = 1.87500000000000000000000000000000000e-01 3ffc8000000000000000000000000000 */ +/* cos(x) = 0.fb835efcf670dd2ce6fe7924697eea13ea358867e9cdb3899b783f4f9f43aa5626e8b67b3bc */ + 9.82473313101255257487327683243622495e-01L, /* 3ffef706bdf9ece1ba59cdfcf248d2fe */ + -1.64924358891557584625463868014230342e-35L, /* bf8b5ec15ca779816324c766487c0b06 */ +/* sin(x) = 0.2fb8205f75e56a2b56a1c4792f856258769af396e0189ef72c05e4df59a6b00e4b44a6ea515 */ + 1.86403296762269884552379983103205261e-01L, /* 3ffc7dc102fbaf2b515ab50e23c97c2b */ + 1.76460304806826780010586715975331753e-36L, /* 3f882c3b4d79cb700c4f7b9602f26fad */ + +/* x = 1.95312500000000000000000000000000000e-01 3ffc9000000000000000000000000000 */ +/* cos(x) = 0.fb21f7f5c156696b00ac1fe28ac5fd76674a92b4df80d9c8a46c684399005deccc41386257c */ + 9.80987069605669190469329896435309665e-01L, /* 3ffef643efeb82acd2d601583fc5158c */ + -1.90899259410096419886996331536278461e-36L, /* bf8844cc5ab6a5903f931badc9cbde34 */ +/* sin(x) = 0.31aec65df552876f82ece9a2356713246eba6799983d7011b0b3698d6e1da919c15d57c30c1 */ + 1.94073102892909791156055200214145404e-01L, /* 3ffc8d7632efaa943b7c17674d11ab39 */ + -9.67304741051998267208945242944928999e-36L, /* bf8a9b7228b30cccf851fdc9e992ce52 */ + +/* x = 2.03125000000000000000000000000000000e-01 3ffca000000000000000000000000000 */ +/* cos(x) = 0.fabca467fb3cb8f1d069f01d8ea33ade5bfd68296ecd1cc9f7b7609bbcf3676e726c3301334 */ + 9.79440951715548359998530954502987493e-01L, /* 3ffef57948cff67971e3a0d3e03b1d46 */ + 4.42878056591560757066844797290067990e-35L, /* 3f8cd6f2dfeb414b7668e64fbdbb04de */ +/* sin(x) = 0.33a4a5a19d86246710f602c44df4fa513f4639ce938477aeeabb82e8e0a7ed583a188879fd4 */ + 2.01731063801638804725038151164000971e-01L, /* 3ffc9d252d0cec31233887b016226fa8 */ + -4.27513434754966978435151290617384120e-36L, /* bf896bb02e718c5b1ee21445511f45c8 */ + +/* x = 2.10937500000000000000000000000000000e-01 3ffcb000000000000000000000000000 */ +/* cos(x) = 0.fa5365e8f1d3ca27be1db5d76ae64d983d7470a4ab0f4ccf65a2b8c67a380df949953a09bc1 */ + 9.77835053797959793331971572944454549e-01L, /* 3ffef4a6cbd1e3a7944f7c3b6baed5cd */ + -3.79207422905180416937210853779192702e-35L, /* bf8c933e145c7adaa7859984d2ea39cc */ +/* sin(x) = 0.3599b652f40ec999df12a0a4c8561de159c98d4e54555de518b97f48886f715d8df5f4f093e */ + 2.09376712085993643711890752724881652e-01L, /* 3ffcaccdb297a0764ccef895052642b1 */ + -1.59470287344329449965314638482515925e-36L, /* bf880f531b3958d5d5510d73a3405bbc */ + +/* x = 2.18750000000000000000000000000000000e-01 3ffcc000000000000000000000000000 */ +/* cos(x) = 0.f9e63e1d9e8b6f6f2e296bae5b5ed9c11fd7fa2fe11e09fc7bde901abed24b6365e72f7db4e */ + 9.76169473868635276723989035435135534e-01L, /* 3ffef3cc7c3b3d16dede5c52d75cb6be */ + -2.87727974249481583047944860626985460e-35L, /* bf8c31f701402e80f70fb01c210b7f2a */ +/* sin(x) = 0.378df09db8c332ce0d2b53d865582e4526ea336c768f68c32b496c6d11c1cd241bb9f1da523 */ + 2.17009581095010156760578095826055396e-01L, /* 3ffcbc6f84edc6199670695a9ec32ac1 */ + 1.07356488794216831812829549198201194e-35L, /* 3f8ac8a4dd466d8ed1ed1865692d8da2 */ + +/* x = 2.26562500000000000000000000000000000e-01 3ffcd000000000000000000000000000 */ +/* cos(x) = 0.f9752eba9fff6b98842beadab054a932fb0f8d5b875ae63d6b2288d09b148921aeb6e52f61b */ + 9.74444313585988980349711056045434344e-01L, /* 3ffef2ea5d753ffed7310857d5b560a9 */ + 3.09947905955053419304514538592548333e-35L, /* 3f8c4997d87c6adc3ad731eb59144685 */ +/* sin(x) = 0.39814cb10513453cb97b21bc1ca6a337b150c21a675ab85503bc09a436a10ab1473934e20c8 */ + 2.24629204957705292350428549796424820e-01L, /* 3ffccc0a6588289a29e5cbd90de0e535 */ + 2.42061510849297469844695751870058679e-36L, /* 3f889bd8a8610d33ad5c2a81de04d21b */ + +/* x = 2.34375000000000000000000000000000000e-01 3ffce000000000000000000000000000 */ +/* cos(x) = 0.f90039843324f9b940416c1984b6cbed1fc733d97354d4265788a86150493ce657cae032674 */ + 9.72659678244912752670913058267565260e-01L, /* 3ffef20073086649f3728082d833096e */ + -3.91759231819314904966076958560252735e-35L, /* bf8ca09701c6613465595ecd43babcf5 */ +/* sin(x) = 0.3b73c2bf6b4b9f668ef9499c81f0d965087f1753fa64b086e58cb8470515c18c1412f8c2e02 */ + 2.32235118611511462413930877746235872e-01L, /* 3ffcdb9e15fb5a5cfb3477ca4ce40f87 */ + -4.96930483364191020075024624332928910e-36L, /* bf89a6bde03a2b0166d3de469cd1ee3f */ + +/* x = 2.42187500000000000000000000000000000e-01 3ffcf000000000000000000000000000 */ +/* cos(x) = 0.f887604e2c39dbb20e4ec5825059a789ffc95b275ad9954078ba8a28d3fcfe9cc2c1d49697b */ + 9.70815676770349462947490545785046027e-01L, /* 3ffef10ec09c5873b7641c9d8b04a0b3 */ + 2.97458820972393859125277682021202860e-35L, /* 3f8c3c4ffe4ad93ad6ccaa03c5d45147 */ +/* sin(x) = 0.3d654aff15cb457a0fca854698aba33039a8a40626609204472d9d40309b626eccc6dff0ffa */ + 2.39826857830661564441369251810886574e-01L, /* 3ffceb2a57f8ae5a2bd07e542a34c55d */ + 2.39867036569896287240938444445071448e-36L, /* 3f88981cd45203133049022396cea018 */ + +/* x = 2.50000000000000000000000000000000000e-01 3ffd0000000000000000000000000000 */ +/* cos(x) = 0.f80aa4fbef750ba783d33cb95f94f8a41426dbe79edc4a023ef9ec13c944551c0795b84fee1 */ + 9.68912421710644784144595449494189205e-01L, /* 3ffef01549f7deea174f07a67972bf2a */ + -5.53634706113461989398873287749326500e-36L, /* bf89d6faf649061848ed7f704184fb0e */ +/* sin(x) = 0.3f55dda9e62aed7513bd7b8e6a3d1635dd5676648d7db525898d7086af9330f03c7f285442a */ + 2.47403959254522929596848704849389203e-01L, /* 3ffcfaaeed4f31576ba89debdc7351e9 */ + -7.36487001108599532943597115275811618e-36L, /* bf8a39445531336e50495b4ece51ef2a */ + +/* x = 2.57812500000000000000000000000000000e-01 3ffd0800000000000000000000000000 */ +/* cos(x) = 0.f78a098069792daabc9ee42591b7c5a68cb1ab822aeb446b3311b4ba5371b8970e2c1547ad7 */ + 9.66950029230677822008341623610531503e-01L, /* 3ffeef141300d2f25b55793dc84b2370 */ + -4.38972214432792412062088059990480514e-35L, /* bf8cd2cb9a72a3eea8a5dca667725a2d */ +/* sin(x) = 0.414572fd94556e6473d620271388dd47c0ba050cdb5270112e3e370e8c4705ae006426fb5d5 */ + 2.54965960415878467487556574864872628e-01L, /* 3ffd0515cbf65155b991cf58809c4e23 */ + 2.20280377918534721005071688328074154e-35L, /* 3f8bd47c0ba050cdb5270112e3e370e9 */ + +/* x = 2.65625000000000000000000000000000000e-01 3ffd1000000000000000000000000000 */ +/* cos(x) = 0.f7058fde0788dfc805b8fe88789e4f4253e3c50afe8b22f41159620ab5940ff7df9557c0d1f */ + 9.64928619104771009581074665315748371e-01L, /* 3ffeee0b1fbc0f11bf900b71fd10f13d */ + -3.66685832670820775002475545602761113e-35L, /* bf8c85ed60e1d7a80ba6e85f7534efaa */ +/* sin(x) = 0.4334033bcd90d6604f5f36c1d4b84451a87150438275b77470b50e5b968fa7962b5ffb379b7 */ + 2.62512399769153281450949626395692931e-01L, /* 3ffd0cd00cef364359813d7cdb0752e1 */ + 3.24923677072031064673177178571821843e-36L, /* 3f89146a1c5410e09d6ddd1c2d4396e6 */ + +/* x = 2.73437500000000000000000000000000000e-01 3ffd1800000000000000000000000000 */ +/* cos(x) = 0.f67d3a26af7d07aa4bd6d42af8c0067fefb96d5b46c031eff53627f215ea3242edc3f2e13eb */ + 9.62848314709379699899701093480214365e-01L, /* 3ffeecfa744d5efa0f5497ada855f180 */ + 4.88986966383343450799422013051821394e-36L, /* 3f899ffbee5b56d1b00c7bfd4d89fc85 */ +/* sin(x) = 0.452186aa5377ab20bbf2524f52e3a06a969f47166ab88cf88c111ad12c55941021ef3317a1a */ + 2.70042816718585031552755063618827102e-01L, /* 3ffd14861aa94ddeac82efc9493d4b8f */ + -2.37608892440611310321138680065803162e-35L, /* bf8bf956960b8e99547730773eee52ed */ + +/* x = 2.81250000000000000000000000000000000e-01 3ffd2000000000000000000000000000 */ +/* cos(x) = 0.f5f10a7bb77d3dfa0c1da8b57842783280d01ce3c0f82bae3b9d623c168d2e7c29977994451 */ + 9.60709243015561903066659350581313472e-01L, /* 3ffeebe214f76efa7bf4183b516af085 */ + -5.87011558231583960712013351601221840e-36L, /* bf89f35fcbf8c70fc1f5147118a770fa */ +/* sin(x) = 0.470df5931ae1d946076fe0dcff47fe31bb2ede618ebc607821f8462b639e1f4298b5ae87fd3 */ + 2.77556751646336325922023446828128568e-01L, /* 3ffd1c37d64c6b8765181dbf8373fd20 */ + -1.35848595468998128214344668770082997e-36L, /* bf87ce44d1219e71439f87de07b9d49c */ + +/* x = 2.89062500000000000000000000000000000e-01 3ffd2800000000000000000000000000 */ +/* cos(x) = 0.f561030ddd7a78960ea9f4a32c6521554995667f5547bafee9ec48b3155cdb0f7fd00509713 */ + 9.58511534581228627301969408154919822e-01L, /* 3ffeeac2061bbaf4f12c1d53e94658ca */ + 2.50770779371636481145735089393154404e-35L, /* 3f8c0aaa4cab33faaa3dd7f74f624599 */ +/* sin(x) = 0.48f948446abcd6b0f7fccb100e7a1b26eccad880b0d24b59948c7cdd49514d44b933e6985c2 */ + 2.85053745940547424587763033323252561e-01L, /* 3ffd23e52111aaf35ac3dff32c4039e8 */ + 2.04269325885902918802700123680403749e-35L, /* 3f8bb26eccad880b0d24b59948c7cdd5 */ + +/* x = 2.96875000000000000000000000000000000e-01 3ffd3000000000000000000000000000 */ +/* cos(x) = 0.f4cd261d3e6c15bb369c8758630d2ac00b7ace2a51c0631bfeb39ed158ba924cc91e259c195 */ + 9.56255323543175296975599942263028361e-01L, /* 3ffee99a4c3a7cd82b766d390eb0c61a */ + 3.21616572190865997051103645135837207e-35L, /* 3f8c56005bd671528e0318dff59cf68b */ +/* sin(x) = 0.4ae37710fad27c8aa9c4cf96c03519b9ce07dc08a1471775499f05c29f86190aaebaeb9716e */ + 2.92533342023327543624702326493913423e-01L, /* 3ffd2b8ddc43eb49f22aa7133e5b00d4 */ + 1.93539408668704450308003687950685128e-35L, /* 3f8b9b9ce07dc08a1471775499f05c2a */ + +/* x = 3.04687500000000000000000000000000000e-01 3ffd3800000000000000000000000000 */ +/* cos(x) = 0.f43575f94d4f6b272f5fb76b14d2a64ab52df1ee8ddf7c651034e5b2889305a9ea9015d758a */ + 9.53940747608894733981324795987611623e-01L, /* 3ffee86aebf29a9ed64e5ebf6ed629a5 */ + 2.88075689052478602008395972924657164e-35L, /* 3f8c3255a96f8f746efbe32881a72d94 */ +/* sin(x) = 0.4ccc7a50127e1de0cb6b40c302c651f7bded4f9e7702b0471ae0288d091a37391950907202f */ + 2.99995083378683051163248282011699944e-01L, /* 3ffd3331e94049f877832dad030c0b19 */ + 1.35174265535697850139283361475571050e-35L, /* 3f8b1f7bded4f9e7702b0471ae0288d1 */ + +/* x = 3.12500000000000000000000000000000000e-01 3ffd4000000000000000000000000000 */ +/* cos(x) = 0.f399f500c9e9fd37ae9957263dab8877102beb569f101ee4495350868e5847d181d50d3cca2 */ + 9.51567948048172202145488217364270962e-01L, /* 3ffee733ea0193d3fa6f5d32ae4c7b57 */ + 6.36842628598115658308749288799884606e-36L, /* 3f8a0ee2057d6ad3e203dc892a6a10d2 */ +/* sin(x) = 0.4eb44a5da74f600207aaa090f0734e288603ffadb3eb2542a46977b105f8547128036dcf7f0 */ + 3.07438514580380850670502958201982091e-01L, /* 3ffd3ad129769d3d80081eaa8243c1cd */ + 1.06515172423204645839241099453417152e-35L, /* 3f8ac510c07ff5b67d64a8548d2ef621 */ + +/* x = 3.20312500000000000000000000000000000e-01 3ffd4800000000000000000000000000 */ +/* cos(x) = 0.f2faa5a1b74e82fd61fa05f9177380e8e69b7b15a945e8e5ae1124bf3d12b0617e03af4fab5 */ + 9.49137069684463027665847421762105623e-01L, /* 3ffee5f54b436e9d05fac3f40bf22ee7 */ + 6.84433965991637152250309190468859701e-37L, /* 3f86d1cd36f62b528bd1cb5c22497e7a */ +/* sin(x) = 0.509adf9a7b9a5a0f638a8fa3a60a199418859f18b37169a644fdb986c21ecb00133853bc35b */ + 3.14863181319745250865036315126939016e-01L, /* 3ffd426b7e69ee69683d8e2a3e8e9828 */ + 1.92431240212432926993057705062834160e-35L, /* 3f8b99418859f18b37169a644fdb986c */ + +/* x = 3.28125000000000000000000000000000000e-01 3ffd5000000000000000000000000000 */ +/* cos(x) = 0.f2578a595224dd2e6bfa2eb2f99cc674f5ea6f479eae2eb580186897ae3f893df1113ca06b8 */ + 9.46648260886053321846099507295532976e-01L, /* 3ffee4af14b2a449ba5cd7f45d65f33a */ + -4.32906339663000890941529420498824645e-35L, /* bf8ccc5850ac85c30a8e8a53ff3cbb43 */ +/* sin(x) = 0.5280326c3cf481823ba6bb08eac82c2093f2bce3c4eb4ee3dec7df41c92c8a4226098616075 */ + 3.22268630433386625687745919893188031e-01L, /* 3ffd4a00c9b0f3d20608ee9aec23ab21 */ + -1.49505897804759263483853908335500228e-35L, /* bf8b3df6c0d431c3b14b11c213820be3 */ + +/* x = 3.35937500000000000000000000000000000e-01 3ffd5800000000000000000000000000 */ +/* cos(x) = 0.f1b0a5b406b526d886c55feadc8d0dcc8eb9ae2ac707051771b48e05b25b000009660bdb3e3 */ + 9.44101673557004345630017691253124860e-01L, /* 3ffee3614b680d6a4db10d8abfd5b91a */ + 1.03812535240120229609822461172145584e-35L, /* 3f8ab991d735c558e0e0a2ee3691c0b6 */ +/* sin(x) = 0.54643b3da29de9b357155eef0f332fb3e66c83bf4dddd9491c5eb8e103ccd92d6175220ed51 */ + 3.29654409930860171914317725126463176e-01L, /* 3ffd5190ecf68a77a6cd5c557bbc3ccd */ + -1.22606996784743214973082192294232854e-35L, /* bf8b04c19937c40b22226b6e3a1471f0 */ + +/* x = 3.43750000000000000000000000000000000e-01 3ffd6000000000000000000000000000 */ +/* cos(x) = 0.f105fa4d66b607a67d44e042725204435142ac8ad54dfb0907a4f6b56b06d98ee60f19e557a */ + 9.41497463127881068644511236053670815e-01L, /* 3ffee20bf49acd6c0f4cfa89c084e4a4 */ + 3.20709366603165602071590241054884900e-36L, /* 3f8910d450ab22b5537ec241e93dad5b */ +/* sin(x) = 0.5646f27e8bd65cbe3a5d61ff06572290ee826d9674a00246b05ae26753cdfc90d9ce81a7d02 */ + 3.37020069022253076261281754173810024e-01L, /* 3ffd591bc9fa2f5972f8e97587fc195d */ + -2.21435756148839473677777545049890664e-35L, /* bf8bd6f117d92698b5ffdb94fa51d98b */ + +/* x = 3.51562500000000000000000000000000000e-01 3ffd6800000000000000000000000000 */ +/* cos(x) = 0.f0578ad01ede707fa39c09dc6b984afef74f3dc8d0efb0f4c5a6b13771145b3e0446fe33887 */ + 9.38835788546265488632578305984712554e-01L, /* 3ffee0af15a03dbce0ff473813b8d731 */ + -3.98758068773974031348585072752245458e-35L, /* bf8ca808458611b978827859d2ca7644 */ +/* sin(x) = 0.582850a41e1dd46c7f602ea244cdbbbfcdfa8f3189be794dda427ce090b5f85164f1f80ac13 */ + 3.44365158145698408207172046472223747e-01L, /* 3ffd60a14290787751b1fd80ba891337 */ + -3.19791885005480924937758467594051927e-36L, /* bf89100c815c339d9061ac896f60c7dc */ + +/* x = 3.59375000000000000000000000000000000e-01 3ffd7000000000000000000000000000 */ +/* cos(x) = 0.efa559f5ec3aec3a4eb03319278a2d41fcf9189462261125fe6147b078f1daa0b06750a1654 */ + 9.36116812267055290294237411019508588e-01L, /* 3ffedf4ab3ebd875d8749d6066324f14 */ + 3.40481591236710658435409862439032162e-35L, /* 3f8c6a0fe7c8c4a31130892ff30a3d84 */ +/* sin(x) = 0.5a084e28e35fda2776dfdbbb5531d74ced2b5d17c0b1afc4647529d50c295e36d8ceec126c1 */ + 3.51689228994814059222584896955547016e-01L, /* 3ffd682138a38d7f689ddb7f6eed54c7 */ + 1.75293433418270210567525412802083294e-35L, /* 3f8b74ced2b5d17c0b1afc4647529d51 */ + +/* x = 3.67187500000000000000000000000000000e-01 3ffd7800000000000000000000000000 */ +/* cos(x) = 0.eeef6a879146af0bf9b95ea2ea0ac0d3e2e4d7e15d93f48cbd41bf8e4fded40bef69e19eafa */ + 9.33340700242548435655299229469995527e-01L, /* 3ffeddded50f228d5e17f372bd45d416 */ + -4.75255707251679831124800898831382223e-35L, /* bf8cf960e8d940f513605b9a15f2038e */ +/* sin(x) = 0.5be6e38ce8095542bc14ee9da0d36483e6734bcab2e07624188af5653f114eeb46738fa899d */ + 3.58991834546065053677710299152868941e-01L, /* 3ffd6f9b8e33a025550af053ba76834e */ + -2.06772389262723368139416970257112089e-35L, /* bf8bb7c198cb4354d1f89dbe7750a9ac */ + +/* x = 3.75000000000000000000000000000000000e-01 3ffd8000000000000000000000000000 */ +/* cos(x) = 0.ee35bf5ccac89052cd91ddb734d3a47e262e3b609db604e217053803be0091e76daf28a89b7 */ + 9.30507621912314291149476792229555481e-01L, /* 3ffedc6b7eb9959120a59b23bb6e69a7 */ + 2.74541088551732982573335285685416092e-35L, /* 3f8c23f13171db04edb02710b829c01e */ +/* sin(x) = 0.5dc40955d9084f48a94675a2498de5d851320ff5528a6afb3f2e24de240fce6cbed1ba0ccd6 */ + 3.66272529086047561372909351716264177e-01L, /* 3ffd7710255764213d22a519d6892638 */ + -1.96768433534936592675897818253108989e-35L, /* bf8ba27aecdf00aad759504c0d1db21e */ + +/* x = 3.82812500000000000000000000000000000e-01 3ffd8800000000000000000000000000 */ +/* cos(x) = 0.ed785b5c44741b4493c56bcb9d338a151c6f6b85d8f8aca658b28572c162b199680eb9304da */ + 9.27617750192851909628030798799961350e-01L, /* 3ffedaf0b6b888e83689278ad7973a67 */ + 7.58520371916345756281201167126854712e-36L, /* 3f8a42a38ded70bb1f1594cb1650ae58 */ +/* sin(x) = 0.5f9fb80f21b53649c432540a50e22c53057ff42ae0fdf1307760dc0093f99c8efeb2fbd7073 */ + 3.73530868238692946416839752660848112e-01L, /* 3ffd7e7ee03c86d4d92710c950294389 */ + -1.48023494778986556048879113411517128e-35L, /* bf8b3acfa800bd51f020ecf889f23ff7 */ + +/* x = 3.90625000000000000000000000000000000e-01 3ffd9000000000000000000000000000 */ +/* cos(x) = 0.ecb7417b8d4ee3fec37aba4073aa48f1f14666006fb431d9671303c8100d10190ec8179c41d */ + 9.24671261467036098502113014560138771e-01L, /* 3ffed96e82f71a9dc7fd86f57480e755 */ + -4.14187124860031825108649347251175815e-35L, /* bf8cb87075cccffc825e7134c767e1bf */ +/* sin(x) = 0.6179e84a09a5258a40e9b5face03e525f8b5753cd0105d93fe6298010c3458e84d75fe420e9 */ + 3.80766408992390192057200703388896675e-01L, /* 3ffd85e7a1282694962903a6d7eb3810 */ + -2.02009541175208636336924533372496107e-35L, /* bf8bada074a8ac32fefa26c019d67fef */ + +/* x = 3.98437500000000000000000000000000000e-01 3ffd9800000000000000000000000000 */ +/* cos(x) = 0.ebf274bf0bda4f62447e56a093626798d3013b5942b1abfd155aacc9dc5c6d0806a20d6b9c1 */ + 9.21668335573351918175411368202712714e-01L, /* 3ffed7e4e97e17b49ec488fcad4126c5 */ + -1.83587995433957622948710263541479322e-35L, /* bf8b8672cfec4a6bd4e5402eaa553362 */ +/* sin(x) = 0.6352929dd264bd44a02ea766325d8aa8bd9695fc8def3caefba5b94c9a3c873f7b2d3776ead */ + 3.87978709727025046051079690813741960e-01L, /* 3ffd8d4a4a774992f51280ba9d98c976 */ + 8.01904783870935075844443278617586301e-36L, /* 3f8a5517b2d2bf91bde795df74b72993 */ + +/* x = 4.06250000000000000000000000000000000e-01 3ffda000000000000000000000000000 */ +/* cos(x) = 0.eb29f839f201fd13b93796827916a78f15c85230a4e8ea4b21558265a14367e1abb4c30695a */ + 9.18609155794918267837824977718549863e-01L, /* 3ffed653f073e403fa27726f2d04f22d */ + 2.97608282778274433460057745798409849e-35L, /* 3f8c3c78ae429185274752590aac132d */ +/* sin(x) = 0.6529afa7d51b129631ec197c0a840a11d7dc5368b0a47956feb285caa8371c4637ef17ef01b */ + 3.95167330240934236244832640419653657e-01L, /* 3ffd94a6be9f546c4a58c7b065f02a10 */ + 7.57560031388312550940040194042627704e-36L, /* 3f8a423afb8a6d16148f2adfd650b955 */ + +/* x = 4.14062500000000000000000000000000000e-01 3ffda800000000000000000000000000 */ +/* cos(x) = 0.ea5dcf0e30cf03e6976ef0b1ec26515fba47383855c3b4055a99b5e86824b2cd1a691fdca7b */ + 9.15493908848301228563917732180221882e-01L, /* 3ffed4bb9e1c619e07cd2edde163d84d */ + -3.50775517955306954815090901168305659e-35L, /* bf8c75022dc63e3d51e25fd52b3250bd */ +/* sin(x) = 0.66ff380ba0144109e39a320b0a3fa5fd65ea0585bcbf9b1a769a9b0334576c658139e1a1cbe */ + 4.02331831777773111217105598880982387e-01L, /* 3ffd9bfce02e805104278e68c82c28ff */ + -1.95678722882848174723569916504871563e-35L, /* bf8ba029a15fa7a434064e5896564fcd */ + +/* x = 4.21875000000000000000000000000000000e-01 3ffdb000000000000000000000000000 */ +/* cos(x) = 0.e98dfc6c6be031e60dd3089cbdd18a75b1f6b2c1e97f79225202f03dbea45b07a5ec4efc062 */ + 9.12322784872117846492029542047341734e-01L, /* 3ffed31bf8d8d7c063cc1ba611397ba3 */ + 7.86903886556373674267948132178845568e-36L, /* 3f8a4eb63ed6583d2fef244a405e07b8 */ +/* sin(x) = 0.68d32473143327973bc712bcc4ccddc47630d755850c0655243b205934dc49ffed8eb76adcb */ + 4.09471777053295066122694027011452236e-01L, /* 3ffda34c91cc50cc9e5cef1c4af31333 */ + 2.23945241468457597921655785729821354e-35L, /* 3f8bdc47630d755850c0655243b20593 */ + +/* x = 4.29687500000000000000000000000000000e-01 3ffdb800000000000000000000000000 */ +/* cos(x) = 0.e8ba8393eca7821aa563d83491b6101189b3b101c3677f73d7bad7c10f9ee02b7ab4009739a */ + 9.09095977415431051650381735684476417e-01L, /* 3ffed1750727d94f04354ac7b069236c */ + 1.20886014028444155733776025085677953e-35L, /* 3f8b01189b3b101c3677f73d7bad7c11 */ +/* sin(x) = 0.6aa56d8e8249db4eb60a761fe3f9e559be456b9e13349ca99b0bfb787f22b95db3b70179615 */ + 4.16586730282041119259112448831069657e-01L, /* 3ffdaa95b63a09276d3ad829d87f8fe8 */ + -2.00488106831998813675438269796963612e-35L, /* bf8baa641ba9461eccb635664f404878 */ + +/* x = 4.37500000000000000000000000000000000e-01 3ffdc000000000000000000000000000 */ +/* cos(x) = 0.e7e367d2956cfb16b6aa11e5419cd0057f5c132a6455bf064297e6a76fe2b72bb630d6d50ff */ + 9.05813683425936420744516660652700258e-01L, /* 3ffecfc6cfa52ad9f62d6d5423ca833a */ + -3.60950307605941169775676563004467163e-35L, /* bf8c7fd4051f66acdd5207cdeb40cac5 */ +/* sin(x) = 0.6c760c14c8585a51dbd34660ae6c52ac7036a0b40887a0b63724f8b4414348c3063a637f457 */ + 4.23676257203938010361683988031102480e-01L, /* 3ffdb1d83053216169476f4d1982b9b1 */ + 1.40484456388654470329473096579312595e-35L, /* 3f8b2ac7036a0b40887a0b63724f8b44 */ + +/* x = 4.45312500000000000000000000000000000e-01 3ffdc800000000000000000000000000 */ +/* cos(x) = 0.e708ac84d4172a3e2737662213429e14021074d7e702e77d72a8f1101a7e70410df8273e9aa */ + 9.02476103237941504925183272675895999e-01L, /* 3ffece115909a82e547c4e6ecc442685 */ + 2.26282899501344419018306295680210602e-35L, /* 3f8be14021074d7e702e77d72a8f1102 */ +/* sin(x) = 0.6e44f8c36eb10a1c752d093c00f4d47ba446ac4c215d26b0316442f168459e677d06e7249e3 */ + 4.30739925110803197216321517850849190e-01L, /* 3ffdb913e30dbac42871d4b424f003d3 */ + 1.54096780001629398850891218396761548e-35L, /* 3f8b47ba446ac4c215d26b0316442f17 */ + +/* x = 4.53125000000000000000000000000000000e-01 3ffdd000000000000000000000000000 */ +/* cos(x) = 0.e62a551594b970a770b15d41d4c0e483e47aca550111df6966f9e7ac3a94ae49e6a71eb031e */ + 8.99083440560138456216544929209379307e-01L, /* 3ffecc54aa2b2972e14ee162ba83a982 */ + -2.06772615490904370666670275154751976e-35L, /* bf8bb7c1b8535aafeee209699061853c */ +/* sin(x) = 0.70122c5ec5028c8cff33abf4fd340ccc382e038379b09cf04f9a52692b10b72586060cbb001 */ + 4.37777302872755132861618974702796680e-01L, /* 3ffdc048b17b140a3233fcceafd3f4d0 */ + 9.62794364503442612477117426033922467e-36L, /* 3f8a998705c0706f36139e09f34a4d25 */ + +/* x = 4.60937500000000000000000000000000000e-01 3ffdd800000000000000000000000000 */ +/* cos(x) = 0.e54864fe33e8575cabf5bd0e5cf1b1a8bc7c0d5f61702450fa6b6539735820dd2603ae355d5 */ + 8.95635902463170698900570000446256350e-01L, /* 3ffeca90c9fc67d0aeb957eb7a1cb9e3 */ + 3.73593741659866883088620495542311808e-35L, /* 3f8c8d45e3e06afb0b812287d35b29cc */ +/* sin(x) = 0.71dd9fb1ff4677853acb970a9f6729c6e3aac247b1c57cea66c77413f1f98e8b9e98e49d851 */ + 4.44787960964527211433056012529525211e-01L, /* 3ffdc7767ec7fd19de14eb2e5c2a7d9d */ + -1.67187936511493678007508371613954899e-35L, /* bf8b6391c553db84e3a831599388bec1 */ + +/* x = 4.68750000000000000000000000000000000e-01 3ffde000000000000000000000000000 */ +/* cos(x) = 0.e462dfc670d421ab3d1a15901228f146a0547011202bf5ab01f914431859aef577966bc4fa4 */ + 8.92133699366994404723900253723788575e-01L, /* 3ffec8c5bf8ce1a843567a342b202452 */ + -1.10771937602567314732693079264692504e-35L, /* bf8ad72bf571fddbfa814a9fc0dd779d */ +/* sin(x) = 0.73a74b8f52947b681baf6928eb3fb021769bf4779bad0e3aa9b1cdb75ec60aad9fc63ff19d5 */ + 4.51771471491683776581688750134062870e-01L, /* 3ffdce9d2e3d4a51eda06ebda4a3acff */ + -1.19387223016472295893794387275284505e-35L, /* bf8afbd12c81710c8a5e38aac9c64914 */ + +/* x = 4.76562500000000000000000000000000000e-01 3ffde800000000000000000000000000 */ +/* cos(x) = 0.e379c9045f29d517c4808aa497c2057b2b3d109e76c0dc302d4d0698b36e3f0bdbf33d8e952 */ + 8.88577045028035543317609023116020980e-01L, /* 3ffec6f39208be53aa2f890115492f84 */ + 4.12354278954664731443813655177022170e-36L, /* 3f895ecacf44279db0370c0b5341a62d */ +/* sin(x) = 0.756f28d011d98528a44a75fc29c779bd734ecdfb582fdb74b68a4c4c4be54cfd0b2d3ad292f */ + 4.58727408216736592377295028972874773e-01L, /* 3ffdd5bca340476614a29129d7f0a71e */ + -4.70946994194182908929251719575431779e-36L, /* bf8990a32c4c8129f40922d25d6ceced */ + +/* x = 4.84375000000000000000000000000000000e-01 3ffdf000000000000000000000000000 */ +/* cos(x) = 0.e28d245c58baef72225e232abc003c4366acd9eb4fc2808c2ab7fe7676cf512ac7f945ae5fb */ + 8.84966156526143291697296536966647926e-01L, /* 3ffec51a48b8b175dee444bc46557800 */ + 4.53370570288325630442037826313462165e-35L, /* 3f8ce21b3566cf5a7e14046155bff3b4 */ +/* sin(x) = 0.77353054ca72690d4c6e171fd99e6b39fa8e1ede5f052fd2964534c75340970a3a9cd3c5c32 */ + 4.65655346585160182681199512507546779e-01L, /* 3ffddcd4c15329c9a43531b85c7f667a */ + -1.56282598978971872478619772155305961e-35L, /* bf8b4c60571e121a0fad02d69bacb38b */ + +/* x = 4.92187500000000000000000000000000000e-01 3ffdf800000000000000000000000000 */ +/* cos(x) = 0.e19cf580eeec046aa1422fa74807ecefb2a1911c94e7b5f20a00f70022d940193691e5bd790 */ + 8.81301254251340599140161908298100173e-01L, /* 3ffec339eb01ddd808d542845f4e9010 */ + -1.43419192312116687783945619009629445e-35L, /* bf8b3104d5e6ee36b184a0df5ff08ffe */ +/* sin(x) = 0.78f95b0560a9a3bd6df7bd981dc38c61224d08bc20631ea932e605e53b579e9e0767dfcbbcb */ + 4.72554863751304451146551317808516942e-01L, /* 3ffde3e56c1582a68ef5b7def660770e */ + 9.31324774957768018850224267625371204e-36L, /* 3f8a8c2449a117840c63d5265cc0bca7 */ + +/* x = 5.00000000000000000000000000000000000e-01 3ffe0000000000000000000000000000 */ +/* cos(x) = 0.e0a94032dbea7cedbddd9da2fafad98556566b3a89f43eabd72350af3e8b19e801204d8fe2e */ + 8.77582561890372716116281582603829681e-01L, /* 3ffec1528065b7d4f9db7bbb3b45f5f6 */ + -2.89484960181363924855192538540698851e-35L, /* bf8c33d54d4ca62bb05e0aa146e57a86 */ +/* sin(x) = 0.7abba1d12c17bfa1d92f0d93f60ded9992f45b4fcaf13cd58b303693d2a0db47db35ae8a3a9 */ + 4.79425538604203000273287935215571402e-01L, /* 3ffdeaee8744b05efe8764bc364fd838 */ + -1.38426977616718318950175848639381926e-35L, /* bf8b2666d0ba4b0350ec32a74cfc96c3 */ + +/* x = 5.07812500000000000000000000000000000e-01 3ffe0400000000000000000000000000 */ +/* cos(x) = 0.dfb20840f3a9b36f7ae2c515342890b5ec583b8366cc2b55029e95094d31112383f2553498b */ + 8.73810306413054508282556837071377159e-01L, /* 3ffebf641081e75366def5c58a2a6851 */ + 1.25716864497849302237218128599994785e-35L, /* 3f8b0b5ec583b8366cc2b55029e95095 */ +/* sin(x) = 0.7c7bfdaf13e5ed17212f8a7525bfb113aba6c0741b5362bb8d59282a850b63716bca0c910f0 */ + 4.86266951793275574311011306895834993e-01L, /* 3ffdf1eff6bc4f97b45c84be29d496ff */ + -1.12269393250914752644352376448094271e-35L, /* bf8add8a8b27f17c9593a88e54dafaaf */ + +/* x = 5.15625000000000000000000000000000000e-01 3ffe0800000000000000000000000000 */ +/* cos(x) = 0.deb7518814a7a931bbcc88c109cd41c50bf8bb48f20ae8c36628d1d3d57574f7dc58f27d91c */ + 8.69984718058417388828915599901466243e-01L, /* 3ffebd6ea310294f526377991182139b */ + -4.68168638300575626782741319792183837e-35L, /* bf8cf1d7a03a25b86fa8b9e4ceb97161 */ +/* sin(x) = 0.7e3a679daaf25c676542bcb4028d0964172961c921823a4ef0c3a9070d886dbd073f6283699 */ + 4.93078685753923057265136552753487121e-01L, /* 3ffdf8e99e76abc9719d950af2d00a34 */ + 7.06498693112535056352301101088624950e-36L, /* 3f8a2c82e52c3924304749de187520e2 */ + +/* x = 5.23437500000000000000000000000000000e-01 3ffe0c00000000000000000000000000 */ +/* cos(x) = 0.ddb91ff318799172bd2452d0a3889f5169c64a0094bcf0b8aa7dcf0d7640a2eba68955a80be */ + 8.66106030320656714696616831654267220e-01L, /* 3ffebb723fe630f322e57a48a5a14711 */ + 2.35610597588322493119667003904687628e-35L, /* 3f8bf5169c64a0094bcf0b8aa7dcf0d7 */ +/* sin(x) = 0.7ff6d8a34bd5e8fa54c97482db5159df1f24e8038419c0b448b9eea8939b5d4dfcf40900257 */ + 4.99860324733013463819556536946425724e-01L, /* 3ffdffdb628d2f57a3e95325d20b6d45 */ + 1.94636052312235297538564591686645139e-35L, /* 3f8b9df1f24e8038419c0b448b9eea89 */ + +/* x = 5.31250000000000000000000000000000000e-01 3ffe1000000000000000000000000000 */ +/* cos(x) = 0.dcb7777ac420705168f31e3eb780ce9c939ecada62843b54522f5407eb7f21e556059fcd734 */ + 8.62174479934880504367162510253324274e-01L, /* 3ffeb96eeef58840e0a2d1e63c7d6f02 */ + -3.71556818317533582234562471835771823e-35L, /* bf8c8b1b6309a92cebde255d6e855fc1 */ +/* sin(x) = 0.81b149ce34caa5a4e650f8d09fd4d6aa74206c32ca951a93074c83b2d294d25dbb0f7fdfad2 */ + 5.06611454814257367642296000893867192e-01L, /* 3ffe0362939c69954b49cca1f1a13faa */ + -3.10963699824274155702706043065967062e-35L, /* bf8c4aac5efc9e69ab572b67c59be269 */ + +/* x = 5.39062500000000000000000000000000000e-01 3ffe1400000000000000000000000000 */ +/* cos(x) = 0.dbb25c25b8260c14f6e7bc98ec991b70c65335198b0ab628bad20cc7b229d4dd62183cfa055 */ + 8.58190306862660347046629564970494649e-01L, /* 3ffeb764b84b704c1829edcf7931d932 */ + 2.06439574601190798155563653000684861e-35L, /* 3f8bb70c65335198b0ab628bad20cc7b */ +/* sin(x) = 0.8369b434a372da7eb5c8a71fe36ce1e0b2b493f6f5cb2e38bcaec2a556b3678c401940d1c3c */ + 5.13331663943471218288801270215706878e-01L, /* 3ffe06d3686946e5b4fd6b914e3fc6da */ + -2.26614796466671970772244932848067224e-35L, /* bf8be1f4d4b6c090a34d1c743513d5ab */ + +/* x = 5.46875000000000000000000000000000000e-01 3ffe1800000000000000000000000000 */ +/* cos(x) = 0.daa9d20860827063fde51c09e855e9932e1b17143e7244fd267a899d41ae1f3bc6a0ec42e27 */ + 8.54153754277385385143451785105103176e-01L, /* 3ffeb553a410c104e0c7fbca3813d0ac */ + -1.68707534013095152873222061722573172e-35L, /* bf8b66cd1e4e8ebc18dbb02d9857662c */ +/* sin(x) = 0.852010f4f0800521378bd8dd614753d080c2e9e0775ffc609947b9132f5357404f464f06a58 */ + 5.20020541953727004760213699874674730e-01L, /* 3ffe0a4021e9e1000a426f17b1bac28f */ + -3.32415021330884924833711842866896734e-35L, /* bf8c617bf9e8b0fc45001cfb35c23767 */ + +/* x = 5.54687500000000000000000000000000000e-01 3ffe1c00000000000000000000000000 */ +/* cos(x) = 0.d99ddd44e44a43d4d4a3a3ed95204106fd54d78e8c7684545c0da0b7c2c72be7a89b7c182ad */ + 8.50065068549420263957072899177793617e-01L, /* 3ffeb33bba89c89487a9a94747db2a41 */ + -4.73753917078785974356016104842568442e-35L, /* bf8cf7c81559438b9c4bdd5d1f92fa42 */ +/* sin(x) = 0.86d45935ab396cb4e421e822dee54f3562dfcefeaa782184c23401d231f5ad981a1cc195b18 */ + 5.26677680590386730710789410624833901e-01L, /* 3ffe0da8b26b5672d969c843d045bdcb */ + -3.67066148195515214077582496518566735e-35L, /* bf8c8654e901880aac3ef3d9ee5ff16e */ + +/* x = 5.62500000000000000000000000000000000e-01 3ffe2000000000000000000000000000 */ +/* cos(x) = 0.d88e820b1526311dd561efbc0c1a9a5375eb26f65d246c5744b13ca26a7e0fd42556da843c8 */ + 8.45924499231067954459723078597493262e-01L, /* 3ffeb11d04162a4c623baac3df781835 */ + 1.98054947141989878179164342925274053e-35L, /* 3f8ba5375eb26f65d246c5744b13ca27 */ +/* sin(x) = 0.88868625b4e1dbb2313310133022527200c143a5cb16637cb7daf8ade82459ff2e98511f40f */ + 5.33302673536020173329131103308161529e-01L, /* 3ffe110d0c4b69c3b764626620266045 */ + -3.42715291319551615996993795226755157e-35L, /* bf8c6c6ff9f5e2d1a74ce41a41283a91 */ + +/* x = 5.70312500000000000000000000000000000e-01 3ffe2400000000000000000000000000 */ +/* cos(x) = 0.d77bc4985e93a607c9d868b906bbc6bbe3a04258814acb0358468b826fc91bd4d814827f65e */ + 8.41732299041338366963111794309701085e-01L, /* 3ffeaef78930bd274c0f93b0d1720d78 */ + -4.30821936750410026005408345400225948e-35L, /* bf8cca20e2fded3bf5a9a7e53dcba3ed */ +/* sin(x) = 0.8a3690fc5bfc11bf9535e2739a8512f448a41251514bbed7fc18d530f9b4650fcbb2861b0aa */ + 5.39895116435204405041660709903993340e-01L, /* 3ffe146d21f8b7f8237f2a6bc4e7350a */ + 1.42595803521626714477253741404712093e-35L, /* 3f8b2f448a41251514bbed7fc18d5310 */ + +/* x = 5.78125000000000000000000000000000000e-01 3ffe2800000000000000000000000000 */ +/* cos(x) = 0.d665a937b4ef2b1f6d51bad6d988a4419c1d7051faf31a9efa151d7631117efac03713f950a */ + 8.37488723850523685315353348917240617e-01L, /* 3ffeaccb526f69de563edaa375adb311 */ + 2.72761997872084533045777718677326179e-35L, /* 3f8c220ce0eb828fd798d4f7d0a8ebb2 */ +/* sin(x) = 0.8be472f9776d809af2b88171243d63d66dfceeeb739cc894e023fbc165a0e3f26ff729c5d57 */ + 5.46454606919203564403349553749411001e-01L, /* 3ffe17c8e5f2eedb0135e57102e2487b */ + -2.11870230730160315420936523771864858e-35L, /* bf8bc29920311148c63376b1fdc043ea */ + +/* x = 5.85937500000000000000000000000000000e-01 3ffe2c00000000000000000000000000 */ +/* cos(x) = 0.d54c3441844897fc8f853f0655f1ba695eba9fbfd7439dbb1171d862d9d9146ca5136f825ac */ + 8.33194032664581363070224042208032321e-01L, /* 3ffeaa98688308912ff91f0a7e0cabe3 */ + 4.39440050052045486567668031751259899e-35L, /* 3f8cd34af5d4fdfeba1cedd88b8ec317 */ +/* sin(x) = 0.8d902565817ee7839bce3cd128060119492cd36d42d82ada30d7f8bde91324808377ddbf5d4 */ + 5.52980744630527369849695082681623667e-01L, /* 3ffe1b204acb02fdcf07379c79a2500c */ + 8.26624790417342895897164123189984127e-37L, /* 3f8719492cd36d42d82ada30d7f8bde9 */ + +/* x = 5.93750000000000000000000000000000000e-01 3ffe3000000000000000000000000000 */ +/* cos(x) = 0.d42f6a1b9f0168cdf031c2f63c8d9304d86f8d34cb1d5fccb68ca0f2241427fc18d1fd5bbdf */ + 8.28848487609325734810171790119116638e-01L, /* 3ffea85ed4373e02d19be06385ec791b */ + 1.43082508100496581719048175506239770e-35L, /* 3f8b304d86f8d34cb1d5fccb68ca0f22 */ +/* sin(x) = 0.8f39a191b2ba6122a3fa4f41d5a3ffd421417d46f19a22230a14f7fcc8fce5c75b4b28b29d1 */ + 5.59473131247366877384844006003116688e-01L, /* 3ffe1e7343236574c24547f49e83ab48 */ + -1.28922620524163922306886952100992796e-37L, /* bf845ef5f415c8732eeee7af584019b8 */ + +/* x = 6.01562500000000000000000000000000000e-01 3ffe3400000000000000000000000000 */ +/* cos(x) = 0.d30f4f392c357ab0661c5fa8a7d9b26627846fef214b1d19a22379ff9eddba087cf410eb097 */ + 8.24452353914429207485643598212356053e-01L, /* 3ffea61e9e72586af560cc38bf514fb3 */ + 3.79160239225080026987031418939026741e-35L, /* 3f8c93313c237f790a58e8cd111bcffd */ +/* sin(x) = 0.90e0e0d81ca678796cc92c8ea8c2815bc72ca78abe571bfa8576aacc571e096a33237e0e830 */ + 5.65931370507905990773159095689276114e-01L, /* 3ffe21c1c1b0394cf0f2d992591d5185 */ + 1.02202775968053982310991962521535027e-36L, /* 3f875bc72ca78abe571bfa8576aacc57 */ + +/* x = 6.09375000000000000000000000000000000e-01 3ffe3800000000000000000000000000 */ +/* cos(x) = 0.d1ebe81a95ee752e48a26bcd32d6e922d7eb44b8ad2232f6930795e84b56317269b9dd1dfa6 */ + 8.20005899897234008255550633876556043e-01L, /* 3ffea3d7d0352bdcea5c9144d79a65ae */ + -1.72008811955230823416724332297991247e-35L, /* bf8b6dd2814bb4752ddcd096cf86a17b */ +/* sin(x) = 0.9285dc9bc45dd9ea3d02457bcce59c4175aab6ff7929a8d287195525fdace200dba032874fb */ + 5.72355068234507240384953706824503608e-01L, /* 3ffe250bb93788bbb3d47a048af799cb */ + 2.12572273479933123944580199464514529e-35L, /* 3f8bc4175aab6ff7929a8d2871955260 */ + +/* x = 6.17187500000000000000000000000000000e-01 3ffe3c00000000000000000000000000 */ +/* cos(x) = 0.d0c5394d772228195e25736c03574707de0af1ca344b13bd3914bfe27518e9e426f5deff1e1 */ + 8.15509396946375476876345384201386217e-01L, /* 3ffea18a729aee445032bc4ae6d806af */ + -4.28589138410712954051679139949341961e-35L, /* bf8cc7c10fa871ae5da76216375a00ec */ +/* sin(x) = 0.94288e48bd0335fc41c4cbd2920497a8f5d1d8185c99fa0081f90c27e2a53ffdd208a0dbe69 */ + 5.78743832357770354521111378581385347e-01L, /* 3ffe28511c917a066bf8838997a52409 */ + 1.77998063432551282609698670002456093e-35L, /* 3f8b7a8f5d1d8185c99fa0081f90c27e */ + +/* x = 6.25000000000000000000000000000000000e-01 3ffe4000000000000000000000000000 */ +/* cos(x) = 0.cf9b476c897c25c5bfe750dd3f308eaf7bcc1ed00179a256870f4200445043dcdb1974b5878 */ + 8.10963119505217902189534803941080724e-01L, /* 3ffe9f368ed912f84b8b7fcea1ba7e61 */ + 1.10481292856794436426051402418804358e-35L, /* 3f8ad5ef7983da002f344ad0e1e84009 */ +/* sin(x) = 0.95c8ef544210ec0b91c49bd2aa09e8515fa61a156ebb10f5f8c232a6445b61ebf3c2ec268f9 */ + 5.85097272940462154805399314150080459e-01L, /* 3ffe2b91dea88421d817238937a55414 */ + -1.78164576278056195136525335403380464e-35L, /* bf8b7aea059e5ea9144ef0a073dcd59c */ + +/* x = 6.32812500000000000000000000000000000e-01 3ffe4400000000000000000000000000 */ +/* cos(x) = 0.ce6e171f92f2e27f32225327ec440ddaefae248413efc0e58ceee1ae369aabe73f88c87ed1a */ + 8.06367345055103913698795406077297399e-01L, /* 3ffe9cdc2e3f25e5c4fe6444a64fd888 */ + 1.04235088143133625463876245029180850e-35L, /* 3f8abb5df5c490827df81cb19ddc35c7 */ +/* sin(x) = 0.9766f93cd18413a6aafc1cfc6fc28abb6817bf94ce349901ae3f48c3215d3eb60acc5f78903 */ + 5.91415002201316315087000225758031236e-01L, /* 3ffe2ecdf279a308274d55f839f8df85 */ + 8.07390238063560077355762466502569603e-36L, /* 3f8a576d02f7f299c6932035c7e91864 */ + +/* x = 6.40625000000000000000000000000000000e-01 3ffe4800000000000000000000000000 */ +/* cos(x) = 0.cd3dad1b5328a2e459f993f4f5108819faccbc4eeba9604e81c7adad51cc8a2561631a06826 */ + 8.01722354098418450607492605652964208e-01L, /* 3ffe9a7b5a36a65145c8b3f327e9ea21 */ + 6.09487851305233089325627939458963741e-36L, /* 3f8a033f599789dd752c09d038f5b5aa */ +/* sin(x) = 0.9902a58a45e27bed68412b426b675ed503f54d14c8172e0d373f42cadf04daf67319a7f94be */ + 5.97696634538701531238647618967334337e-01L, /* 3ffe32054b148bc4f7dad0825684d6cf */ + -2.49527608940873714527427941350461554e-35L, /* bf8c0957e0559759bf468f964605e9a9 */ + +/* x = 6.48437500000000000000000000000000000e-01 3ffe4c00000000000000000000000000 */ +/* cos(x) = 0.cc0a0e21709883a3ff00911e11a07ee3bd7ea2b04e081be99be0264791170761ae64b8b744a */ + 7.97028430141468342004642741431945296e-01L, /* 3ffe98141c42e1310747fe01223c2341 */ + -8.35364432831812599727083251866305534e-37L, /* bf871c42815d4fb1f7e416641fd9b86f */ +/* sin(x) = 0.9a9bedcdf01b38d993f3d7820781de292033ead73b89e28f39313dbe3a6e463f845b5fa8490 */ + 6.03941786554156657267270287527367726e-01L, /* 3ffe3537db9be03671b327e7af040f04 */ + -2.54578992328947177770363936132309779e-35L, /* bf8c0eb6fe60a94623b0eb863676120e */ + +/* x = 6.56250000000000000000000000000000000e-01 3ffe5000000000000000000000000000 */ +/* cos(x) = 0.cad33f00658fe5e8204bbc0f3a66a0e6a773f87987a780b243d7be83b3db1448ca0e0e62787 */ + 7.92285859677178543141501323781709399e-01L, /* 3ffe95a67e00cb1fcbd04097781e74cd */ + 2.47519558228473167879248891673807645e-35L, /* 3f8c07353b9fc3cc3d3c05921ebdf41e */ +/* sin(x) = 0.9c32cba2b14156ef05256c4f857991ca6a547cd7ceb1ac8a8e62a282bd7b9183648a462bd04 */ + 6.10150077075791371273742393566183220e-01L, /* 3ffe386597456282adde0a4ad89f0af3 */ + 1.33842237929938963780969418369150532e-35L, /* 3f8b1ca6a547cd7ceb1ac8a8e62a282c */ + +/* x = 6.64062500000000000000000000000000000e-01 3ffe5400000000000000000000000000 */ +/* cos(x) = 0.c99944936cf48c8911ff93fe64b3ddb7981e414bdaf6aae1203577de44878c62bc3bc9cf7b9 */ + 7.87494932167606083931328295965533034e-01L, /* 3ffe93328926d9e9191223ff27fcc968 */ + -2.57915385618070637156514241185180920e-35L, /* bf8c12433f0df5a1284aa8f6fe54410e */ +/* sin(x) = 0.9dc738ad14204e689ac582d0f85826590feece34886cfefe2e08cf2bb8488d55424dc9d3525 */ + 6.16321127181550943005700433761731837e-01L, /* 3ffe3b8e715a28409cd1358b05a1f0b0 */ + 2.88497530050197716298085892460478666e-35L, /* 3f8c32c87f7671a44367f7f17046795e */ + +/* x = 6.71875000000000000000000000000000000e-01 3ffe5800000000000000000000000000 */ +/* cos(x) = 0.c85c23c26ed7b6f014ef546c47929682122876bfbf157de0aff3c4247d820c746e32cd4174f */ + 7.82655940026272796930787447428139026e-01L, /* 3ffe90b84784ddaf6de029dea8d88f25 */ + 1.69332045679237919427807771288506254e-35L, /* 3f8b682122876bfbf157de0aff3c4248 */ +/* sin(x) = 0.9f592e9b66a9cf906a3c7aa3c10199849040c45ec3f0a747597311038101780c5f266059dbf */ + 6.22454560222343683041926705090443330e-01L, /* 3ffe3eb25d36cd539f20d478f5478203 */ + 1.91974786921147072717621236192269859e-35L, /* 3f8b9849040c45ec3f0a747597311038 */ + +/* x = 6.79687500000000000000000000000000000e-01 3ffe5c00000000000000000000000000 */ +/* cos(x) = 0.c71be181ecd6875ce2da5615a03cca207d9adcb9dfb0a1d6c40a4f0056437f1a59ccddd06ee */ + 7.77769178600317903122203513685412863e-01L, /* 3ffe8e37c303d9ad0eb9c5b4ac2b407a */ + -4.05296033424632846931240580239929672e-35L, /* bf8caefc13291a31027af149dfad87fd */ +/* sin(x) = 0.a0e8a725d33c828c11fa50fd9e9a15ffecfad43f3e534358076b9b0f6865694842b1e8c67dc */ + 6.28550001845029662028004327939032867e-01L, /* 3ffe41d14e4ba679051823f4a1fb3d34 */ + 1.65507421184028099672784511397428852e-35L, /* 3f8b5ffecfad43f3e534358076b9b0f7 */ + +/* x = 6.87500000000000000000000000000000000e-01 3ffe6000000000000000000000000000 */ +/* cos(x) = 0.c5d882d2ee48030c7c07d28e981e34804f82ed4cf93655d2365389b716de6ad44676a1cc5da */ + 7.72834946152471544810851845913425178e-01L, /* 3ffe8bb105a5dc900618f80fa51d303c */ + 3.94975229341211664237241534741146939e-35L, /* 3f8ca4027c176a67c9b2ae91b29c4db9 */ +/* sin(x) = 0.a2759c0e79c35582527c32b55f5405c182c66160cb1d9eb7bb0b7cdf4ad66f317bda4332914 */ + 6.34607080015269296850309914203671436e-01L, /* 3ffe44eb381cf386ab04a4f8656abea8 */ + 4.33025916939968369326060156455927002e-36L, /* 3f897060b1985832c767adeec2df37d3 */ + +/* x = 6.95312500000000000000000000000000000e-01 3ffe6400000000000000000000000000 */ +/* cos(x) = 0.c4920cc2ec38fb891b38827db08884fc66371ac4c2052ca8885b981bbcfd3bb7b093ee31515 */ + 7.67853543842850365879920759114193964e-01L, /* 3ffe89241985d871f712367104fb6111 */ + 3.75100035267325597157244776081706979e-36L, /* 3f893f198dc6b130814b2a2216e606ef */ +/* sin(x) = 0.a400072188acf49cd6b173825e038346f105e1301afe642bcc364cea455e21e506e3e927ed8 */ + 6.40625425040230409188409779413961021e-01L, /* 3ffe48000e431159e939ad62e704bc07 */ + 2.46542747294664049615806500747173281e-36L, /* 3f88a37882f0980d7f3215e61b267523 */ + +/* x = 7.03125000000000000000000000000000000e-01 3ffe6800000000000000000000000000 */ +/* cos(x) = 0.c348846bbd3631338ffe2bfe9dd1381a35b4e9c0c51b4c13fe376bad1bf5caacc4542be0aa9 */ + 7.62825275710576250507098753625429792e-01L, /* 3ffe869108d77a6c62671ffc57fd3ba2 */ + 4.22067411888601505004748939382325080e-35L, /* 3f8cc0d1ada74e0628da609ff1bb5d69 */ +/* sin(x) = 0.a587e23555bb08086d02b9c662cdd29316c3e9bd08d93793634a21b1810cce73bdb97a99b9e */ + 6.46604669591152370524042159882800763e-01L, /* 3ffe4b0fc46aab761010da05738cc59c */ + -3.41742981816219412415674365946079826e-35L, /* bf8c6b6749e0b217b9364364e5aef274 */ + +/* x = 7.10937500000000000000000000000000000e-01 3ffe6c00000000000000000000000000 */ +/* cos(x) = 0.c1fbeef380e4ffdd5a613ec8722f643ffe814ec2343e53adb549627224fdc9f2a7b77d3d69f */ + 7.57750448655219342240234832230493361e-01L, /* 3ffe83f7dde701c9ffbab4c27d90e45f */ + -2.08767968311222650582659938787920125e-35L, /* bf8bbc0017eb13dcbc1ac524ab69d8de */ +/* sin(x) = 0.a70d272a76a8d4b6da0ec90712bb748b96dabf88c3079246f3db7eea6e58ead4ed0e2843303 */ + 6.52544448725765956407573982284767763e-01L, /* 3ffe4e1a4e54ed51a96db41d920e2577 */ + -8.61758060284379660697102362141557170e-36L, /* bf8a6e8d24a80ee79f0db721849022b2 */ + +/* x = 7.18750000000000000000000000000000000e-01 3ffe7000000000000000000000000000 */ +/* cos(x) = 0.c0ac518c8b6ae710ba37a3eeb90cb15aebcb8bed4356fb507a48a6e97de9aa6d9660116b436 */ + 7.52629372418066476054541324847143116e-01L, /* 3ffe8158a31916d5ce21746f47dd7219 */ + 3.71306958657663189665450864311104571e-35L, /* 3f8c8ad75e5c5f6a1ab7da83d245374c */ +/* sin(x) = 0.a88fcfebd9a8dd47e2f3c76ef9e2439920f7e7fbe735f8bcc985491ec6f12a2d4214f8cfa99 */ + 6.58444399910567541589583954884041989e-01L, /* 3ffe511f9fd7b351ba8fc5e78eddf3c5 */ + -4.54412944084300330523721391865787219e-35L, /* bf8ce336f840c020c6503a19b3d5b70a */ + +/* x = 7.26562500000000000000000000000000000e-01 3ffe7400000000000000000000000000 */ +/* cos(x) = 0.bf59b17550a4406875969296567cf3e3b4e483061877c02811c6cae85fad5a6c3da58f49292 */ + 7.47462359563216166669700384714767552e-01L, /* 3ffe7eb362eaa14880d0eb2d252cacfa */ + -9.11094340926220027288083639048016945e-36L, /* bf8a8389636f9f3cf107fafdc726a2f4 */ +/* sin(x) = 0.aa0fd66eddb921232c28520d3911b8a03193b47f187f1471ac216fbcd5bb81029294d3a73f1 */ + 6.64304163042946276515506587432846246e-01L, /* 3ffe541facddbb7242465850a41a7223 */ + 4.26004843895378210155889028714676019e-35L, /* 3f8cc5018c9da3f8c3f8a38d610b7de7 */ + +/* x = 7.34375000000000000000000000000000000e-01 3ffe7800000000000000000000000000 */ +/* cos(x) = 0.be0413f84f2a771c614946a88cbf4da1d75a5560243de8f2283fefa0ea4a48468a52d51d8b3 */ + 7.42249725458501306991347253449610537e-01L, /* 3ffe7c0827f09e54ee38c2928d51197f */ + -3.78925270049800913539923473871287550e-35L, /* bf8c92f1452d54fede10b86ebe0082f9 */ +/* sin(x) = 0.ab8d34b36acd987210ed343ec65d7e3adc2e7109fce43d55c8d57dfdf55b9e01d2cc1f1b9ec */ + 6.70123380473162894654531583500648495e-01L, /* 3ffe571a6966d59b30e421da687d8cbb */ + -1.33165852952743729897634069393684656e-36L, /* bf87c523d18ef6031bc2aa372a82020b */ + +/* x = 7.42187500000000000000000000000000000e-01 3ffe7c00000000000000000000000000 */ +/* cos(x) = 0.bcab7e6bfb2a14a9b122c574a376bec98ab14808c64a4e731b34047e217611013ac99c0f25d */ + 7.36991788256240741057089385586450844e-01L, /* 3ffe7956fcd7f654295362458ae946ed */ + 4.72358938637974850573747497460125519e-35L, /* 3f8cf64c558a404632527398d9a023f1 */ +/* sin(x) = 0.ad07e4c409d08c4fa3a9057bb0ac24b8636e74e76f51e09bd6b2319707cbd9f5e254643897a */ + 6.75901697026178809189642203142423973e-01L, /* 3ffe5a0fc98813a1189f47520af76158 */ + 2.76252586616364878801928456702948857e-35L, /* 3f8c25c31b73a73b7a8f04deb5918cb8 */ + +/* x = 7.50000000000000000000000000000000000e-01 3ffe8000000000000000000000000000 */ +/* cos(x) = 0.bb4ff632a908f73ec151839cb9d993b4e0bfb8f20e7e44e6e4aee845e35575c3106dbe6fd06 */ + 7.31688868873820886311838753000084529e-01L, /* 3ffe769fec655211ee7d82a3073973b3 */ + 1.48255637548931697184991710293198620e-35L, /* 3f8b3b4e0bfb8f20e7e44e6e4aee845e */ +/* sin(x) = 0.ae7fe0b5fc786b2d966e1d6af140a488476747c2646425fc7533f532cd044cb10a971a49a6a */ + 6.81638760023334166733241952779893908e-01L, /* 3ffe5cffc16bf8f0d65b2cdc3ad5e281 */ + 2.74838775935027549024224114338667371e-35L, /* 3f8c24423b3a3e1323212fe3a99fa996 */ + +/* x = 7.57812500000000000000000000000000000e-01 3ffe8400000000000000000000000000 */ +/* cos(x) = 0.b9f180ba77dd0751628e135a9508299012230f14becacdd14c3f8862d122de5b56d55b53360 */ + 7.26341290974108590410147630237598973e-01L, /* 3ffe73e30174efba0ea2c51c26b52a10 */ + 3.12683579338351123545814364980658990e-35L, /* 3f8c4c80911878a5f6566e8a61fc4317 */ +/* sin(x) = 0.aff522a954f2ba16d9defdc416e33f5e9a5dfd5a6c228e0abc4d521327ff6e2517a7b3851dd */ + 6.87334219303873534951703613035647220e-01L, /* 3ffe5fea4552a9e5742db3bdfb882dc6 */ + 4.76739454455410744997012795035529128e-35L, /* 3f8cfaf4d2efead361147055e26a9099 */ + +/* x = 7.65625000000000000000000000000000000e-01 3ffe8800000000000000000000000000 */ +/* cos(x) = 0.b890237d3bb3c284b614a0539016bfa1053730bbdf940fa895e185f8e58884d3dda15e63371 */ + 7.20949380945696418043812784148447688e-01L, /* 3ffe712046fa776785096c2940a7202d */ + 4.78691285733673379499536326050811832e-35L, /* 3f8cfd0829b985defca07d44af0c2fc7 */ +/* sin(x) = 0.b167a4c90d63c4244cf5493b7cc23bd3c3c1225e078baa0c53d6d400b926281f537a1a260e6 */ + 6.92987727246317910281815490823048210e-01L, /* 3ffe62cf49921ac7884899ea9276f984 */ + 4.50089871077663557180849219529189918e-35L, /* 3f8cde9e1e0912f03c5d50629eb6a006 */ + +/* x = 7.73437500000000000000000000000000000e-01 3ffe8c00000000000000000000000000 */ +/* cos(x) = 0.b72be40067aaf2c050dbdb7a14c3d7d4f203f6b3f0224a4afe55d6ec8e92b508fd5c5984b3b */ + 7.15513467882981573520620561289896903e-01L, /* 3ffe6e57c800cf55e580a1b7b6f42988 */ + -3.02191815581445336509438104625489192e-35L, /* bf8c41586fe04a607eedada80d51489c */ +/* sin(x) = 0.b2d7614b1f3aaa24df2d6e20a77e1ca3e6d838c03e29c1bcb026e6733324815fadc9eb89674 */ + 6.98598938789681741301929277107891591e-01L, /* 3ffe65aec2963e755449be5adc414efc */ + 2.15465226809256290914423429408722521e-35L, /* 3f8bca3e6d838c03e29c1bcb026e6733 */ + +/* x = 7.81250000000000000000000000000000000e-01 3ffe9000000000000000000000000000 */ +/* cos(x) = 0.b5c4c7d4f7dae915ac786ccf4b1a498d3e73b6e5e74fe7519d9c53ee6d6b90e881bddfc33e1 */ + 7.10033883566079674974121643959490219e-01L, /* 3ffe6b898fa9efb5d22b58f0d99e9635 */ + -4.09623224763692443220896752907902465e-35L, /* bf8cb3960c6248d0c580c573131d608d */ +/* sin(x) = 0.b44452709a59752905913765434a59d111f0433eb2b133f7d103207e2aeb4aae111ddc385b3 */ + 7.04167511454533672780059509973942844e-01L, /* 3ffe6888a4e134b2ea520b226eca8695 */ + -2.87259372740393348676633610275598640e-35L, /* bf8c3177707de60a6a76604177e6fc0f */ + +/* x = 7.89062500000000000000000000000000000e-01 3ffe9400000000000000000000000000 */ +/* cos(x) = 0.b45ad4975b1294cadca4cf40ec8f22a68cd14b175835239a37e63acb85e8e9505215df18140 */ + 7.04510962440574606164129481545916976e-01L, /* 3ffe68b5a92eb6252995b9499e81d91e */ + 2.60682037357042658395360726992048803e-35L, /* 3f8c1534668a58bac1a91cd1bf31d65c */ +/* sin(x) = 0.b5ae7285bc10cf515753847e8f8b7a30e0a580d929d770103509880680f7b8b0e8ad23b65d8 */ + 7.09693105363899724959669028139035515e-01L, /* 3ffe6b5ce50b78219ea2aea708fd1f17 */ + -4.37026016974122945368562319136420097e-36L, /* bf8973c7d69fc9b58a23fbf2bd9dfe60 */ +};
new file mode 100644 --- /dev/null +++ b/lib/sinl.c @@ -0,0 +1,101 @@ +/* s_sinl.c -- long double version of s_sin.c. + * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* sinl(x) + * Return sine function of x. + * + * kernel function: + * __kernel_sinl ... sine function on [-pi/4,pi/4] + * __kernel_cosl ... cose function on [-pi/4,pi/4] + * __ieee754_rem_pio2l ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ + +#include <math.h> + +#include "mathl.h" + +#include "trigl.h" +#include "trigl.c" +#include "sincosl.c" + +long double +sinl (long double x) +{ + long double y[2], z = 0.0L; + int n; + + /* |x| ~< pi/4 */ + if (x >= -0.7853981633974483096156608458198757210492 && + x <= 0.7853981633974483096156608458198757210492) + return kernel_sinl (x, z, 0); + + /* sinl(Inf or NaN) is NaN, sinl(0) is 0 */ + else if (x + x == x || x != x) + return x - x; /* NaN */ + + /* argument reduction needed */ + else + { + n = ieee754_rem_pio2l (x, y); + switch (n & 3) + { + case 0: + return kernel_sinl (y[0], y[1], 1); + case 1: + return kernel_cosl (y[0], y[1]); + case 2: + return -kernel_sinl (y[0], y[1], 1); + default: + return -kernel_cosl (y[0], y[1]); + } + } +} + +#if 0 +int +main () +{ + printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492)); + printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *29)); + printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *2)); + printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *30)); + printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *4)); + printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *32)); + printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *2/3)); + printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *4/3)); +} +#endif
new file mode 100644 --- /dev/null +++ b/lib/sqrtl.c @@ -0,0 +1,56 @@ +/* Emulation for sqrtl. + Contributed by Paolo Bonzini + + Copyright 2002, 2003 Free Software Foundation, Inc. + + This file is part of gnulib. + + gnulib is free software; you can redistribute it and/or modify it + under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation; either version 2.1, or (at your option) + any later version. + + gnulib is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public + License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with gnulib; see the file COPYING.LIB. If not, write to the Free + Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, + USA. + */ + +#include <float.h> +#include <math.h> + +#include "mathl.h" + +/* A simple Newton-Raphson method. */ +long double +sqrtl(long double x) +{ + long double delta, y; + int exponent; + + /* Check for negative numbers */ + if (x < 0.0L) + return (long double) sqrt(-1); + + /* Check for zero, NANs and infinites */ + if (x + x == x || x != x) + return x; + + frexpl (x, &exponent); + y = ldexpl (x, -exponent / 2); + + do + { + delta = y; + y = (y + x / y) * 0.5L; + delta -= y; + } + while (delta != 0.0L); + + return y; +}
new file mode 100644 --- /dev/null +++ b/lib/tanl.c @@ -0,0 +1,223 @@ +/* s_tanl.c -- long double version of s_tan.c. + * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. + */ + +/* @(#)s_tan.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* tanl(x) + * Return tangent function of x. + * + * kernel function: + * __kernel_tanl ... tangent function on [-pi/4,pi/4] + * __ieee754_rem_pio2l ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ + +#include <math.h> + +#include "mathl.h" + +#include "trigl.h" +#ifdef HAVE_SINL +#ifdef HAVE_COSL +#include "trigl.c" +#endif +#endif + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + Long double expansions contributed by + Stephen L. Moshier <moshier@na-net.ornl.gov> +*/ + +/* __kernel_tanl( x, y, k ) + * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854 + * Input x is assumed to be bounded by ~pi/4 in magnitude. + * Input y is the tail of x. + * Input k indicates whether tan (if k=1) or + * -1/tan (if k= -1) is returned. + * + * Algorithm + * 1. Since tan(-x) = -tan(x), we need only to consider positive x. + * 2. if x < 2^-57, return x with inexact if x!=0. + * 3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2) + * on [0,0.67433]. + * + * Note: tan(x+y) = tan(x) + tan'(x)*y + * ~ tan(x) + (1+x*x)*y + * Therefore, for better accuracy in computing tan(x+y), let + * r = x^3 * R(x^2) + * then + * tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y)) + * + * 4. For x in [0.67433,pi/4], let y = pi/4 - x, then + * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) + * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) + */ + + +static const long double + pio4hi = 7.8539816339744830961566084581987569936977E-1L, + pio4lo = 2.1679525325309452561992610065108379921906E-35L, + + /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2) + 0 <= x <= 0.6743316650390625 + Peak relative error 8.0e-36 */ + TH = 3.333333333333333333333333333333333333333E-1L, + T0 = -1.813014711743583437742363284336855889393E7L, + T1 = 1.320767960008972224312740075083259247618E6L, + T2 = -2.626775478255838182468651821863299023956E4L, + T3 = 1.764573356488504935415411383687150199315E2L, + T4 = -3.333267763822178690794678978979803526092E-1L, + + U0 = -1.359761033807687578306772463253710042010E8L, + U1 = 6.494370630656893175666729313065113194784E7L, + U2 = -4.180787672237927475505536849168729386782E6L, + U3 = 8.031643765106170040139966622980914621521E4L, + U4 = -5.323131271912475695157127875560667378597E2L; + /* 1.000000000000000000000000000000000000000E0 */ + + +long double +kernel_tanl (long double x, long double y, int iy) +{ + long double z, r, v, w, s, u, u1; + int flag, sign; + + sign = 1; + if (x < 0) + { + x = -x; + y = -y; + sign = -1; + } + + if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* x < 2**-57 */ + { + if ((int) x == 0) + { /* generate inexact */ + if (iy == -1 && x == 0.0) + return 1.0L / fabs (x); + else + return (iy == 1) ? x : -1.0L / x; + } + } + if (x >= 0.6743316650390625) /* |x| >= 0.6743316650390625 */ + { + flag = 1; + + z = pio4hi - x; + w = pio4lo - y; + x = z + w; + y = 0.0; + } + z = x * x; + r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4))); + v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z)))); + r = r / v; + + s = z * x; + r = y + z * (s * r + y); + r += TH * s; + w = x + r; + if (flag) + { + v = (long double) iy; + w = (v - 2.0 * (x - (w * w / (w + v) - r))); + if (sign < 0) + w = -w; + return w; + } + if (iy == 1) + return w; + else + { /* if allow error up to 2 ulp, + simply return -1.0/(x+r) here */ + /* compute -1.0/(x+r) accurately */ + u1 = (double) w; + v = r - (u1 - x); + z = -1.0 / w; + u = (double) z; + s = 1.0 + u * u1; + return u + z * (s + u * v); + } +} + +long double +tanl (long double x) +{ + long double y[2], z = 0.0L; + int n; + + /* |x| ~< pi/4 */ + if (x >= -0.7853981633974483096156608458198757210492 && + x <= 0.7853981633974483096156608458198757210492) + return kernel_tanl (x, z, 1); + + /* tanl(Inf or NaN) is NaN, tanl(0) is 0 */ + else if (x + x == x || x != x) + return x - x; /* NaN */ + + /* argument reduction needed */ + else + { + n = ieee754_rem_pio2l (x, y); + /* 1 -- n even, -1 -- n odd */ + return kernel_tanl (y[0], y[1], 1 - ((n & 1) << 1)); + } +} + +#if 0 +int +main () +{ + printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492)); + printf ("%.16Lg\n", tanl(-0.7853981633974483096156608458198757210492)); + printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 *3)); + printf ("%.16Lg\n", tanl(-0.7853981633974483096156608458198757210492 *31)); + printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 / 2)); + printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 * 3/2)); + printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 * 5/2)); +} +#endif
new file mode 100644 --- /dev/null +++ b/lib/trigl.c @@ -0,0 +1,636 @@ +/* Quad-precision floating point argument reduction. + Copyright (C) 1999 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Jakub Jelinek <jj@ultra.linux.cz> + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include <math.h> +#include <float.h> + +#include "mathl.h" + +/* Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi */ +static const int two_over_pi[] = { + 0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62, + 0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a, + 0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129, + 0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41, + 0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8, + 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf, + 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5, + 0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08, + 0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3, + 0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880, + 0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b, + 0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6, + 0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2, + 0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35, + 0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30, + 0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c, + 0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4, + 0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770, + 0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7, + 0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19, + 0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522, + 0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16, + 0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6, + 0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e, + 0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48, + 0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3, + 0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf, + 0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55, + 0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612, + 0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929, + 0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec, + 0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b, + 0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c, + 0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4, + 0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb, + 0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc, + 0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c, + 0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f, + 0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5, + 0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437, + 0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b, + 0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea, + 0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad, + 0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3, + 0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3, + 0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717, + 0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f, + 0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61, + 0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db, + 0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51, + 0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0, + 0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c, + 0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6, + 0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc, + 0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed, + 0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328, + 0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d, + 0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0, + 0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b, + 0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4, + 0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3, + 0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f, + 0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad, + 0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b, + 0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4, + 0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761, + 0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31, + 0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30, + 0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262, + 0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e, + 0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1, + 0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c, + 0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4, + 0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08, + 0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196, + 0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9, + 0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4, + 0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc, + 0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c, + 0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0, + 0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c, + 0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0, + 0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac, + 0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22, + 0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893, + 0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7, + 0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5, + 0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f, + 0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4, + 0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf, + 0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b, + 0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2, + 0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138, + 0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e, + 0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569, + 0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34, + 0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9, + 0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d, + 0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f, + 0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855, + 0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569, + 0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b, + 0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe, + 0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41, + 0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49, + 0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f, + 0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110, + 0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8, + 0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365, + 0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a, + 0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270, + 0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5, + 0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616, + 0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b, + 0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0, + 0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb, + 0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a, + 0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e, + 0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa, + 0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5, + 0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0, + 0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2, + 0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886, + 0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142, + 0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba, + 0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4, + 0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708, + 0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555, + 0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3, + 0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55, + 0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58, + 0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5, + 0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c, + 0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe, + 0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b, + 0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8, + 0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005, + 0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7, + 0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50, + 0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604, + 0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643, + 0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485, + 0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d, + 0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6, + 0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2, + 0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02, + 0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3, + 0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412, + 0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274, + 0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755, + 0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849, + 0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce, + 0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5, + 0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba, + 0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6, + 0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d, + 0x7b7b89, 0x483d38, +}; + +static const long double c[] = { +/* 93 bits of pi/2 */ +#define PI_2_1 c[0] + 1.57079632679489661923132169155131424e+00L, /* 3fff921fb54442d18469898cc5100000 */ + +/* pi/2 - PI_2_1 */ +#define PI_2_1t c[1] + 8.84372056613570112025531863263659260e-29L, /* 3fa1c06e0e68948127044533e63a0106 */ +}; + +static int kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec, + const int *ipio2); + +int +ieee754_rem_pio2l (long double x, long double *y) +{ + long double z, w, t; + double tx[8]; + int exp, n; + + if (x >= -0.78539816339744830961566084581987572104929234984377 + && x < 0.78539816339744830961566084581987572104929234984377) + /* x in <-pi/4, pi/4> */ + { + y[0] = x; + y[1] = 0; + return 0; + } + + if (x >= 2.35619449019234492884698253745962716314787704953131 + && x < 2.35619449019234492884698253745962716314787704953131) + if (x > 0) + { + /* 113 + 93 bit PI is ok */ + z = x - PI_2_1; + y[0] = z - PI_2_1t; + y[1] = (z - y[0]) - PI_2_1t; + return 1; + } + else + { + /* 113 + 93 bit PI is ok */ + z = x + PI_2_1; + y[0] = z + PI_2_1t; + y[1] = (z - y[0]) + PI_2_1t; + return -1; + } + + if (x + x == x || x != x) /* x is +=oo or NaN */ + { + y[0] = x - x; + y[1] = y[0]; + return 0; + } + + /* Handle large arguments. + We split the 113 bits of the mantissa into 5 24bit integers + stored in a double array. */ + z = frexp (x, &exp); + tx[0] = floorl (z *= 16777216.0); + z -= tx[0]; + tx[1] = floorl (z *= 16777216.0); + z -= tx[1]; + tx[2] = floorl (z *= 16777216.0); + z -= tx[2]; + tx[3] = floorl (z *= 16777216.0); + z -= tx[3]; + tx[4] = floorl (z *= 16777216.0); + + n = kernel_rem_pio2 (tx, tx + 5, exp - 24, tx[4] ? 5 : 4, 3, two_over_pi); + + /* The result is now stored in 3 double values, we need to convert it into + two long double values. */ + t = (long double) tx[6] + (long double) tx[7]; + w = (long double) tx[5]; + + if (x > 0) + { + y[0] = w + t; + y[1] = t - (y[0] - w); + return n; + } + else + { + y[0] = -(w + t); + y[1] = -t - (y[0] + w); + return -n; + } +} + +/* @(#)k_rem_pio2.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = + "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $"; +#endif + +/* + * kernel_rem_pio2(x,y,e0,nx,prec,ipio2) + * double x[],y[]; int e0,nx,prec; int ipio2[]; + * + * kernel_rem_pio2 return the last three digits of N with + * y = x - N*pi/2 + * so that |y| < pi/2. + * + * The method is to compute the integer (mod 8) and fraction parts of + * (2/pi)*x without doing the full multiplication. In general we + * skip the part of the product that are known to be a huge integer ( + * more accurately, = 0 mod 8 ). Thus the number of operations are + * independent of the exponent of the input. + * + * (2/pi) is represented by an array of 24-bit integers in ipio2[]. + * + * Input parameters: + * x[] The input value (must be positive) is broken into nx + * pieces of 24-bit integers in double precision format. + * x[i] will be the i-th 24 bit of x. The scaled exponent + * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 + * match x's up to 24 bits. + * + * Example of breaking a double positive z into x[0]+x[1]+x[2]: + * e0 = ilogb(z)-23 + * z = scalbn(z,-e0) + * for i = 0,1,2 + * x[i] = floor(z) + * z = (z-x[i])*2**24 + * + * + * y[] ouput result in an array of double precision numbers. + * The dimension of y[] is: + * 24-bit precision 1 + * 53-bit precision 2 + * 64-bit precision 2 + * 113-bit precision 3 + * The actual value is the sum of them. Thus for 113-bit + * precision, one may have to do something like: + * + * long double t,w,r_head, r_tail; + * t = (long double)y[2] + (long double)y[1]; + * w = (long double)y[0]; + * r_head = t+w; + * r_tail = w - (r_head - t); + * + * e0 The exponent of x[0] + * + * nx dimension of x[] + * + * prec an integer indicating the precision: + * 0 24 bits (single) + * 1 53 bits (double) + * 2 64 bits (extended) + * 3 113 bits (quad) + * + * ipio2[] + * integer array, contains the (24*i)-th to (24*i+23)-th + * bit of 2/pi after binary point. The corresponding + * floating value is + * + * ipio2[i] * 2^(-24(i+1)). + * + * External function: + * double scalbn(), floor(); + * + * + * Here is the description of some local variables: + * + * jk jk+1 is the initial number of terms of ipio2[] needed + * in the computation. The recommended value is 2,3,4, + * 6 for single, double, extended,and quad. + * + * jz local integer variable indicating the number of + * terms of ipio2[] used. + * + * jx nx - 1 + * + * jv index for pointing to the suitable ipio2[] for the + * computation. In general, we want + * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 + * is an integer. Thus + * e0-3-24*jv >= 0 or (e0-3)/24 >= jv + * Hence jv = max(0,(e0-3)/24). + * + * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. + * + * q[] double array with integral value, representing the + * 24-bits chunk of the product of x and 2/pi. + * + * q0 the corresponding exponent of q[0]. Note that the + * exponent for q[i] would be q0-24*i. + * + * PIo2[] double precision array, obtained by cutting pi/2 + * into 24 bits chunks. + * + * f[] ipio2[] in floating point + * + * iq[] integer array by breaking up q[] in 24-bits chunk. + * + * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] + * + * ih integer. If >0 it indicates q[] is >= 0.5, hence + * it also indicates the *sign* of the result. + * + */ + + +/* + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +static const int init_jk[] = { 2, 3, 4, 6 }; /* initial value for jk */ +static const double PIo2[] = { + 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ + 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ + 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ + 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */ + 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */ + 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */ + 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */ + 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ +}; + +static const double zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ + twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ + +int +kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec, + const int *ipio2) +{ + int jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih; + double z, fw, f[20], fq[20], q[20]; + + /* initialize jk */ + jk = init_jk[prec]; + jp = jk; + + /* determine jx,jv,q0, note that 3>q0 */ + jx = nx - 1; + jv = (e0 - 3) / 24; + if (jv < 0) + jv = 0; + q0 = e0 - 24 * (jv + 1); + + /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ + j = jv - jx; + m = jx + jk; + for (i = 0; i <= m; i++, j++) + f[i] = (j < 0) ? zero : (double) ipio2[j]; + + /* compute q[0],q[1],...q[jk] */ + for (i = 0; i <= jk; i++) + { + for (j = 0, fw = 0.0; j <= jx; j++) + fw += x[j] * f[jx + i - j]; + q[i] = fw; + } + + jz = jk; +recompute: + /* distill q[] into iq[] reversingly */ + for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) + { + fw = (double) ((int) (twon24 * z)); + iq[i] = (int) (z - two24 * fw); + z = q[j - 1] + fw; + } + + /* compute n */ + z = ldexp (z, q0); /* actual value of z */ + z -= 8.0 * floor (z * 0.125); /* trim off integer >= 8 */ + n = (int) z; + z -= (double) n; + ih = 0; + if (q0 > 0) + { /* need iq[jz-1] to determine n */ + i = (iq[jz - 1] >> (24 - q0)); + n += i; + iq[jz - 1] -= i << (24 - q0); + ih = iq[jz - 1] >> (23 - q0); + } + else if (q0 == 0) + ih = iq[jz - 1] >> 23; + else if (z >= 0.5) + ih = 2; + + if (ih > 0) + { /* q > 0.5 */ + n += 1; + carry = 0; + for (i = 0; i < jz; i++) + { /* compute 1-q */ + j = iq[i]; + if (carry == 0) + { + if (j != 0) + { + carry = 1; + iq[i] = 0x1000000 - j; + } + } + else + iq[i] = 0xffffff - j; + } + if (q0 > 0) + { /* rare case: chance is 1 in 12 */ + switch (q0) + { + case 1: + iq[jz - 1] &= 0x7fffff; + break; + case 2: + iq[jz - 1] &= 0x3fffff; + break; + } + } + if (ih == 2) + { + z = one - z; + if (carry != 0) + z -= ldexp (one, q0); + } + } + + /* check if recomputation is needed */ + if (z == zero) + { + j = 0; + for (i = jz - 1; i >= jk; i--) + j |= iq[i]; + if (j == 0) + { /* need recomputation */ + for (k = 1; iq[jk - k] == 0; k++); /* k = no. of terms needed */ + + for (i = jz + 1; i <= jz + k; i++) + { /* add q[jz+1] to q[jz+k] */ + f[jx + i] = (double) ipio2[jv + i]; + for (j = 0, fw = 0.0; j <= jx; j++) + fw += x[j] * f[jx + i - j]; + q[i] = fw; + } + jz += k; + goto recompute; + } + } + + /* chop off zero terms */ + if (z == 0.0) + { + jz -= 1; + q0 -= 24; + while (iq[jz] == 0) + { + jz--; + q0 -= 24; + } + } + else + { /* break z into 24-bit if necessary */ + z = ldexp (z, -q0); + if (z >= two24) + { + fw = (double) ((int) (twon24 * z)); + iq[jz] = (int) (z - two24 * fw); + jz += 1; + q0 += 24; + iq[jz] = (int) fw; + } + else + iq[jz] = (int) z; + } + + /* convert integer "bit" chunk to floating-point value */ + fw = ldexp (one, q0); + for (i = jz; i >= 0; i--) + { + q[i] = fw * (double) iq[i]; + fw *= twon24; + } + + /* compute PIo2[0,...,jp]*q[jz,...,0] */ + for (i = jz; i >= 0; i--) + { + for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++) + fw += PIo2[k] * q[i + k]; + fq[jz - i] = fw; + } + + /* compress fq[] into y[] */ + switch (prec) + { + case 0: + fw = 0.0; + for (i = jz; i >= 0; i--) + fw += fq[i]; + y[0] = (ih == 0) ? fw : -fw; + break; + case 1: + case 2: + fw = 0.0; + for (i = jz; i >= 0; i--) + fw += fq[i]; + y[0] = (ih == 0) ? fw : -fw; + fw = fq[0] - fw; + for (i = 1; i <= jz; i++) + fw += fq[i]; + y[1] = (ih == 0) ? fw : -fw; + break; + case 3: /* painful */ + for (i = jz; i > 0; i--) + { + fw = fq[i - 1] + fq[i]; + fq[i] += fq[i - 1] - fw; + fq[i - 1] = fw; + } + for (i = jz; i > 1; i--) + { + fw = fq[i - 1] + fq[i]; + fq[i] += fq[i - 1] - fw; + fq[i - 1] = fw; + } + for (fw = 0.0, i = jz; i >= 2; i--) + fw += fq[i]; + if (ih == 0) + { + y[0] = fq[0]; + y[1] = fq[1]; + y[2] = fw; + } + else + { + y[0] = -fq[0]; + y[1] = -fq[1]; + y[2] = -fw; + } + } + return n & 7; +}
new file mode 100644 --- /dev/null +++ b/lib/trigl.h @@ -0,0 +1,27 @@ +/* Declarations for sinl, cosl, tanl internal functions + Contributed by Paolo Bonzini + + Copyright 2002, 2003 Free Software Foundation, Inc. + + This file is part of gnulib. + + gnulib is free software; you can redistribute it and/or modify it + under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation; either version 2.1, or (at your option) + any later version. + + gnulib is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public + License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with gnulib; see the file COPYING.LIB. If not, write to the Free + Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, + USA. + */ + +extern int ieee754_rem_pio2l (long double x, long double *y); +extern long double kernel_sinl (long double x, long double y, int iy); +extern long double kernel_cosl (long double x, long double y); +
--- a/m4/ChangeLog +++ b/m4/ChangeLog @@ -1,3 +1,7 @@ +2003-02-18 Paolo Bonzino <bonzini@gnu.org> + + * mathl.m4: New file. + 2003-02-17 Bruno Haible <bruno@clisp.org> * mkdtemp.m4: New file, from GNU gettext with modifications.
new file mode 100644 --- /dev/null +++ b/m4/mathl.m4 @@ -0,0 +1,15 @@ +# mathl.m4 serial 1 +dnl Copyright (c) 2003 Free Software Foundation, Inc. +dnl This file is free software, distributed under the terms of the GNU +dnl General Public License. As a special exception to the GNU General +dnl Public License, this file may be distributed as part of a program +dnl that contains a configuration script generated by Autoconf, under +dnl the same distribution terms as the rest of the program. + +AC_DEFUN([gl_FUNC_LONG_DOUBLE_MATH], [ + +AC_CHECK_LIB(m, atan) +AC_REPLACE_FUNCS(floorl ceill sqrtl frexpl ldexpl asinl acosl atanl \ + logl expl tanl sinl cosl) + +])
new file mode 100644 --- /dev/null +++ b/modules/mathl @@ -0,0 +1,37 @@ +Description: +C99 functions for transcendental functions with long double arguments. + +Files: +lib/mathl.h +lib/acosl.c +lib/asinl.c +lib/atanl.c +lib/ceill.c +lib/cosl.c +lib/expl.c +lib/floorl.c +lib/frexpl.c +lib/ldexpl.c +lib/logl.c +lib/sincosl.c +lib/sinl.c +lib/sqrtl.c +lib/tanl.c +lib/trigl.c +lib/trigl.h +m4/mathl.m4 + +Depends-on: + +configure.ac: +gl_FUNC_LONG_DOUBLE_MATH + +Makefile.am: +noinst_HEADERS += mathl.h trigl.h trigl.c sincosl.c + +Include: +"mathl.h" + +Maintainer: +Paolo Bonzini <bonzini@gnu.org> +