lib/atanl.c

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