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