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