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    This program is free software; you can redistribute it and/or modify
  48    it under the terms of the GNU General Public License as published by
  49    the Free Software Foundation; either version 2, or (at your option)
  50    any later version.
  51
  52    This program is distributed in the hope that it will be useful,
  53    but WITHOUT ANY WARRANTY; without even the implied warranty of
  54    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  55    GNU General Public License for more details.
  56
  57    You should have received a copy of the GNU General Public License
  58    along with this program; see the file COPYING.
  59    If not, write to the Free Software Foundation,
  60    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.  */
  61
  62 #include "mathl.h"
  63
  64 #include <math.h>
  65
  66 /* arctan(k/8), k = 0, ..., 82 */
  67 static const long double atantbl[84] = {
  68   0.0000000000000000000000000000000000000000E0L,
  69   1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125)  */
  70   2.4497866312686415417208248121127581091414E-1L,
  71   3.5877067027057222039592006392646049977698E-1L,
  72   4.6364760900080611621425623146121440202854E-1L,
  73   5.5859931534356243597150821640166127034645E-1L,
  74   6.4350110879328438680280922871732263804151E-1L,
  75   7.1882999962162450541701415152590465395142E-1L,
  76   7.8539816339744830961566084581987572104929E-1L,
  77   8.4415398611317100251784414827164750652594E-1L,
  78   8.9605538457134395617480071802993782702458E-1L,
  79   9.4200004037946366473793717053459358607166E-1L,
  80   9.8279372324732906798571061101466601449688E-1L,
  81   1.0191413442663497346383429170230636487744E0L,
  82   1.0516502125483736674598673120862998296302E0L,
  83   1.0808390005411683108871567292171998202703E0L,
  84   1.1071487177940905030170654601785370400700E0L,
  85   1.1309537439791604464709335155363278047493E0L,
  86   1.1525719972156675180401498626127513797495E0L,
  87   1.1722738811284763866005949441337046149712E0L,
  88   1.1902899496825317329277337748293183376012E0L,
  89   1.2068173702852525303955115800565576303133E0L,
  90   1.2220253232109896370417417439225704908830E0L,
  91   1.2360594894780819419094519711090786987027E0L,
  92   1.2490457723982544258299170772810901230778E0L,
  93   1.2610933822524404193139408812473357720101E0L,
  94   1.2722973952087173412961937498224804940684E0L,
  95   1.2827408797442707473628852511364955306249E0L,
  96   1.2924966677897852679030914214070816845853E0L,
  97   1.3016288340091961438047858503666855921414E0L,
  98   1.3101939350475556342564376891719053122733E0L,
  99   1.3182420510168370498593302023271362531155E0L,
 100   1.3258176636680324650592392104284756311844E0L,
 101   1.3329603993374458675538498697331558093700E0L,
 102   1.3397056595989995393283037525895557411039E0L,
 103   1.3460851583802539310489409282517796256512E0L,
 104   1.3521273809209546571891479413898128509842E0L,
 105   1.3578579772154994751124898859640585287459E0L,
 106   1.3633001003596939542892985278250991189943E0L,
 107   1.3684746984165928776366381936948529556191E0L,
 108   1.3734007669450158608612719264449611486510E0L,
 109   1.3780955681325110444536609641291551522494E0L,
 110   1.3825748214901258580599674177685685125566E0L,
 111   1.3868528702577214543289381097042486034883E0L,
 112   1.3909428270024183486427686943836432060856E0L,
 113   1.3948567013423687823948122092044222644895E0L,
 114   1.3986055122719575950126700816114282335732E0L,
 115   1.4021993871854670105330304794336492676944E0L,
 116   1.4056476493802697809521934019958079881002E0L,
 117   1.4089588955564736949699075250792569287156E0L,
 118   1.4121410646084952153676136718584891599630E0L,
 119   1.4152014988178669079462550975833894394929E0L,
 120   1.4181469983996314594038603039700989523716E0L,
 121   1.4209838702219992566633046424614466661176E0L,
 122   1.4237179714064941189018190466107297503086E0L,
 123   1.4263547484202526397918060597281265695725E0L,
 124   1.4288992721907326964184700745371983590908E0L,
 125   1.4313562697035588982240194668401779312122E0L,
 126   1.4337301524847089866404719096698873648610E0L,
 127   1.4360250423171655234964275337155008780675E0L,
 128   1.4382447944982225979614042479354815855386E0L,
 129   1.4403930189057632173997301031392126865694E0L,
 130   1.4424730991091018200252920599377292525125E0L,
 131   1.4444882097316563655148453598508037025938E0L,
 132   1.4464413322481351841999668424758804165254E0L,
 133   1.4483352693775551917970437843145232637695E0L,
 134   1.4501726582147939000905940595923466567576E0L,
 135   1.4519559822271314199339700039142990228105E0L,
 136   1.4536875822280323362423034480994649820285E0L,
 137   1.4553696664279718992423082296859928222270E0L,
 138   1.4570043196511885530074841089245667532358E0L,
 139   1.4585935117976422128825857356750737658039E0L,
 140   1.4601391056210009726721818194296893361233E0L,
 141   1.4616428638860188872060496086383008594310E0L,
 142   1.4631064559620759326975975316301202111560E0L,
 143   1.4645314639038178118428450961503371619177E0L,
 144   1.4659193880646627234129855241049975398470E0L,
 145   1.4672716522843522691530527207287398276197E0L,
 146   1.4685896086876430842559640450619880951144E0L,
 147   1.4698745421276027686510391411132998919794E0L,
 148   1.4711276743037345918528755717617308518553E0L,
 149   1.4723501675822635384916444186631899205983E0L,
 150   1.4735431285433308455179928682541563973416E0L, /* arctan(10.25) */
 151   1.5707963267948966192313216916397514420986E0L  /* pi/2 */
 152 };
 153
 154
 155 /* arctan t = t + t^3 p(t^2) / q(t^2)
 156    |t| <= 0.09375
 157    peak relative error 5.3e-37 */
 158
 159 static const long double
 160   p0 = -4.283708356338736809269381409828726405572E1L,
 161   p1 = -8.636132499244548540964557273544599863825E1L,
 162   p2 = -5.713554848244551350855604111031839613216E1L,
 163   p3 = -1.371405711877433266573835355036413750118E1L,
 164   p4 = -8.638214309119210906997318946650189640184E-1L,
 165   q0 = 1.285112506901621042780814422948906537959E2L,
 166   q1 = 3.361907253914337187957855834229672347089E2L,
 167   q2 = 3.180448303864130128268191635189365331680E2L,
 168   q3 = 1.307244136980865800160844625025280344686E2L,
 169   q4 = 2.173623741810414221251136181221172551416E1L;
 170   /* q5 = 1.000000000000000000000000000000000000000E0 */
 171
 172
 173 long double
 174 atanl (long double x)
 175 {
 176   int k, sign;
 177   long double t, u, p, q;
 178
 179   sign = x < 0.0;
 180
 181   /* Check for zero or NaN.  */
 182   if (x != x || x == 0.0)
 183     return x + x;
 184
 185   if (x + x == x)
 186     {
 187       /* Infinity. */
 188       if (sign)
 189         return -atantbl[83];
 190       else
 191         return atantbl[83];
 192     }
 193
 194   if (sign)
 195       x = -x;
 196
 197   if (x >= 10.25)
 198     {
 199       k = 83;
 200       t = -1.0/x;
 201     }
 202   else
 203     {
 204       /* Index of nearest table element.
 205          Roundoff to integer is asymmetrical to avoid cancellation when t < 0
 206          (cf. fdlibm). */
 207       k = 8.0 * x + 0.25;
 208       u = 0.125 * k;
 209       /* Small arctan argument.  */
 210       t = (x - u) / (1.0 + x * u);
 211     }
 212
 213   /* Arctan of small argument t.  */
 214   u = t * t;
 215   p =     ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
 216   q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
 217   u = t * u * p / q  +  t;
 218
 219   /* arctan x = arctan u  +  arctan t */
 220   u = atantbl[k] + u;
 221   if (sign)
 222     return (-u);
 223   else
 224     return u;
 225 }