sysdeps/ieee754/ldbl-128/s_atanl.c

   1 /*                                                      s_atanl.c
   2  *
   3  *      Inverse circular tangent for 128-bit long double precision
   4  *      (arctangent)
   5  *
   6  *
   7  *
   8  * SYNOPSIS:
   9  *
  10  * long double x, y, atanl();
  11  *
  12  * y = atanl( x );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
  19  *
  20  * The function uses a rational approximation of the form
  21  * t + t^3 P(t^2)/Q(t^2), optimized for |t| < 0.09375.
  22  *
  23  * The argument is reduced using the identity
  24  *    arctan x - arctan u  =  arctan ((x-u)/(1 + ux))
  25  * and an 83-entry lookup table for arctan u, with u = 0, 1/8, ..., 10.25.
  26  * Use of the table improves the execution speed of the routine.
  27  *
  28  *
  29  *
  30  * ACCURACY:
  31  *
  32  *                      Relative error:
  33  * arithmetic   domain     # trials      peak         rms
  34  *    IEEE      -19, 19       4e5       1.7e-34     5.4e-35
  35  *
  36  *
  37  * WARNING:
  38  *
  39  * This program uses integer operations on bit fields of floating-point
  40  * numbers.  It does not work with data structures other than the
  41  * structure assumed.
  42  *
  43  */
  44
  45 /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
  46
  47     This library is free software; you can redistribute it and/or
  48     modify it under the terms of the GNU Lesser General Public
  49     License as published by the Free Software Foundation; either
  50     version 2.1 of the License, or (at your option) any later version.
  51
  52     This library is distributed in the hope that it will be useful,
  53     but WITHOUT ANY WARRANTY; without even the implied warranty of
  54     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  55     Lesser General Public License for more details.
  56
  57     You should have received a copy of the GNU Lesser General Public
  58     License along with this library; if not, see
  59     <http://www.gnu.org/licenses/>.  */
  60
  61
  62 #include <float.h>
  63 #include <math.h>
  64 #include <math_private.h>
  65
  66 /* arctan(k/8), k = 0, ..., 82 */
  67 static const long double atantbl[84] = {
  68   0.0000000000000000000000000000000000000000E0L,
  69   1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125)  */
  70   2.4497866312686415417208248121127581091414E-1L,
  71   3.5877067027057222039592006392646049977698E-1L,
  72   4.6364760900080611621425623146121440202854E-1L,
  73   5.5859931534356243597150821640166127034645E-1L,
  74   6.4350110879328438680280922871732263804151E-1L,
  75   7.1882999962162450541701415152590465395142E-1L,
  76   7.8539816339744830961566084581987572104929E-1L,
  77   8.4415398611317100251784414827164750652594E-1L,
  78   8.9605538457134395617480071802993782702458E-1L,
  79   9.4200004037946366473793717053459358607166E-1L,
  80   9.8279372324732906798571061101466601449688E-1L,
  81   1.0191413442663497346383429170230636487744E0L,
  82   1.0516502125483736674598673120862998296302E0L,
  83   1.0808390005411683108871567292171998202703E0L,
  84   1.1071487177940905030170654601785370400700E0L,
  85   1.1309537439791604464709335155363278047493E0L,
  86   1.1525719972156675180401498626127513797495E0L,
  87   1.1722738811284763866005949441337046149712E0L,
  88   1.1902899496825317329277337748293183376012E0L,
  89   1.2068173702852525303955115800565576303133E0L,
  90   1.2220253232109896370417417439225704908830E0L,
  91   1.2360594894780819419094519711090786987027E0L,
  92   1.2490457723982544258299170772810901230778E0L,
  93   1.2610933822524404193139408812473357720101E0L,
  94   1.2722973952087173412961937498224804940684E0L,
  95   1.2827408797442707473628852511364955306249E0L,
  96   1.2924966677897852679030914214070816845853E0L,
  97   1.3016288340091961438047858503666855921414E0L,
  98   1.3101939350475556342564376891719053122733E0L,
  99   1.3182420510168370498593302023271362531155E0L,
 100   1.3258176636680324650592392104284756311844E0L,
 101   1.3329603993374458675538498697331558093700E0L,
 102   1.3397056595989995393283037525895557411039E0L,
 103   1.3460851583802539310489409282517796256512E0L,
 104   1.3521273809209546571891479413898128509842E0L,
 105   1.3578579772154994751124898859640585287459E0L,
 106   1.3633001003596939542892985278250991189943E0L,
 107   1.3684746984165928776366381936948529556191E0L,
 108   1.3734007669450158608612719264449611486510E0L,
 109   1.3780955681325110444536609641291551522494E0L,
 110   1.3825748214901258580599674177685685125566E0L,
 111   1.3868528702577214543289381097042486034883E0L,
 112   1.3909428270024183486427686943836432060856E0L,
 113   1.3948567013423687823948122092044222644895E0L,
 114   1.3986055122719575950126700816114282335732E0L,
 115   1.4021993871854670105330304794336492676944E0L,
 116   1.4056476493802697809521934019958079881002E0L,
 117   1.4089588955564736949699075250792569287156E0L,
 118   1.4121410646084952153676136718584891599630E0L,
 119   1.4152014988178669079462550975833894394929E0L,
 120   1.4181469983996314594038603039700989523716E0L,
 121   1.4209838702219992566633046424614466661176E0L,
 122   1.4237179714064941189018190466107297503086E0L,
 123   1.4263547484202526397918060597281265695725E0L,
 124   1.4288992721907326964184700745371983590908E0L,
 125   1.4313562697035588982240194668401779312122E0L,
 126   1.4337301524847089866404719096698873648610E0L,
 127   1.4360250423171655234964275337155008780675E0L,
 128   1.4382447944982225979614042479354815855386E0L,
 129   1.4403930189057632173997301031392126865694E0L,
 130   1.4424730991091018200252920599377292525125E0L,
 131   1.4444882097316563655148453598508037025938E0L,
 132   1.4464413322481351841999668424758804165254E0L,
 133   1.4483352693775551917970437843145232637695E0L,
 134   1.4501726582147939000905940595923466567576E0L,
 135   1.4519559822271314199339700039142990228105E0L,
 136   1.4536875822280323362423034480994649820285E0L,
 137   1.4553696664279718992423082296859928222270E0L,
 138   1.4570043196511885530074841089245667532358E0L,
 139   1.4585935117976422128825857356750737658039E0L,
 140   1.4601391056210009726721818194296893361233E0L,
 141   1.4616428638860188872060496086383008594310E0L,
 142   1.4631064559620759326975975316301202111560E0L,
 143   1.4645314639038178118428450961503371619177E0L,
 144   1.4659193880646627234129855241049975398470E0L,
 145   1.4672716522843522691530527207287398276197E0L,
 146   1.4685896086876430842559640450619880951144E0L,
 147   1.4698745421276027686510391411132998919794E0L,
 148   1.4711276743037345918528755717617308518553E0L,
 149   1.4723501675822635384916444186631899205983E0L,
 150   1.4735431285433308455179928682541563973416E0L, /* arctan(10.25) */
 151   1.5707963267948966192313216916397514420986E0L  /* pi/2 */
 152 };
 153
 154
 155 /* arctan t = t + t^3 p(t^2) / q(t^2)
 156    |t| <= 0.09375
 157    peak relative error 5.3e-37 */
 158
 159 static const long double
 160   p0 = -4.283708356338736809269381409828726405572E1L,
 161   p1 = -8.636132499244548540964557273544599863825E1L,
 162   p2 = -5.713554848244551350855604111031839613216E1L,
 163   p3 = -1.371405711877433266573835355036413750118E1L,
 164   p4 = -8.638214309119210906997318946650189640184E-1L,
 165   q0 = 1.285112506901621042780814422948906537959E2L,
 166   q1 = 3.361907253914337187957855834229672347089E2L,
 167   q2 = 3.180448303864130128268191635189365331680E2L,
 168   q3 = 1.307244136980865800160844625025280344686E2L,
 169   q4 = 2.173623741810414221251136181221172551416E1L;
 170   /* q5 = 1.000000000000000000000000000000000000000E0 */
 171
 172 static const long double huge = 1.0e4930L;
 173
 174 long double
 175 __atanl (long double x)
 176 {
 177   int k, sign;
 178   long double t, u, p, q;
 179   ieee854_long_double_shape_type s;
 180
 181   s.value = x;
 182   k = s.parts32.w0;
 183   if (k & 0x80000000)
 184     sign = 1;
 185   else
 186     sign = 0;
 187
 188   /* Check for IEEE special cases.  */
 189   k &= 0x7fffffff;
 190   if (k >= 0x7fff0000)
 191     {
 192       /* NaN. */
 193       if ((k & 0xffff) | s.parts32.w1 | s.parts32.w2 | s.parts32.w3)
 194         return (x + x);
 195
 196       /* Infinity. */
 197       if (sign)
 198         return -atantbl[83];
 199       else
 200         return atantbl[83];
 201     }
 202
 203   if (k <= 0x3fc50000) /* |x| < 2**-58 */
 204     {
 205       if (fabsl (x) < LDBL_MIN)
 206         {
 207           long double force_underflow = x * x;
 208           math_force_eval (force_underflow);
 209         }
 210       /* Raise inexact. */
 211       if (huge + x > 0.0)
 212         return x;
 213     }
 214
 215   if (k >= 0x40720000) /* |x| > 2**115 */
 216     {
 217       /* Saturate result to {-,+}pi/2 */
 218       if (sign)
 219         return -atantbl[83];
 220       else
 221         return atantbl[83];
 222     }
 223
 224   if (sign)
 225       x = -x;
 226
 227   if (k >= 0x40024800) /* 10.25 */
 228     {
 229       k = 83;
 230       t = -1.0/x;
 231     }
 232   else
 233     {
 234       /* Index of nearest table element.
 235          Roundoff to integer is asymmetrical to avoid cancellation when t < 0
 236          (cf. fdlibm). */
 237       k = 8.0 * x + 0.25;
 238       u = 0.125 * k;
 239       /* Small arctan argument.  */
 240       t = (x - u) / (1.0 + x * u);
 241     }
 242
 243   /* Arctan of small argument t.  */
 244   u = t * t;
 245   p =     ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
 246   q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
 247   u = t * u * p / q  +  t;
 248
 249   /* arctan x = arctan u  +  arctan t */
 250   u = atantbl[k] + u;
 251   if (sign)
 252     return (-u);
 253   else
 254     return u;
 255 }
 256
 257 weak_alias (__atanl, atanl)