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