Blame sysdeps/ieee754/ldbl-128ibm/s_atanl.c

Packit Service 82fcde
/*							s_atanl.c
Packit Service 82fcde
 *
Packit Service 82fcde
 *	Inverse circular tangent for 128-bit long double precision
Packit Service 82fcde
 *      (arctangent)
Packit Service 82fcde
 *
Packit Service 82fcde
 *
Packit Service 82fcde
 *
Packit Service 82fcde
 * SYNOPSIS:
Packit Service 82fcde
 *
Packit Service 82fcde
 * long double x, y, atanl();
Packit Service 82fcde
 *
Packit Service 82fcde
 * y = atanl( x );
Packit Service 82fcde
 *
Packit Service 82fcde
 *
Packit Service 82fcde
 *
Packit Service 82fcde
 * DESCRIPTION:
Packit Service 82fcde
 *
Packit Service 82fcde
 * Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
Packit Service 82fcde
 *
Packit Service 82fcde
 * The function uses a rational approximation of the form
Packit Service 82fcde
 * t + t^3 P(t^2)/Q(t^2), optimized for |t| < 0.09375.
Packit Service 82fcde
 *
Packit Service 82fcde
 * The argument is reduced using the identity
Packit Service 82fcde
 *    arctan x - arctan u  =  arctan ((x-u)/(1 + ux))
Packit Service 82fcde
 * and an 83-entry lookup table for arctan u, with u = 0, 1/8, ..., 10.25.
Packit Service 82fcde
 * Use of the table improves the execution speed of the routine.
Packit Service 82fcde
 *
Packit Service 82fcde
 *
Packit Service 82fcde
 *
Packit Service 82fcde
 * ACCURACY:
Packit Service 82fcde
 *
Packit Service 82fcde
 *                      Relative error:
Packit Service 82fcde
 * arithmetic   domain     # trials      peak         rms
Packit Service 82fcde
 *    IEEE      -19, 19       4e5       1.7e-34     5.4e-35
Packit Service 82fcde
 *
Packit Service 82fcde
 *
Packit Service 82fcde
 * WARNING:
Packit Service 82fcde
 *
Packit Service 82fcde
 * This program uses integer operations on bit fields of floating-point
Packit Service 82fcde
 * numbers.  It does not work with data structures other than the
Packit Service 82fcde
 * structure assumed.
Packit Service 82fcde
 *
Packit Service 82fcde
 */
Packit Service 82fcde
Packit Service 82fcde
/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
Packit Service 82fcde
Packit Service 82fcde
    This library is free software; you can redistribute it and/or
Packit Service 82fcde
    modify it under the terms of the GNU Lesser General Public
Packit Service 82fcde
    License as published by the Free Software Foundation; either
Packit Service 82fcde
    version 2.1 of the License, or (at your option) any later version.
Packit Service 82fcde
Packit Service 82fcde
    This library is distributed in the hope that it will be useful,
Packit Service 82fcde
    but WITHOUT ANY WARRANTY; without even the implied warranty of
Packit Service 82fcde
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Packit Service 82fcde
    Lesser General Public License for more details.
Packit Service 82fcde
Packit Service 82fcde
    You should have received a copy of the GNU Lesser General Public
Packit Service 82fcde
    License along with this library; if not, see
Packit Service 82fcde
    <http://www.gnu.org/licenses/>.  */
Packit Service 82fcde
Packit Service 82fcde
Packit Service 82fcde
#include <float.h>
Packit Service 82fcde
#include <math.h>
Packit Service 82fcde
#include <math_private.h>
Packit Service 82fcde
#include <math-underflow.h>
Packit Service 82fcde
#include <math_ldbl_opt.h>
Packit Service 82fcde
Packit Service 82fcde
/* arctan(k/8), k = 0, ..., 82 */
Packit Service 82fcde
static const long double atantbl[84] = {
Packit Service 82fcde
  0.0000000000000000000000000000000000000000E0L,
Packit Service 82fcde
  1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125)  */
Packit Service 82fcde
  2.4497866312686415417208248121127581091414E-1L,
Packit Service 82fcde
  3.5877067027057222039592006392646049977698E-1L,
Packit Service 82fcde
  4.6364760900080611621425623146121440202854E-1L,
Packit Service 82fcde
  5.5859931534356243597150821640166127034645E-1L,
Packit Service 82fcde
  6.4350110879328438680280922871732263804151E-1L,
Packit Service 82fcde
  7.1882999962162450541701415152590465395142E-1L,
Packit Service 82fcde
  7.8539816339744830961566084581987572104929E-1L,
Packit Service 82fcde
  8.4415398611317100251784414827164750652594E-1L,
Packit Service 82fcde
  8.9605538457134395617480071802993782702458E-1L,
Packit Service 82fcde
  9.4200004037946366473793717053459358607166E-1L,
Packit Service 82fcde
  9.8279372324732906798571061101466601449688E-1L,
Packit Service 82fcde
  1.0191413442663497346383429170230636487744E0L,
Packit Service 82fcde
  1.0516502125483736674598673120862998296302E0L,
Packit Service 82fcde
  1.0808390005411683108871567292171998202703E0L,
Packit Service 82fcde
  1.1071487177940905030170654601785370400700E0L,
Packit Service 82fcde
  1.1309537439791604464709335155363278047493E0L,
Packit Service 82fcde
  1.1525719972156675180401498626127513797495E0L,
Packit Service 82fcde
  1.1722738811284763866005949441337046149712E0L,
Packit Service 82fcde
  1.1902899496825317329277337748293183376012E0L,
Packit Service 82fcde
  1.2068173702852525303955115800565576303133E0L,
Packit Service 82fcde
  1.2220253232109896370417417439225704908830E0L,
Packit Service 82fcde
  1.2360594894780819419094519711090786987027E0L,
Packit Service 82fcde
  1.2490457723982544258299170772810901230778E0L,
Packit Service 82fcde
  1.2610933822524404193139408812473357720101E0L,
Packit Service 82fcde
  1.2722973952087173412961937498224804940684E0L,
Packit Service 82fcde
  1.2827408797442707473628852511364955306249E0L,
Packit Service 82fcde
  1.2924966677897852679030914214070816845853E0L,
Packit Service 82fcde
  1.3016288340091961438047858503666855921414E0L,
Packit Service 82fcde
  1.3101939350475556342564376891719053122733E0L,
Packit Service 82fcde
  1.3182420510168370498593302023271362531155E0L,
Packit Service 82fcde
  1.3258176636680324650592392104284756311844E0L,
Packit Service 82fcde
  1.3329603993374458675538498697331558093700E0L,
Packit Service 82fcde
  1.3397056595989995393283037525895557411039E0L,
Packit Service 82fcde
  1.3460851583802539310489409282517796256512E0L,
Packit Service 82fcde
  1.3521273809209546571891479413898128509842E0L,
Packit Service 82fcde
  1.3578579772154994751124898859640585287459E0L,
Packit Service 82fcde
  1.3633001003596939542892985278250991189943E0L,
Packit Service 82fcde
  1.3684746984165928776366381936948529556191E0L,
Packit Service 82fcde
  1.3734007669450158608612719264449611486510E0L,
Packit Service 82fcde
  1.3780955681325110444536609641291551522494E0L,
Packit Service 82fcde
  1.3825748214901258580599674177685685125566E0L,
Packit Service 82fcde
  1.3868528702577214543289381097042486034883E0L,
Packit Service 82fcde
  1.3909428270024183486427686943836432060856E0L,
Packit Service 82fcde
  1.3948567013423687823948122092044222644895E0L,
Packit Service 82fcde
  1.3986055122719575950126700816114282335732E0L,
Packit Service 82fcde
  1.4021993871854670105330304794336492676944E0L,
Packit Service 82fcde
  1.4056476493802697809521934019958079881002E0L,
Packit Service 82fcde
  1.4089588955564736949699075250792569287156E0L,
Packit Service 82fcde
  1.4121410646084952153676136718584891599630E0L,
Packit Service 82fcde
  1.4152014988178669079462550975833894394929E0L,
Packit Service 82fcde
  1.4181469983996314594038603039700989523716E0L,
Packit Service 82fcde
  1.4209838702219992566633046424614466661176E0L,
Packit Service 82fcde
  1.4237179714064941189018190466107297503086E0L,
Packit Service 82fcde
  1.4263547484202526397918060597281265695725E0L,
Packit Service 82fcde
  1.4288992721907326964184700745371983590908E0L,
Packit Service 82fcde
  1.4313562697035588982240194668401779312122E0L,
Packit Service 82fcde
  1.4337301524847089866404719096698873648610E0L,
Packit Service 82fcde
  1.4360250423171655234964275337155008780675E0L,
Packit Service 82fcde
  1.4382447944982225979614042479354815855386E0L,
Packit Service 82fcde
  1.4403930189057632173997301031392126865694E0L,
Packit Service 82fcde
  1.4424730991091018200252920599377292525125E0L,
Packit Service 82fcde
  1.4444882097316563655148453598508037025938E0L,
Packit Service 82fcde
  1.4464413322481351841999668424758804165254E0L,
Packit Service 82fcde
  1.4483352693775551917970437843145232637695E0L,
Packit Service 82fcde
  1.4501726582147939000905940595923466567576E0L,
Packit Service 82fcde
  1.4519559822271314199339700039142990228105E0L,
Packit Service 82fcde
  1.4536875822280323362423034480994649820285E0L,
Packit Service 82fcde
  1.4553696664279718992423082296859928222270E0L,
Packit Service 82fcde
  1.4570043196511885530074841089245667532358E0L,
Packit Service 82fcde
  1.4585935117976422128825857356750737658039E0L,
Packit Service 82fcde
  1.4601391056210009726721818194296893361233E0L,
Packit Service 82fcde
  1.4616428638860188872060496086383008594310E0L,
Packit Service 82fcde
  1.4631064559620759326975975316301202111560E0L,
Packit Service 82fcde
  1.4645314639038178118428450961503371619177E0L,
Packit Service 82fcde
  1.4659193880646627234129855241049975398470E0L,
Packit Service 82fcde
  1.4672716522843522691530527207287398276197E0L,
Packit Service 82fcde
  1.4685896086876430842559640450619880951144E0L,
Packit Service 82fcde
  1.4698745421276027686510391411132998919794E0L,
Packit Service 82fcde
  1.4711276743037345918528755717617308518553E0L,
Packit Service 82fcde
  1.4723501675822635384916444186631899205983E0L,
Packit Service 82fcde
  1.4735431285433308455179928682541563973416E0L, /* arctan(10.25) */
Packit Service 82fcde
  1.5707963267948966192313216916397514420986E0L  /* pi/2 */
Packit Service 82fcde
};
Packit Service 82fcde
Packit Service 82fcde
Packit Service 82fcde
/* arctan t = t + t^3 p(t^2) / q(t^2)
Packit Service 82fcde
   |t| <= 0.09375
Packit Service 82fcde
   peak relative error 5.3e-37 */
Packit Service 82fcde
Packit Service 82fcde
static const long double
Packit Service 82fcde
  p0 = -4.283708356338736809269381409828726405572E1L,
Packit Service 82fcde
  p1 = -8.636132499244548540964557273544599863825E1L,
Packit Service 82fcde
  p2 = -5.713554848244551350855604111031839613216E1L,
Packit Service 82fcde
  p3 = -1.371405711877433266573835355036413750118E1L,
Packit Service 82fcde
  p4 = -8.638214309119210906997318946650189640184E-1L,
Packit Service 82fcde
  q0 = 1.285112506901621042780814422948906537959E2L,
Packit Service 82fcde
  q1 = 3.361907253914337187957855834229672347089E2L,
Packit Service 82fcde
  q2 = 3.180448303864130128268191635189365331680E2L,
Packit Service 82fcde
  q3 = 1.307244136980865800160844625025280344686E2L,
Packit Service 82fcde
  q4 = 2.173623741810414221251136181221172551416E1L;
Packit Service 82fcde
  /* q5 = 1.000000000000000000000000000000000000000E0 */
Packit Service 82fcde
Packit Service 82fcde
Packit Service 82fcde
long double
Packit Service 82fcde
__atanl (long double x)
Packit Service 82fcde
{
Packit Service 82fcde
  int32_t k, sign, lx;
Packit Service 82fcde
  long double t, u, p, q;
Packit Service 82fcde
  double xhi;
Packit Service 82fcde
Packit Service 82fcde
  xhi = ldbl_high (x);
Packit Service 82fcde
  EXTRACT_WORDS (k, lx, xhi);
Packit Service 82fcde
  sign = k & 0x80000000;
Packit Service 82fcde
Packit Service 82fcde
  /* Check for IEEE special cases.  */
Packit Service 82fcde
  k &= 0x7fffffff;
Packit Service 82fcde
  if (k >= 0x7ff00000)
Packit Service 82fcde
    {
Packit Service 82fcde
      /* NaN. */
Packit Service 82fcde
      if (((k - 0x7ff00000) | lx) != 0)
Packit Service 82fcde
	return (x + x);
Packit Service 82fcde
Packit Service 82fcde
      /* Infinity. */
Packit Service 82fcde
      if (sign)
Packit Service 82fcde
	return -atantbl[83];
Packit Service 82fcde
      else
Packit Service 82fcde
	return atantbl[83];
Packit Service 82fcde
    }
Packit Service 82fcde
Packit Service 82fcde
  if (k <= 0x3c800000) /* |x| <= 2**-55.  */
Packit Service 82fcde
    {
Packit Service 82fcde
      math_check_force_underflow (x);
Packit Service 82fcde
      /* Raise inexact.  */
Packit Service 82fcde
      if (1e300L + x > 0.0)
Packit Service 82fcde
	return x;
Packit Service 82fcde
    }
Packit Service 82fcde
Packit Service 82fcde
  if (k >= 0x46c00000) /* |x| >= 2**109.  */
Packit Service 82fcde
    {
Packit Service 82fcde
      /* Saturate result to {-,+}pi/2.  */
Packit Service 82fcde
      if (sign)
Packit Service 82fcde
	return -atantbl[83];
Packit Service 82fcde
      else
Packit Service 82fcde
	return atantbl[83];
Packit Service 82fcde
    }
Packit Service 82fcde
Packit Service 82fcde
  if (sign)
Packit Service 82fcde
      x = -x;
Packit Service 82fcde
Packit Service 82fcde
  if (k >= 0x40248000) /* 10.25 */
Packit Service 82fcde
    {
Packit Service 82fcde
      k = 83;
Packit Service 82fcde
      t = -1.0/x;
Packit Service 82fcde
    }
Packit Service 82fcde
  else
Packit Service 82fcde
    {
Packit Service 82fcde
      /* Index of nearest table element.
Packit Service 82fcde
	 Roundoff to integer is asymmetrical to avoid cancellation when t < 0
Packit Service 82fcde
         (cf. fdlibm). */
Packit Service 82fcde
      k = 8.0 * x + 0.25;
Packit Service 82fcde
      u = 0.125 * k;
Packit Service 82fcde
      /* Small arctan argument.  */
Packit Service 82fcde
      t = (x - u) / (1.0 + x * u);
Packit Service 82fcde
    }
Packit Service 82fcde
Packit Service 82fcde
  /* Arctan of small argument t.  */
Packit Service 82fcde
  u = t * t;
Packit Service 82fcde
  p =     ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
Packit Service 82fcde
  q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
Packit Service 82fcde
  u = t * u * p / q  +  t;
Packit Service 82fcde
Packit Service 82fcde
  /* arctan x = arctan u  +  arctan t */
Packit Service 82fcde
  u = atantbl[k] + u;
Packit Service 82fcde
  if (sign)
Packit Service 82fcde
    return (-u);
Packit Service 82fcde
  else
Packit Service 82fcde
    return u;
Packit Service 82fcde
}
Packit Service 82fcde
Packit Service 82fcde
long_double_symbol (libm, __atanl, atanl);