maint/gnulib/lib/atanl.c

/* */
This source file includes following definitions.
atanl
atanl
   1 /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
   2 
   3    This file is free software: you can redistribute it and/or modify
   4    it under the terms of the GNU Lesser General Public License as
   5    published by the Free Software Foundation; either version 3 of the
   6    License, or (at your option) any later version.
   7 
   8    This file 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 Lesser General Public License for more details.
  12 
  13    You should have received a copy of the GNU Lesser General Public License
  14    along with this program.  If not, see <https://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
/* */
root/maint/gnulib/lib/atanl.c

DEFINITIONS
