libquadmath/math/atanq.c

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