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     <https://www.gnu.org/licenses/>.  */
  60
  61
  62 #include <float.h>
  63 #include <math.h>
  64 #include <math_private.h>
  65 #include <math-underflow.h>
  66 #include <libm-alias-ldouble.h>
  67
  68 /* arctan(k/8), k = 0, ..., 82 */
  69 static const _Float128 atantbl[84] = {
  70   L(0.0000000000000000000000000000000000000000E0),
  71   L(1.2435499454676143503135484916387102557317E-1), /* arctan(0.125)  */
  72   L(2.4497866312686415417208248121127581091414E-1),
  73   L(3.5877067027057222039592006392646049977698E-1),
  74   L(4.6364760900080611621425623146121440202854E-1),
  75   L(5.5859931534356243597150821640166127034645E-1),
  76   L(6.4350110879328438680280922871732263804151E-1),
  77   L(7.1882999962162450541701415152590465395142E-1),
  78   L(7.8539816339744830961566084581987572104929E-1),
  79   L(8.4415398611317100251784414827164750652594E-1),
  80   L(8.9605538457134395617480071802993782702458E-1),
  81   L(9.4200004037946366473793717053459358607166E-1),
  82   L(9.8279372324732906798571061101466601449688E-1),
  83   L(1.0191413442663497346383429170230636487744E0),
  84   L(1.0516502125483736674598673120862998296302E0),
  85   L(1.0808390005411683108871567292171998202703E0),
  86   L(1.1071487177940905030170654601785370400700E0),
  87   L(1.1309537439791604464709335155363278047493E0),
  88   L(1.1525719972156675180401498626127513797495E0),
  89   L(1.1722738811284763866005949441337046149712E0),
  90   L(1.1902899496825317329277337748293183376012E0),
  91   L(1.2068173702852525303955115800565576303133E0),
  92   L(1.2220253232109896370417417439225704908830E0),
  93   L(1.2360594894780819419094519711090786987027E0),
  94   L(1.2490457723982544258299170772810901230778E0),
  95   L(1.2610933822524404193139408812473357720101E0),
  96   L(1.2722973952087173412961937498224804940684E0),
  97   L(1.2827408797442707473628852511364955306249E0),
  98   L(1.2924966677897852679030914214070816845853E0),
  99   L(1.3016288340091961438047858503666855921414E0),
 100   L(1.3101939350475556342564376891719053122733E0),
 101   L(1.3182420510168370498593302023271362531155E0),
 102   L(1.3258176636680324650592392104284756311844E0),
 103   L(1.3329603993374458675538498697331558093700E0),
 104   L(1.3397056595989995393283037525895557411039E0),
 105   L(1.3460851583802539310489409282517796256512E0),
 106   L(1.3521273809209546571891479413898128509842E0),
 107   L(1.3578579772154994751124898859640585287459E0),
 108   L(1.3633001003596939542892985278250991189943E0),
 109   L(1.3684746984165928776366381936948529556191E0),
 110   L(1.3734007669450158608612719264449611486510E0),
 111   L(1.3780955681325110444536609641291551522494E0),
 112   L(1.3825748214901258580599674177685685125566E0),
 113   L(1.3868528702577214543289381097042486034883E0),
 114   L(1.3909428270024183486427686943836432060856E0),
 115   L(1.3948567013423687823948122092044222644895E0),
 116   L(1.3986055122719575950126700816114282335732E0),
 117   L(1.4021993871854670105330304794336492676944E0),
 118   L(1.4056476493802697809521934019958079881002E0),
 119   L(1.4089588955564736949699075250792569287156E0),
 120   L(1.4121410646084952153676136718584891599630E0),
 121   L(1.4152014988178669079462550975833894394929E0),
 122   L(1.4181469983996314594038603039700989523716E0),
 123   L(1.4209838702219992566633046424614466661176E0),
 124   L(1.4237179714064941189018190466107297503086E0),
 125   L(1.4263547484202526397918060597281265695725E0),
 126   L(1.4288992721907326964184700745371983590908E0),
 127   L(1.4313562697035588982240194668401779312122E0),
 128   L(1.4337301524847089866404719096698873648610E0),
 129   L(1.4360250423171655234964275337155008780675E0),
 130   L(1.4382447944982225979614042479354815855386E0),
 131   L(1.4403930189057632173997301031392126865694E0),
 132   L(1.4424730991091018200252920599377292525125E0),
 133   L(1.4444882097316563655148453598508037025938E0),
 134   L(1.4464413322481351841999668424758804165254E0),
 135   L(1.4483352693775551917970437843145232637695E0),
 136   L(1.4501726582147939000905940595923466567576E0),
 137   L(1.4519559822271314199339700039142990228105E0),
 138   L(1.4536875822280323362423034480994649820285E0),
 139   L(1.4553696664279718992423082296859928222270E0),
 140   L(1.4570043196511885530074841089245667532358E0),
 141   L(1.4585935117976422128825857356750737658039E0),
 142   L(1.4601391056210009726721818194296893361233E0),
 143   L(1.4616428638860188872060496086383008594310E0),
 144   L(1.4631064559620759326975975316301202111560E0),
 145   L(1.4645314639038178118428450961503371619177E0),
 146   L(1.4659193880646627234129855241049975398470E0),
 147   L(1.4672716522843522691530527207287398276197E0),
 148   L(1.4685896086876430842559640450619880951144E0),
 149   L(1.4698745421276027686510391411132998919794E0),
 150   L(1.4711276743037345918528755717617308518553E0),
 151   L(1.4723501675822635384916444186631899205983E0),
 152   L(1.4735431285433308455179928682541563973416E0), /* arctan(10.25) */
 153   L(1.5707963267948966192313216916397514420986E0)  /* pi/2 */
 154 };
 155
 156
 157 /* arctan t = t + t^3 p(t^2) / q(t^2)
 158    |t| <= 0.09375
 159    peak relative error 5.3e-37 */
 160
 161 static const _Float128
 162   p0 = L(-4.283708356338736809269381409828726405572E1),
 163   p1 = L(-8.636132499244548540964557273544599863825E1),
 164   p2 = L(-5.713554848244551350855604111031839613216E1),
 165   p3 = L(-1.371405711877433266573835355036413750118E1),
 166   p4 = L(-8.638214309119210906997318946650189640184E-1),
 167   q0 = L(1.285112506901621042780814422948906537959E2),
 168   q1 = L(3.361907253914337187957855834229672347089E2),
 169   q2 = L(3.180448303864130128268191635189365331680E2),
 170   q3 = L(1.307244136980865800160844625025280344686E2),
 171   q4 = L(2.173623741810414221251136181221172551416E1);
 172   /* q5 = 1.000000000000000000000000000000000000000E0 */
 173
 174 static const _Float128 huge = L(1.0e4930);
 175
 176 _Float128
 177 __atanl (_Float128 x)
 178 {
 179   int k, sign;
 180   _Float128 t, u, p, q;
 181   ieee854_long_double_shape_type s;
 182
 183   s.value = x;
 184   k = s.parts32.w0;
 185   if (k & 0x80000000)
 186     sign = 1;
 187   else
 188     sign = 0;
 189
 190   /* Check for IEEE special cases.  */
 191   k &= 0x7fffffff;
 192   if (k >= 0x7fff0000)
 193     {
 194       /* NaN. */
 195       if ((k & 0xffff) | s.parts32.w1 | s.parts32.w2 | s.parts32.w3)
 196         return (x + x);
 197
 198       /* Infinity. */
 199       if (sign)
 200         return -atantbl[83];
 201       else
 202         return atantbl[83];
 203     }
 204
 205   if (k <= 0x3fc50000) /* |x| < 2**-58 */
 206     {
 207       math_check_force_underflow (x);
 208       /* Raise inexact. */
 209       if (huge + x > 0.0)
 210         return x;
 211     }
 212
 213   if (k >= 0x40720000) /* |x| > 2**115 */
 214     {
 215       /* Saturate result to {-,+}pi/2 */
 216       if (sign)
 217         return -atantbl[83];
 218       else
 219         return atantbl[83];
 220     }
 221
 222   if (sign)
 223       x = -x;
 224
 225   if (k >= 0x40024800) /* 10.25 */
 226     {
 227       k = 83;
 228       t = -1.0/x;
 229     }
 230   else
 231     {
 232       /* Index of nearest table element.
 233          Roundoff to integer is asymmetrical to avoid cancellation when t < 0
 234          (cf. fdlibm). */
 235       k = 8.0 * x + 0.25;
 236       u = L(0.125) * k;
 237       /* Small arctan argument.  */
 238       t = (x - u) / (1.0 + x * u);
 239     }
 240
 241   /* Arctan of small argument t.  */
 242   u = t * t;
 243   p =     ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
 244   q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
 245   u = t * u * p / q  +  t;
 246
 247   /* arctan x = arctan u  +  arctan t */
 248   u = atantbl[k] + u;
 249   if (sign)
 250     return (-u);
 251   else
 252     return u;
 253 }
 254
 255 libm_alias_ldouble (__atan, atan)