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