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 #include <math_ldbl_opt.h>
00064 
00065 /* arctan(k/8), k = 0, ..., 82 */
00066 static const long double atantbl[84] = {
00067   0.0000000000000000000000000000000000000000E0L,
00068   1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125)  */
00069   2.4497866312686415417208248121127581091414E-1L,
00070   3.5877067027057222039592006392646049977698E-1L,
00071   4.6364760900080611621425623146121440202854E-1L,
00072   5.5859931534356243597150821640166127034645E-1L,
00073   6.4350110879328438680280922871732263804151E-1L,
00074   7.1882999962162450541701415152590465395142E-1L,
00075   7.8539816339744830961566084581987572104929E-1L,
00076   8.4415398611317100251784414827164750652594E-1L,
00077   8.9605538457134395617480071802993782702458E-1L,
00078   9.4200004037946366473793717053459358607166E-1L,
00079   9.8279372324732906798571061101466601449688E-1L,
00080   1.0191413442663497346383429170230636487744E0L,
00081   1.0516502125483736674598673120862998296302E0L,
00082   1.0808390005411683108871567292171998202703E0L,
00083   1.1071487177940905030170654601785370400700E0L,
00084   1.1309537439791604464709335155363278047493E0L,
00085   1.1525719972156675180401498626127513797495E0L,
00086   1.1722738811284763866005949441337046149712E0L,
00087   1.1902899496825317329277337748293183376012E0L,
00088   1.2068173702852525303955115800565576303133E0L,
00089   1.2220253232109896370417417439225704908830E0L,
00090   1.2360594894780819419094519711090786987027E0L,
00091   1.2490457723982544258299170772810901230778E0L,
00092   1.2610933822524404193139408812473357720101E0L,
00093   1.2722973952087173412961937498224804940684E0L,
00094   1.2827408797442707473628852511364955306249E0L,
00095   1.2924966677897852679030914214070816845853E0L,
00096   1.3016288340091961438047858503666855921414E0L,
00097   1.3101939350475556342564376891719053122733E0L,
00098   1.3182420510168370498593302023271362531155E0L,
00099   1.3258176636680324650592392104284756311844E0L,
00100   1.3329603993374458675538498697331558093700E0L,
00101   1.3397056595989995393283037525895557411039E0L,
00102   1.3460851583802539310489409282517796256512E0L,
00103   1.3521273809209546571891479413898128509842E0L,
00104   1.3578579772154994751124898859640585287459E0L,
00105   1.3633001003596939542892985278250991189943E0L,
00106   1.3684746984165928776366381936948529556191E0L,
00107   1.3734007669450158608612719264449611486510E0L,
00108   1.3780955681325110444536609641291551522494E0L,
00109   1.3825748214901258580599674177685685125566E0L,
00110   1.3868528702577214543289381097042486034883E0L,
00111   1.3909428270024183486427686943836432060856E0L,
00112   1.3948567013423687823948122092044222644895E0L,
00113   1.3986055122719575950126700816114282335732E0L,
00114   1.4021993871854670105330304794336492676944E0L,
00115   1.4056476493802697809521934019958079881002E0L,
00116   1.4089588955564736949699075250792569287156E0L,
00117   1.4121410646084952153676136718584891599630E0L,
00118   1.4152014988178669079462550975833894394929E0L,
00119   1.4181469983996314594038603039700989523716E0L,
00120   1.4209838702219992566633046424614466661176E0L,
00121   1.4237179714064941189018190466107297503086E0L,
00122   1.4263547484202526397918060597281265695725E0L,
00123   1.4288992721907326964184700745371983590908E0L,
00124   1.4313562697035588982240194668401779312122E0L,
00125   1.4337301524847089866404719096698873648610E0L,
00126   1.4360250423171655234964275337155008780675E0L,
00127   1.4382447944982225979614042479354815855386E0L,
00128   1.4403930189057632173997301031392126865694E0L,
00129   1.4424730991091018200252920599377292525125E0L,
00130   1.4444882097316563655148453598508037025938E0L,
00131   1.4464413322481351841999668424758804165254E0L,
00132   1.4483352693775551917970437843145232637695E0L,
00133   1.4501726582147939000905940595923466567576E0L,
00134   1.4519559822271314199339700039142990228105E0L,
00135   1.4536875822280323362423034480994649820285E0L,
00136   1.4553696664279718992423082296859928222270E0L,
00137   1.4570043196511885530074841089245667532358E0L,
00138   1.4585935117976422128825857356750737658039E0L,
00139   1.4601391056210009726721818194296893361233E0L,
00140   1.4616428638860188872060496086383008594310E0L,
00141   1.4631064559620759326975975316301202111560E0L,
00142   1.4645314639038178118428450961503371619177E0L,
00143   1.4659193880646627234129855241049975398470E0L,
00144   1.4672716522843522691530527207287398276197E0L,
00145   1.4685896086876430842559640450619880951144E0L,
00146   1.4698745421276027686510391411132998919794E0L,
00147   1.4711276743037345918528755717617308518553E0L,
00148   1.4723501675822635384916444186631899205983E0L,
00149   1.4735431285433308455179928682541563973416E0L, /* arctan(10.25) */
00150   1.5707963267948966192313216916397514420986E0L  /* pi/2 */
00151 };
00152 
00153 
00154 /* arctan t = t + t^3 p(t^2) / q(t^2)
00155    |t| <= 0.09375
00156    peak relative error 5.3e-37 */
00157 
00158 static const long double
00159   p0 = -4.283708356338736809269381409828726405572E1L,
00160   p1 = -8.636132499244548540964557273544599863825E1L,
00161   p2 = -5.713554848244551350855604111031839613216E1L,
00162   p3 = -1.371405711877433266573835355036413750118E1L,
00163   p4 = -8.638214309119210906997318946650189640184E-1L,
00164   q0 = 1.285112506901621042780814422948906537959E2L,
00165   q1 = 3.361907253914337187957855834229672347089E2L,
00166   q2 = 3.180448303864130128268191635189365331680E2L,
00167   q3 = 1.307244136980865800160844625025280344686E2L,
00168   q4 = 2.173623741810414221251136181221172551416E1L;
00169   /* q5 = 1.000000000000000000000000000000000000000E0 */
00170 
00171 
00172 long double
00173 __atanl (long double x)
00174 {
00175   int k, sign;
00176   long double t, u, p, q;
00177   ieee854_long_double_shape_type s;
00178 
00179   s.value = x;
00180   k = s.parts32.w0;
00181   if (k & 0x80000000)
00182     sign = 1;
00183   else
00184     sign = 0;
00185 
00186   /* Check for IEEE special cases.  */
00187   k &= 0x7fffffff;
00188   if (k >= 0x7ff00000)
00189     {
00190       /* NaN. */
00191       if ((k & 0xfffff) | s.parts32.w1 )
00192        return (x + x);
00193 
00194       /* Infinity. */
00195       if (sign)
00196        return -atantbl[83];
00197       else
00198        return atantbl[83];
00199     }
00200 
00201   if (sign)
00202       x = -x;
00203 
00204   if (k >= 0x40248000) /* 10.25 */
00205     {
00206       k = 83;
00207       t = -1.0/x;
00208     }
00209   else
00210     {
00211       /* Index of nearest table element.
00212         Roundoff to integer is asymmetrical to avoid cancellation when t < 0
00213          (cf. fdlibm). */
00214       k = 8.0 * x + 0.25;
00215       u = 0.125 * k;
00216       /* Small arctan argument.  */
00217       t = (x - u) / (1.0 + x * u);
00218     }
00219 
00220   /* Arctan of small argument t.  */
00221   u = t * t;
00222   p =     ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
00223   q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
00224   u = t * u * p / q  +  t;
00225 
00226   /* arctan x = arctan u  +  arctan t */
00227   u = atantbl[k] + u;
00228   if (sign)
00229     return (-u);
00230   else
00231     return u;
00232 }
00233 
00234 long_double_symbol (libm, __atanl, atanl);