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 #include "isnanl.h"
  66
  67 /* arctan(k/8), k = 0, ..., 82 */
  68 static const long double atantbl[84] = {
  69   0.0000000000000000000000000000000000000000E0L,
  70   1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125)  */
  71   2.4497866312686415417208248121127581091414E-1L,
  72   3.5877067027057222039592006392646049977698E-1L,
  73   4.6364760900080611621425623146121440202854E-1L,
  74   5.5859931534356243597150821640166127034645E-1L,
  75   6.4350110879328438680280922871732263804151E-1L,
  76   7.1882999962162450541701415152590465395142E-1L,
  77   7.8539816339744830961566084581987572104929E-1L,
  78   8.4415398611317100251784414827164750652594E-1L,
  79   8.9605538457134395617480071802993782702458E-1L,
  80   9.4200004037946366473793717053459358607166E-1L,
  81   9.8279372324732906798571061101466601449688E-1L,
  82   1.0191413442663497346383429170230636487744E0L,
  83   1.0516502125483736674598673120862998296302E0L,
  84   1.0808390005411683108871567292171998202703E0L,
  85   1.1071487177940905030170654601785370400700E0L,
  86   1.1309537439791604464709335155363278047493E0L,
  87   1.1525719972156675180401498626127513797495E0L,
  88   1.1722738811284763866005949441337046149712E0L,
  89   1.1902899496825317329277337748293183376012E0L,
  90   1.2068173702852525303955115800565576303133E0L,
  91   1.2220253232109896370417417439225704908830E0L,
  92   1.2360594894780819419094519711090786987027E0L,
  93   1.2490457723982544258299170772810901230778E0L,
  94   1.2610933822524404193139408812473357720101E0L,
  95   1.2722973952087173412961937498224804940684E0L,
  96   1.2827408797442707473628852511364955306249E0L,
  97   1.2924966677897852679030914214070816845853E0L,
  98   1.3016288340091961438047858503666855921414E0L,
  99   1.3101939350475556342564376891719053122733E0L,
 100   1.3182420510168370498593302023271362531155E0L,
 101   1.3258176636680324650592392104284756311844E0L,
 102   1.3329603993374458675538498697331558093700E0L,
 103   1.3397056595989995393283037525895557411039E0L,
 104   1.3460851583802539310489409282517796256512E0L,
 105   1.3521273809209546571891479413898128509842E0L,
 106   1.3578579772154994751124898859640585287459E0L,
 107   1.3633001003596939542892985278250991189943E0L,
 108   1.3684746984165928776366381936948529556191E0L,
 109   1.3734007669450158608612719264449611486510E0L,
 110   1.3780955681325110444536609641291551522494E0L,
 111   1.3825748214901258580599674177685685125566E0L,
 112   1.3868528702577214543289381097042486034883E0L,
 113   1.3909428270024183486427686943836432060856E0L,
 114   1.3948567013423687823948122092044222644895E0L,
 115   1.3986055122719575950126700816114282335732E0L,
 116   1.4021993871854670105330304794336492676944E0L,
 117   1.4056476493802697809521934019958079881002E0L,
 118   1.4089588955564736949699075250792569287156E0L,
 119   1.4121410646084952153676136718584891599630E0L,
 120   1.4152014988178669079462550975833894394929E0L,
 121   1.4181469983996314594038603039700989523716E0L,
 122   1.4209838702219992566633046424614466661176E0L,
 123   1.4237179714064941189018190466107297503086E0L,
 124   1.4263547484202526397918060597281265695725E0L,
 125   1.4288992721907326964184700745371983590908E0L,
 126   1.4313562697035588982240194668401779312122E0L,
 127   1.4337301524847089866404719096698873648610E0L,
 128   1.4360250423171655234964275337155008780675E0L,
 129   1.4382447944982225979614042479354815855386E0L,
 130   1.4403930189057632173997301031392126865694E0L,
 131   1.4424730991091018200252920599377292525125E0L,
 132   1.4444882097316563655148453598508037025938E0L,
 133   1.4464413322481351841999668424758804165254E0L,
 134   1.4483352693775551917970437843145232637695E0L,
 135   1.4501726582147939000905940595923466567576E0L,
 136   1.4519559822271314199339700039142990228105E0L,
 137   1.4536875822280323362423034480994649820285E0L,
 138   1.4553696664279718992423082296859928222270E0L,
 139   1.4570043196511885530074841089245667532358E0L,
 140   1.4585935117976422128825857356750737658039E0L,
 141   1.4601391056210009726721818194296893361233E0L,
 142   1.4616428638860188872060496086383008594310E0L,
 143   1.4631064559620759326975975316301202111560E0L,
 144   1.4645314639038178118428450961503371619177E0L,
 145   1.4659193880646627234129855241049975398470E0L,
 146   1.4672716522843522691530527207287398276197E0L,
 147   1.4685896086876430842559640450619880951144E0L,
 148   1.4698745421276027686510391411132998919794E0L,
 149   1.4711276743037345918528755717617308518553E0L,
 150   1.4723501675822635384916444186631899205983E0L,
 151   1.4735431285433308455179928682541563973416E0L, /* arctan(10.25) */
 152   1.5707963267948966192313216916397514420986E0L  /* pi/2 */
 153 };
 154
 155
 156 /* arctan t = t + t^3 p(t^2) / q(t^2)
 157    |t| <= 0.09375
 158    peak relative error 5.3e-37 */
 159
 160 static const long double
 161   p0 = -4.283708356338736809269381409828726405572E1L,
 162   p1 = -8.636132499244548540964557273544599863825E1L,
 163   p2 = -5.713554848244551350855604111031839613216E1L,
 164   p3 = -1.371405711877433266573835355036413750118E1L,
 165   p4 = -8.638214309119210906997318946650189640184E-1L,
 166   q0 = 1.285112506901621042780814422948906537959E2L,
 167   q1 = 3.361907253914337187957855834229672347089E2L,
 168   q2 = 3.180448303864130128268191635189365331680E2L,
 169   q3 = 1.307244136980865800160844625025280344686E2L,
 170   q4 = 2.173623741810414221251136181221172551416E1L;
 171   /* q5 = 1.000000000000000000000000000000000000000E0 */
 172
 173
 174 long double
 175 atanl (long double x)
 176 {
 177   int k, sign;
 178   long double t, u, p, q;
 179
 180   /* Check for zero or NaN.  */
 181   if (isnanl (x) || x == 0.0)
 182     return x + x;
 183
 184   sign = x < 0.0;
 185
 186   if (x + x == x)
 187     {
 188       /* Infinity. */
 189       if (sign)
 190         return -atantbl[83];
 191       else
 192         return atantbl[83];
 193     }
 194
 195   if (sign)
 196       x = -x;
 197
 198   if (x >= 10.25)
 199     {
 200       k = 83;
 201       t = -1.0/x;
 202     }
 203   else
 204     {
 205       /* Index of nearest table element.
 206          Roundoff to integer is asymmetrical to avoid cancellation when t < 0
 207          (cf. fdlibm). */
 208       k = 8.0 * x + 0.25;
 209       u = 0.125 * k;
 210       /* Small arctan argument.  */
 211       t = (x - u) / (1.0 + x * u);
 212     }
 213
 214   /* Arctan of small argument t.  */
 215   u = t * t;
 216   p =     ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
 217   q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
 218   u = t * u * p / q  +  t;
 219
 220   /* arctan x = arctan u  +  arctan t */
 221   u = atantbl[k] + u;
 222   if (sign)
 223     return (-u);
 224   else
 225     return u;
 226 }