lib/atanl.c

   1 /*                                                      s_atanl.c
   2  *
   3  *      Inverse circular tangent for 128-bit long double precision
   4  *      (arctangent)
   5  *
   6  *
   7  *
   8  * SYNOPSIS:
   9  *
  10  * long double x, y, atanl();
  11  *
  12  * y = atanl( x );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
  19  *
  20  * The function uses a rational approximation of the form
  21  * t + t^3 P(t^2)/Q(t^2), optimized for |t| < 0.09375.
  22  *
  23  * The argument is reduced using the identity
  24  *    arctan x - arctan u  =  arctan ((x-u)/(1 + ux))
  25  * and an 83-entry lookup table for arctan u, with u = 0, 1/8, ..., 10.25.
  26  * Use of the table improves the execution speed of the routine.
  27  *
  28  *
  29  *
  30  * ACCURACY:
  31  *
  32  *                      Relative error:
  33  * arithmetic   domain     # trials      peak         rms
  34  *    IEEE      -19, 19       4e5       1.7e-34     5.4e-35
  35  *
  36  *
  37  * WARNING:
  38  *
  39  * This program uses integer operations on bit fields of floating-point
  40  * numbers.  It does not work with data structures other than the
  41  * structure assumed.
  42  *
  43  */
  44
  45 /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>  */
  46
  47 #include "mathl.h"
  48
  49 #include <math.h>
  50
  51 /* arctan(k/8), k = 0, ..., 82 */
  52 static const long double atantbl[84] = {
  53   0.0000000000000000000000000000000000000000E0L,
  54   1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125)  */
  55   2.4497866312686415417208248121127581091414E-1L,
  56   3.5877067027057222039592006392646049977698E-1L,
  57   4.6364760900080611621425623146121440202854E-1L,
  58   5.5859931534356243597150821640166127034645E-1L,
  59   6.4350110879328438680280922871732263804151E-1L,
  60   7.1882999962162450541701415152590465395142E-1L,
  61   7.8539816339744830961566084581987572104929E-1L,
  62   8.4415398611317100251784414827164750652594E-1L,
  63   8.9605538457134395617480071802993782702458E-1L,
  64   9.4200004037946366473793717053459358607166E-1L,
  65   9.8279372324732906798571061101466601449688E-1L,
  66   1.0191413442663497346383429170230636487744E0L,
  67   1.0516502125483736674598673120862998296302E0L,
  68   1.0808390005411683108871567292171998202703E0L,
  69   1.1071487177940905030170654601785370400700E0L,
  70   1.1309537439791604464709335155363278047493E0L,
  71   1.1525719972156675180401498626127513797495E0L,
  72   1.1722738811284763866005949441337046149712E0L,
  73   1.1902899496825317329277337748293183376012E0L,
  74   1.2068173702852525303955115800565576303133E0L,
  75   1.2220253232109896370417417439225704908830E0L,
  76   1.2360594894780819419094519711090786987027E0L,
  77   1.2490457723982544258299170772810901230778E0L,
  78   1.2610933822524404193139408812473357720101E0L,
  79   1.2722973952087173412961937498224804940684E0L,
  80   1.2827408797442707473628852511364955306249E0L,
  81   1.2924966677897852679030914214070816845853E0L,
  82   1.3016288340091961438047858503666855921414E0L,
  83   1.3101939350475556342564376891719053122733E0L,
  84   1.3182420510168370498593302023271362531155E0L,
  85   1.3258176636680324650592392104284756311844E0L,
  86   1.3329603993374458675538498697331558093700E0L,
  87   1.3397056595989995393283037525895557411039E0L,
  88   1.3460851583802539310489409282517796256512E0L,
  89   1.3521273809209546571891479413898128509842E0L,
  90   1.3578579772154994751124898859640585287459E0L,
  91   1.3633001003596939542892985278250991189943E0L,
  92   1.3684746984165928776366381936948529556191E0L,
  93   1.3734007669450158608612719264449611486510E0L,
  94   1.3780955681325110444536609641291551522494E0L,
  95   1.3825748214901258580599674177685685125566E0L,
  96   1.3868528702577214543289381097042486034883E0L,
  97   1.3909428270024183486427686943836432060856E0L,
  98   1.3948567013423687823948122092044222644895E0L,
  99   1.3986055122719575950126700816114282335732E0L,
 100   1.4021993871854670105330304794336492676944E0L,
 101   1.4056476493802697809521934019958079881002E0L,
 102   1.4089588955564736949699075250792569287156E0L,
 103   1.4121410646084952153676136718584891599630E0L,
 104   1.4152014988178669079462550975833894394929E0L,
 105   1.4181469983996314594038603039700989523716E0L,
 106   1.4209838702219992566633046424614466661176E0L,
 107   1.4237179714064941189018190466107297503086E0L,
 108   1.4263547484202526397918060597281265695725E0L,
 109   1.4288992721907326964184700745371983590908E0L,
 110   1.4313562697035588982240194668401779312122E0L,
 111   1.4337301524847089866404719096698873648610E0L,
 112   1.4360250423171655234964275337155008780675E0L,
 113   1.4382447944982225979614042479354815855386E0L,
 114   1.4403930189057632173997301031392126865694E0L,
 115   1.4424730991091018200252920599377292525125E0L,
 116   1.4444882097316563655148453598508037025938E0L,
 117   1.4464413322481351841999668424758804165254E0L,
 118   1.4483352693775551917970437843145232637695E0L,
 119   1.4501726582147939000905940595923466567576E0L,
 120   1.4519559822271314199339700039142990228105E0L,
 121   1.4536875822280323362423034480994649820285E0L,
 122   1.4553696664279718992423082296859928222270E0L,
 123   1.4570043196511885530074841089245667532358E0L,
 124   1.4585935117976422128825857356750737658039E0L,
 125   1.4601391056210009726721818194296893361233E0L,
 126   1.4616428638860188872060496086383008594310E0L,
 127   1.4631064559620759326975975316301202111560E0L,
 128   1.4645314639038178118428450961503371619177E0L,
 129   1.4659193880646627234129855241049975398470E0L,
 130   1.4672716522843522691530527207287398276197E0L,
 131   1.4685896086876430842559640450619880951144E0L,
 132   1.4698745421276027686510391411132998919794E0L,
 133   1.4711276743037345918528755717617308518553E0L,
 134   1.4723501675822635384916444186631899205983E0L,
 135   1.4735431285433308455179928682541563973416E0L, /* arctan(10.25) */
 136   1.5707963267948966192313216916397514420986E0L  /* pi/2 */
 137 };
 138
 139
 140 /* arctan t = t + t^3 p(t^2) / q(t^2)
 141    |t| <= 0.09375
 142    peak relative error 5.3e-37 */
 143
 144 static const long double
 145   p0 = -4.283708356338736809269381409828726405572E1L,
 146   p1 = -8.636132499244548540964557273544599863825E1L,
 147   p2 = -5.713554848244551350855604111031839613216E1L,
 148   p3 = -1.371405711877433266573835355036413750118E1L,
 149   p4 = -8.638214309119210906997318946650189640184E-1L,
 150   q0 = 1.285112506901621042780814422948906537959E2L,
 151   q1 = 3.361907253914337187957855834229672347089E2L,
 152   q2 = 3.180448303864130128268191635189365331680E2L,
 153   q3 = 1.307244136980865800160844625025280344686E2L,
 154   q4 = 2.173623741810414221251136181221172551416E1L;
 155   /* q5 = 1.000000000000000000000000000000000000000E0 */
 156
 157
 158 long double
 159 atanl (long double x)
 160 {
 161   int k, sign;
 162   long double t, u, p, q;
 163
 164   sign = (x < 0) ? -1 : 1;
 165
 166   /* Check for zero or NaN.  */
 167   if (x != x || x == 0.0)
 168     return x + x;
 169
 170   /* Check for infinity.  */
 171   if (x + x == x)
 172     return sign * atantbl[83];
 173
 174   x *= sign;
 175
 176   if (x >= 10.25) /* 10.25 */
 177     {
 178       k = 83;
 179       t = -1.0/x;
 180     }
 181   else
 182     {
 183       /* Index of nearest table element.
 184          Roundoff to integer is asymmetrical to avoid cancellation when t < 0
 185          (cf. fdlibm). */
 186       k = 8.0 * x + 0.25;
 187       u = 0.125 * k;
 188       /* Small arctan argument.  */
 189       t = (x - u) / (1.0 + x * u);
 190     }
 191
 192   /* Arctan of small argument t.  */
 193   u = t * t;
 194   p =     ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
 195   q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
 196   u = t * u * p / q  +  t;
 197
 198   /* arctan x = arctan u  +  arctan t */
 199   return sign * (atantbl[k] + u);
 200 }