diff options
Diffstat (limited to 'REORG.TODO/sysdeps/ieee754/ldbl-128/s_atanl.c')
| -rw-r--r-- | REORG.TODO/sysdeps/ieee754/ldbl-128/s_atanl.c | 253 |
1 files changed, 253 insertions, 0 deletions
diff --git a/REORG.TODO/sysdeps/ieee754/ldbl-128/s_atanl.c b/REORG.TODO/sysdeps/ieee754/ldbl-128/s_atanl.c new file mode 100644 index 0000000000..6f2cd549ec --- /dev/null +++ b/REORG.TODO/sysdeps/ieee754/ldbl-128/s_atanl.c @@ -0,0 +1,253 @@ +/* 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> + + This 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. + + This 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 this library; if not, see + <http://www.gnu.org/licenses/>. */ + + +#include <float.h> +#include <math.h> +#include <math_private.h> + +/* arctan(k/8), k = 0, ..., 82 */ +static const _Float128 atantbl[84] = { + L(0.0000000000000000000000000000000000000000E0), + L(1.2435499454676143503135484916387102557317E-1), /* arctan(0.125) */ + L(2.4497866312686415417208248121127581091414E-1), + L(3.5877067027057222039592006392646049977698E-1), + L(4.6364760900080611621425623146121440202854E-1), + L(5.5859931534356243597150821640166127034645E-1), + L(6.4350110879328438680280922871732263804151E-1), + L(7.1882999962162450541701415152590465395142E-1), + L(7.8539816339744830961566084581987572104929E-1), + L(8.4415398611317100251784414827164750652594E-1), + L(8.9605538457134395617480071802993782702458E-1), + L(9.4200004037946366473793717053459358607166E-1), + L(9.8279372324732906798571061101466601449688E-1), + L(1.0191413442663497346383429170230636487744E0), + L(1.0516502125483736674598673120862998296302E0), + L(1.0808390005411683108871567292171998202703E0), + L(1.1071487177940905030170654601785370400700E0), + L(1.1309537439791604464709335155363278047493E0), + L(1.1525719972156675180401498626127513797495E0), + L(1.1722738811284763866005949441337046149712E0), + L(1.1902899496825317329277337748293183376012E0), + L(1.2068173702852525303955115800565576303133E0), + L(1.2220253232109896370417417439225704908830E0), + L(1.2360594894780819419094519711090786987027E0), + L(1.2490457723982544258299170772810901230778E0), + L(1.2610933822524404193139408812473357720101E0), + L(1.2722973952087173412961937498224804940684E0), + L(1.2827408797442707473628852511364955306249E0), + L(1.2924966677897852679030914214070816845853E0), + L(1.3016288340091961438047858503666855921414E0), + L(1.3101939350475556342564376891719053122733E0), + L(1.3182420510168370498593302023271362531155E0), + L(1.3258176636680324650592392104284756311844E0), + L(1.3329603993374458675538498697331558093700E0), + L(1.3397056595989995393283037525895557411039E0), + L(1.3460851583802539310489409282517796256512E0), + L(1.3521273809209546571891479413898128509842E0), + L(1.3578579772154994751124898859640585287459E0), + L(1.3633001003596939542892985278250991189943E0), + L(1.3684746984165928776366381936948529556191E0), + L(1.3734007669450158608612719264449611486510E0), + L(1.3780955681325110444536609641291551522494E0), + L(1.3825748214901258580599674177685685125566E0), + L(1.3868528702577214543289381097042486034883E0), + L(1.3909428270024183486427686943836432060856E0), + L(1.3948567013423687823948122092044222644895E0), + L(1.3986055122719575950126700816114282335732E0), + L(1.4021993871854670105330304794336492676944E0), + L(1.4056476493802697809521934019958079881002E0), + L(1.4089588955564736949699075250792569287156E0), + L(1.4121410646084952153676136718584891599630E0), + L(1.4152014988178669079462550975833894394929E0), + L(1.4181469983996314594038603039700989523716E0), + L(1.4209838702219992566633046424614466661176E0), + L(1.4237179714064941189018190466107297503086E0), + L(1.4263547484202526397918060597281265695725E0), + L(1.4288992721907326964184700745371983590908E0), + L(1.4313562697035588982240194668401779312122E0), + L(1.4337301524847089866404719096698873648610E0), + L(1.4360250423171655234964275337155008780675E0), + L(1.4382447944982225979614042479354815855386E0), + L(1.4403930189057632173997301031392126865694E0), + L(1.4424730991091018200252920599377292525125E0), + L(1.4444882097316563655148453598508037025938E0), + L(1.4464413322481351841999668424758804165254E0), + L(1.4483352693775551917970437843145232637695E0), + L(1.4501726582147939000905940595923466567576E0), + L(1.4519559822271314199339700039142990228105E0), + L(1.4536875822280323362423034480994649820285E0), + L(1.4553696664279718992423082296859928222270E0), + L(1.4570043196511885530074841089245667532358E0), + L(1.4585935117976422128825857356750737658039E0), + L(1.4601391056210009726721818194296893361233E0), + L(1.4616428638860188872060496086383008594310E0), + L(1.4631064559620759326975975316301202111560E0), + L(1.4645314639038178118428450961503371619177E0), + L(1.4659193880646627234129855241049975398470E0), + L(1.4672716522843522691530527207287398276197E0), + L(1.4685896086876430842559640450619880951144E0), + L(1.4698745421276027686510391411132998919794E0), + L(1.4711276743037345918528755717617308518553E0), + L(1.4723501675822635384916444186631899205983E0), + L(1.4735431285433308455179928682541563973416E0), /* arctan(10.25) */ + L(1.5707963267948966192313216916397514420986E0) /* pi/2 */ +}; + + +/* arctan t = t + t^3 p(t^2) / q(t^2) + |t| <= 0.09375 + peak relative error 5.3e-37 */ + +static const _Float128 + p0 = L(-4.283708356338736809269381409828726405572E1), + p1 = L(-8.636132499244548540964557273544599863825E1), + p2 = L(-5.713554848244551350855604111031839613216E1), + p3 = L(-1.371405711877433266573835355036413750118E1), + p4 = L(-8.638214309119210906997318946650189640184E-1), + q0 = L(1.285112506901621042780814422948906537959E2), + q1 = L(3.361907253914337187957855834229672347089E2), + q2 = L(3.180448303864130128268191635189365331680E2), + q3 = L(1.307244136980865800160844625025280344686E2), + q4 = L(2.173623741810414221251136181221172551416E1); + /* q5 = 1.000000000000000000000000000000000000000E0 */ + +static const _Float128 huge = L(1.0e4930); + +_Float128 +__atanl (_Float128 x) +{ + int k, sign; + _Float128 t, u, p, q; + ieee854_long_double_shape_type s; + + s.value = x; + k = s.parts32.w0; + if (k & 0x80000000) + sign = 1; + else + sign = 0; + + /* Check for IEEE special cases. */ + k &= 0x7fffffff; + if (k >= 0x7fff0000) + { + /* NaN. */ + if ((k & 0xffff) | s.parts32.w1 | s.parts32.w2 | s.parts32.w3) + return (x + x); + + /* Infinity. */ + if (sign) + return -atantbl[83]; + else + return atantbl[83]; + } + + if (k <= 0x3fc50000) /* |x| < 2**-58 */ + { + math_check_force_underflow (x); + /* Raise inexact. */ + if (huge + x > 0.0) + return x; + } + + if (k >= 0x40720000) /* |x| > 2**115 */ + { + /* Saturate result to {-,+}pi/2 */ + if (sign) + return -atantbl[83]; + else + return atantbl[83]; + } + + if (sign) + x = -x; + + if (k >= 0x40024800) /* 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 = L(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 */ + u = atantbl[k] + u; + if (sign) + return (-u); + else + return u; +} + +weak_alias (__atanl, atanl) |