Back to index

glibc  2.9
s_atanl.c
Go to the documentation of this file.
00001 /*                                               s_atanl.c
00002  *
00003  *     Inverse circular tangent for 128-bit long double precision
00004  *      (arctangent)
00005  *
00006  *
00007  *
00008  * SYNOPSIS:
00009  *
00010  * long double x, y, atanl();
00011  *
00012  * y = atanl( x );
00013  *
00014  *
00015  *
00016  * DESCRIPTION:
00017  *
00018  * Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
00019  *
00020  * The function uses a rational approximation of the form
00021  * t + t^3 P(t^2)/Q(t^2), optimized for |t| < 0.09375.
00022  *
00023  * The argument is reduced using the identity
00024  *    arctan x - arctan u  =  arctan ((x-u)/(1 + ux))
00025  * and an 83-entry lookup table for arctan u, with u = 0, 1/8, ..., 10.25.
00026  * Use of the table improves the execution speed of the routine.
00027  *
00028  *
00029  *
00030  * ACCURACY:
00031  *
00032  *                      Relative error:
00033  * arithmetic   domain     # trials      peak         rms
00034  *    IEEE      -19, 19       4e5       1.7e-34     5.4e-35
00035  *
00036  *
00037  * WARNING:
00038  *
00039  * This program uses integer operations on bit fields of floating-point
00040  * numbers.  It does not work with data structures other than the
00041  * structure assumed.
00042  *
00043  */
00044 
00045 /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov> 
00046 
00047     This library is free software; you can redistribute it and/or
00048     modify it under the terms of the GNU Lesser General Public
00049     License as published by the Free Software Foundation; either
00050     version 2.1 of the License, or (at your option) any later version.
00051 
00052     This library is distributed in the hope that it will be useful,
00053     but WITHOUT ANY WARRANTY; without even the implied warranty of
00054     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00055     Lesser General Public License for more details.
00056 
00057     You should have received a copy of the GNU Lesser General Public
00058     License along with this library; if not, write to the Free Software
00059     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00060 
00061 
00062 #include "math_private.h"
00063 
00064 /* arctan(k/8), k = 0, ..., 82 */
00065 static const long double atantbl[84] = {
00066   0.0000000000000000000000000000000000000000E0L,
00067   1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125)  */
00068   2.4497866312686415417208248121127581091414E-1L,
00069   3.5877067027057222039592006392646049977698E-1L,
00070   4.6364760900080611621425623146121440202854E-1L,
00071   5.5859931534356243597150821640166127034645E-1L,
00072   6.4350110879328438680280922871732263804151E-1L,
00073   7.1882999962162450541701415152590465395142E-1L,
00074   7.8539816339744830961566084581987572104929E-1L,
00075   8.4415398611317100251784414827164750652594E-1L,
00076   8.9605538457134395617480071802993782702458E-1L,
00077   9.4200004037946366473793717053459358607166E-1L,
00078   9.8279372324732906798571061101466601449688E-1L,
00079   1.0191413442663497346383429170230636487744E0L,
00080   1.0516502125483736674598673120862998296302E0L,
00081   1.0808390005411683108871567292171998202703E0L,
00082   1.1071487177940905030170654601785370400700E0L,
00083   1.1309537439791604464709335155363278047493E0L,
00084   1.1525719972156675180401498626127513797495E0L,
00085   1.1722738811284763866005949441337046149712E0L,
00086   1.1902899496825317329277337748293183376012E0L,
00087   1.2068173702852525303955115800565576303133E0L,
00088   1.2220253232109896370417417439225704908830E0L,
00089   1.2360594894780819419094519711090786987027E0L,
00090   1.2490457723982544258299170772810901230778E0L,
00091   1.2610933822524404193139408812473357720101E0L,
00092   1.2722973952087173412961937498224804940684E0L,
00093   1.2827408797442707473628852511364955306249E0L,
00094   1.2924966677897852679030914214070816845853E0L,
00095   1.3016288340091961438047858503666855921414E0L,
00096   1.3101939350475556342564376891719053122733E0L,
00097   1.3182420510168370498593302023271362531155E0L,
00098   1.3258176636680324650592392104284756311844E0L,
00099   1.3329603993374458675538498697331558093700E0L,
00100   1.3397056595989995393283037525895557411039E0L,
00101   1.3460851583802539310489409282517796256512E0L,
00102   1.3521273809209546571891479413898128509842E0L,
00103   1.3578579772154994751124898859640585287459E0L,
00104   1.3633001003596939542892985278250991189943E0L,
00105   1.3684746984165928776366381936948529556191E0L,
00106   1.3734007669450158608612719264449611486510E0L,
00107   1.3780955681325110444536609641291551522494E0L,
00108   1.3825748214901258580599674177685685125566E0L,
00109   1.3868528702577214543289381097042486034883E0L,
00110   1.3909428270024183486427686943836432060856E0L,
00111   1.3948567013423687823948122092044222644895E0L,
00112   1.3986055122719575950126700816114282335732E0L,
00113   1.4021993871854670105330304794336492676944E0L,
00114   1.4056476493802697809521934019958079881002E0L,
00115   1.4089588955564736949699075250792569287156E0L,
00116   1.4121410646084952153676136718584891599630E0L,
00117   1.4152014988178669079462550975833894394929E0L,
00118   1.4181469983996314594038603039700989523716E0L,
00119   1.4209838702219992566633046424614466661176E0L,
00120   1.4237179714064941189018190466107297503086E0L,
00121   1.4263547484202526397918060597281265695725E0L,
00122   1.4288992721907326964184700745371983590908E0L,
00123   1.4313562697035588982240194668401779312122E0L,
00124   1.4337301524847089866404719096698873648610E0L,
00125   1.4360250423171655234964275337155008780675E0L,
00126   1.4382447944982225979614042479354815855386E0L,
00127   1.4403930189057632173997301031392126865694E0L,
00128   1.4424730991091018200252920599377292525125E0L,
00129   1.4444882097316563655148453598508037025938E0L,
00130   1.4464413322481351841999668424758804165254E0L,
00131   1.4483352693775551917970437843145232637695E0L,
00132   1.4501726582147939000905940595923466567576E0L,
00133   1.4519559822271314199339700039142990228105E0L,
00134   1.4536875822280323362423034480994649820285E0L,
00135   1.4553696664279718992423082296859928222270E0L,
00136   1.4570043196511885530074841089245667532358E0L,
00137   1.4585935117976422128825857356750737658039E0L,
00138   1.4601391056210009726721818194296893361233E0L,
00139   1.4616428638860188872060496086383008594310E0L,
00140   1.4631064559620759326975975316301202111560E0L,
00141   1.4645314639038178118428450961503371619177E0L,
00142   1.4659193880646627234129855241049975398470E0L,
00143   1.4672716522843522691530527207287398276197E0L,
00144   1.4685896086876430842559640450619880951144E0L,
00145   1.4698745421276027686510391411132998919794E0L,
00146   1.4711276743037345918528755717617308518553E0L,
00147   1.4723501675822635384916444186631899205983E0L,
00148   1.4735431285433308455179928682541563973416E0L, /* arctan(10.25) */
00149   1.5707963267948966192313216916397514420986E0L  /* pi/2 */
00150 };
00151 
00152 
00153 /* arctan t = t + t^3 p(t^2) / q(t^2)
00154    |t| <= 0.09375
00155    peak relative error 5.3e-37 */
00156 
00157 static const long double
00158   p0 = -4.283708356338736809269381409828726405572E1L,
00159   p1 = -8.636132499244548540964557273544599863825E1L,
00160   p2 = -5.713554848244551350855604111031839613216E1L,
00161   p3 = -1.371405711877433266573835355036413750118E1L,
00162   p4 = -8.638214309119210906997318946650189640184E-1L,
00163   q0 = 1.285112506901621042780814422948906537959E2L,
00164   q1 = 3.361907253914337187957855834229672347089E2L,
00165   q2 = 3.180448303864130128268191635189365331680E2L,
00166   q3 = 1.307244136980865800160844625025280344686E2L,
00167   q4 = 2.173623741810414221251136181221172551416E1L;
00168   /* q5 = 1.000000000000000000000000000000000000000E0 */
00169 
00170 
00171 long double
00172 __atanl (long double x)
00173 {
00174   int k, sign;
00175   long double t, u, p, q;
00176   ieee854_long_double_shape_type s;
00177 
00178   s.value = x;
00179   k = s.parts32.w0;
00180   if (k & 0x80000000)
00181     sign = 1;
00182   else
00183     sign = 0;
00184 
00185   /* Check for IEEE special cases.  */
00186   k &= 0x7fffffff;
00187   if (k >= 0x7fff0000)
00188     {
00189       /* NaN. */
00190       if ((k & 0xffff) | s.parts32.w1 | s.parts32.w2 | s.parts32.w3)
00191        return (x + x);
00192 
00193       /* Infinity. */
00194       if (sign)
00195        return -atantbl[83];
00196       else
00197        return atantbl[83];
00198     }
00199 
00200   if (sign)
00201       x = -x;
00202 
00203   if (k >= 0x40024800) /* 10.25 */
00204     {
00205       k = 83;
00206       t = -1.0/x;
00207     }
00208   else
00209     {
00210       /* Index of nearest table element.
00211         Roundoff to integer is asymmetrical to avoid cancellation when t < 0
00212          (cf. fdlibm). */
00213       k = 8.0 * x + 0.25;
00214       u = 0.125 * k;
00215       /* Small arctan argument.  */
00216       t = (x - u) / (1.0 + x * u);
00217     }
00218 
00219   /* Arctan of small argument t.  */
00220   u = t * t;
00221   p =     ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
00222   q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
00223   u = t * u * p / q  +  t;
00224 
00225   /* arctan x = arctan u  +  arctan t */
00226   u = atantbl[k] + u;
00227   if (sign)
00228     return (-u);
00229   else
00230     return u;
00231 }
00232 
00233 weak_alias (__atanl, atanl)