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