[thirdparty/glibc.git] / sysdeps / ieee754 / dbl-64 / s_atan.c

/*
 * IBM Accurate Mathematical Library
 * written by International Business Machines Corp.
 * Copyright (C) 2001-2015 Free Software Foundation, Inc.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, see <http://www.gnu.org/licenses/>.
 */
/************************************************************************/
/*  MODULE_NAME: atnat.c                                                */
/*                                                                      */
/*  FUNCTIONS:  uatan                                                   */
/*              atanMp                                                  */
/*              signArctan                                              */
/*                                                                      */
/*                                                                      */
/*  FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat.h                 */
/*                mpatan.c mpatan2.c mpsqrt.c                           */
/*                uatan.tbl                                             */
/*                                                                      */
/* An ultimate atan() routine. Given an IEEE double machine number x    */
/* it computes the correctly rounded (to nearest) value of atan(x).     */
/*                                                                      */
/* Assumption: Machine arithmetic operations are performed in           */
/* round to nearest mode of IEEE 754 standard.                          */
/*                                                                      */
/************************************************************************/

#include <dla.h>
#include "mpa.h"
#include "MathLib.h"
#include "uatan.tbl"
#include "atnat.h"
#include <math.h>
#include <stap-probe.h>

void __mpatan (mp_no *, mp_no *, int);	/* see definition in mpatan.c */
static double atanMp (double, const int[]);

  /* Fix the sign of y and return */
static double
__signArctan (double x, double y)
{
  return __copysign (y, x);
}


/* An ultimate atan() routine. Given an IEEE double machine number x,    */
/* routine computes the correctly rounded (to nearest) value of atan(x). */
double
atan (double x)
{
  double cor, s1, ss1, s2, ss2, t1, t2, t3, t7, t8, t9, t10, u, u2, u3,
	 v, vv, w, ww, y, yy, z, zz;
#ifndef DLA_FMS
  double t4, t5, t6;
#endif
  int i, ux, dx;
  static const int pr[M] = { 6, 8, 10, 32 };
  number num;

  num.d = x;
  ux = num.i[HIGH_HALF];
  dx = num.i[LOW_HALF];

  /* x=NaN */
  if (((ux & 0x7ff00000) == 0x7ff00000)
      && (((ux & 0x000fffff) | dx) != 0x00000000))
    return x + x;

  /* Regular values of x, including denormals +-0 and +-INF */
  u = (x < 0) ? -x : x;
  if (u < C)
    {
      if (u < B)
	{
	  if (u < A)
	    return x;
	  else
	    {			/* A <= u < B */
	      v = x * x;
	      yy = d11.d + v * d13.d;
	      yy = d9.d + v * yy;
	      yy = d7.d + v * yy;
	      yy = d5.d + v * yy;
	      yy = d3.d + v * yy;
	      yy *= x * v;

	      if ((y = x + (yy - U1 * x)) == x + (yy + U1 * x))
		return y;

	      EMULV (x, x, v, vv, t1, t2, t3, t4, t5);	/* v+vv=x^2 */

	      s1 = f17.d + v * f19.d;
	      s1 = f15.d + v * s1;
	      s1 = f13.d + v * s1;
	      s1 = f11.d + v * s1;
	      s1 *= v;

	      ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
	      MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	      ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
	      MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	      ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
	      MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	      ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
	      MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	      MUL2 (x, 0, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7,
		    t8);
	      ADD2 (x, 0, s2, ss2, s1, ss1, t1, t2);
	      if ((y = s1 + (ss1 - U5 * s1)) == s1 + (ss1 + U5 * s1))
		return y;

	      return atanMp (x, pr);
	    }
	}
      else
	{			/* B <= u < C */
	  i = (TWO52 + TWO8 * u) - TWO52;
	  i -= 16;
	  z = u - cij[i][0].d;
	  yy = cij[i][5].d + z * cij[i][6].d;
	  yy = cij[i][4].d + z * yy;
	  yy = cij[i][3].d + z * yy;
	  yy = cij[i][2].d + z * yy;
	  yy *= z;

	  t1 = cij[i][1].d;
	  if (i < 112)
	    {
	      if (i < 48)
		u2 = U21;	/* u < 1/4        */
	      else
		u2 = U22;
	    }			/* 1/4 <= u < 1/2 */
	  else
	    {
	      if (i < 176)
		u2 = U23;	/* 1/2 <= u < 3/4 */
	      else
		u2 = U24;
	    }			/* 3/4 <= u <= 1  */
	  if ((y = t1 + (yy - u2 * t1)) == t1 + (yy + u2 * t1))
	    return __signArctan (x, y);

	  z = u - hij[i][0].d;

	  s1 = hij[i][14].d + z * hij[i][15].d;
	  s1 = hij[i][13].d + z * s1;
	  s1 = hij[i][12].d + z * s1;
	  s1 = hij[i][11].d + z * s1;
	  s1 *= z;

	  ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
	  MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	  ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
	  MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	  ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
	  MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	  ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
	  MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	  ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
	  if ((y = s2 + (ss2 - U6 * s2)) == s2 + (ss2 + U6 * s2))
	    return __signArctan (x, y);

	  return atanMp (x, pr);
	}
    }
  else
    {
      if (u < D)
	{			/* C <= u < D */
	  w = 1 / u;
	  EMULV (w, u, t1, t2, t3, t4, t5, t6, t7);
	  ww = w * ((1 - t1) - t2);
	  i = (TWO52 + TWO8 * w) - TWO52;
	  i -= 16;
	  z = (w - cij[i][0].d) + ww;

	  yy = cij[i][5].d + z * cij[i][6].d;
	  yy = cij[i][4].d + z * yy;
	  yy = cij[i][3].d + z * yy;
	  yy = cij[i][2].d + z * yy;
	  yy = HPI1 - z * yy;

	  t1 = HPI - cij[i][1].d;
	  if (i < 112)
	    u3 = U31;           /* w <  1/2 */
	  else
	    u3 = U32;           /* w >= 1/2 */
	  if ((y = t1 + (yy - u3)) == t1 + (yy + u3))
	    return __signArctan (x, y);

	  DIV2 (1, 0, u, 0, w, ww, t1, t2, t3, t4, t5, t6, t7, t8, t9,
		t10);
	  t1 = w - hij[i][0].d;
	  EADD (t1, ww, z, zz);

	  s1 = hij[i][14].d + z * hij[i][15].d;
	  s1 = hij[i][13].d + z * s1;
	  s1 = hij[i][12].d + z * s1;
	  s1 = hij[i][11].d + z * s1;
	  s1 *= z;

	  ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
	  MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	  ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
	  MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	  ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
	  MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	  ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
	  MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	  ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
	  SUB2 (HPI, HPI1, s2, ss2, s1, ss1, t1, t2);
	  if ((y = s1 + (ss1 - U7)) == s1 + (ss1 + U7))
	    return __signArctan (x, y);

	  return atanMp (x, pr);
	}
      else
	{
	  if (u < E)
	    {                   /* D <= u < E */
	      w = 1 / u;
	      v = w * w;
	      EMULV (w, u, t1, t2, t3, t4, t5, t6, t7);

	      yy = d11.d + v * d13.d;
	      yy = d9.d + v * yy;
	      yy = d7.d + v * yy;
	      yy = d5.d + v * yy;
	      yy = d3.d + v * yy;
	      yy *= w * v;

	      ww = w * ((1 - t1) - t2);
	      ESUB (HPI, w, t3, cor);
	      yy = ((HPI1 + cor) - ww) - yy;
	      if ((y = t3 + (yy - U4)) == t3 + (yy + U4))
		return __signArctan (x, y);

	      DIV2 (1, 0, u, 0, w, ww, t1, t2, t3, t4, t5, t6, t7, t8,
		    t9, t10);
	      MUL2 (w, ww, w, ww, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);

	      s1 = f17.d + v * f19.d;
	      s1 = f15.d + v * s1;
	      s1 = f13.d + v * s1;
	      s1 = f11.d + v * s1;
	      s1 *= v;

	      ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
	      MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	      ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
	      MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	      ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
	      MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	      ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
	      MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	      MUL2 (w, ww, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
	      ADD2 (w, ww, s2, ss2, s1, ss1, t1, t2);
	      SUB2 (HPI, HPI1, s1, ss1, s2, ss2, t1, t2);

	      if ((y = s2 + (ss2 - U8)) == s2 + (ss2 + U8))
		return __signArctan (x, y);

	      return atanMp (x, pr);
	    }
	  else
	    {
	      /* u >= E */
	      if (x > 0)
		return HPI;
	      else
		return MHPI;
	    }
	}
    }
}

 /* Final stages. Compute atan(x) by multiple precision arithmetic */
static double
atanMp (double x, const int pr[])
{
  mp_no mpx, mpy, mpy2, mperr, mpt1, mpy1;
  double y1, y2;
  int i, p;

  for (i = 0; i < M; i++)
    {
      p = pr[i];
      __dbl_mp (x, &mpx, p);
      __mpatan (&mpx, &mpy, p);
      __dbl_mp (u9[i].d, &mpt1, p);
      __mul (&mpy, &mpt1, &mperr, p);
      __add (&mpy, &mperr, &mpy1, p);
      __sub (&mpy, &mperr, &mpy2, p);
      __mp_dbl (&mpy1, &y1, p);
      __mp_dbl (&mpy2, &y2, p);
      if (y1 == y2)
	{
	  LIBC_PROBE (slowatan, 3, &p, &x, &y1);
	  return y1;
	}
    }
  LIBC_PROBE (slowatan_inexact, 3, &p, &x, &y1);
  return y1;			/*if impossible to do exact computing */
}

#ifdef NO_LONG_DOUBLE
weak_alias (atan, atanl)
#endif
Commit	Line	Data
f7eac6eb	1	/*
e4d82761	2	* IBM Accurate Mathematical Library
aeb25823	3	* written by International Business Machines Corp.
b168057a	4	* Copyright (C) 2001-2015 Free Software Foundation, Inc.
f7eac6eb	5	*
e4d82761 UD	6	* This program is free software; you can redistribute it and/or modify
e4d82761 UD	7	* it under the terms of the GNU Lesser General Public License as published by
cc7375ce	8	* the Free Software Foundation; either version 2.1 of the License, or
e4d82761	9	* (at your option) any later version.
f7eac6eb	10	*
e4d82761 UD	11	* This program is distributed in the hope that it will be useful,
	12	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	13	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c6c6dd48	14	* GNU Lesser General Public License for more details.
f7eac6eb	15	*
e4d82761	16	* You should have received a copy of the GNU Lesser General Public License
59ba27a6	17	* along with this program; if not, see <http://www.gnu.org/licenses/>.
f7eac6eb	18	*/
e4d82761 UD	19	/************************************************************************/
	20	/* MODULE_NAME: atnat.c */
	21	/* */
	22	/* FUNCTIONS: uatan */
	23	/* atanMp */
	24	/* signArctan */
	25	/* */
	26	/* */
	27	/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat.h */
	28	/* mpatan.c mpatan2.c mpsqrt.c */
	29	/* uatan.tbl */
	30	/* */
	31	/* An ultimate atan() routine. Given an IEEE double machine number x */
	32	/* it computes the correctly rounded (to nearest) value of atan(x). */
	33	/* */
	34	/* Assumption: Machine arithmetic operations are performed in */
	35	/* round to nearest mode of IEEE 754 standard. */
	36	/* */
	37	/************************************************************************/
f7eac6eb	38
c8b3296b	39	#include <dla.h>
e4d82761 UD	40	#include "mpa.h"
	41	#include "MathLib.h"
	42	#include "uatan.tbl"
	43	#include "atnat.h"
1ed0291c	44	#include <math.h>
10e1cf6b	45	#include <stap-probe.h>
e4d82761	46
d26dd3eb SP	47	void __mpatan (mp_no , mp_no , int); /* see definition in mpatan.c */
d26dd3eb SP	48	static double atanMp (double, const int[]);
af968f62 UD	49
af968f62 UD	50	/* Fix the sign of y and return */
d26dd3eb SP	51	static double
	52	__signArctan (double x, double y)
	53	{
	54	return __copysign (y, x);
af968f62 UD	55	}
	56
	57
e4d82761 UD	58	/* An ultimate atan() routine. Given an IEEE double machine number x, */
e4d82761 UD	59	/* routine computes the correctly rounded (to nearest) value of atan(x). */
d26dd3eb SP	60	double
	61	atan (double x)
	62	{
	63	double cor, s1, ss1, s2, ss2, t1, t2, t3, t7, t8, t9, t10, u, u2, u3,
c5d5d574	64	v, vv, w, ww, y, yy, z, zz;
58985aa9	65	#ifndef DLA_FMS
d26dd3eb	66	double t4, t5, t6;
50944bca	67	#endif
d26dd3eb SP	68	int i, ux, dx;
d26dd3eb SP	69	static const int pr[M] = { 6, 8, 10, 32 };
e4d82761	70	number num;
e4d82761	71
d26dd3eb SP	72	num.d = x;
	73	ux = num.i[HIGH_HALF];
	74	dx = num.i[LOW_HALF];
e4d82761 UD	75
e4d82761 UD	76	/* x=NaN */
d26dd3eb SP	77	if (((ux & 0x7ff00000) == 0x7ff00000)
	78	&& (((ux & 0x000fffff) \| dx) != 0x00000000))
	79	return x + x;
e4d82761 UD	80
e4d82761 UD	81	/* Regular values of x, including denormals +-0 and +-INF */
a64d7e0e	82	u = (x < 0) ? -x : x;
d26dd3eb SP	83	if (u < C)
	84	{
	85	if (u < B)
	86	{
	87	if (u < A)
	88	return x;
	89	else
	90	{ /* A <= u < B */
	91	v = x * x;
	92	yy = d11.d + v * d13.d;
	93	yy = d9.d + v * yy;
	94	yy = d7.d + v * yy;
	95	yy = d5.d + v * yy;
	96	yy = d3.d + v * yy;
	97	yy = x v;
	98
	99	if ((y = x + (yy - U1 * x)) == x + (yy + U1 * x))
	100	return y;
	101
	102	EMULV (x, x, v, vv, t1, t2, t3, t4, t5); /* v+vv=x^2 */
	103
	104	s1 = f17.d + v * f19.d;
	105	s1 = f15.d + v * s1;
	106	s1 = f13.d + v * s1;
	107	s1 = f11.d + v * s1;
	108	s1 *= v;
	109
a64d7e0e	110	ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
d26dd3eb SP	111	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	112	ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
	113	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	114	ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
	115	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	116	ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
	117	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
a64d7e0e	118	MUL2 (x, 0, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7,
d26dd3eb	119	t8);
a64d7e0e	120	ADD2 (x, 0, s2, ss2, s1, ss1, t1, t2);
d26dd3eb SP	121	if ((y = s1 + (ss1 - U5 * s1)) == s1 + (ss1 + U5 * s1))
	122	return y;
	123
	124	return atanMp (x, pr);
	125	}
	126	}
	127	else
	128	{ /* B <= u < C */
	129	i = (TWO52 + TWO8 * u) - TWO52;
	130	i -= 16;
	131	z = u - cij[i][0].d;
	132	yy = cij[i][5].d + z * cij[i][6].d;
	133	yy = cij[i][4].d + z * yy;
	134	yy = cij[i][3].d + z * yy;
	135	yy = cij[i][2].d + z * yy;
	136	yy *= z;
	137
	138	t1 = cij[i][1].d;
	139	if (i < 112)
	140	{
	141	if (i < 48)
	142	u2 = U21; /* u < 1/4 */
	143	else
	144	u2 = U22;
	145	} /* 1/4 <= u < 1/2 */
	146	else
	147	{
	148	if (i < 176)
	149	u2 = U23; /* 1/2 <= u < 3/4 */
	150	else
	151	u2 = U24;
	152	} /* 3/4 <= u <= 1 */
	153	if ((y = t1 + (yy - u2 * t1)) == t1 + (yy + u2 * t1))
	154	return __signArctan (x, y);
	155
	156	z = u - hij[i][0].d;
	157
	158	s1 = hij[i][14].d + z * hij[i][15].d;
	159	s1 = hij[i][13].d + z * s1;
	160	s1 = hij[i][12].d + z * s1;
	161	s1 = hij[i][11].d + z * s1;
	162	s1 *= z;
	163
a64d7e0e SP	164	ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
a64d7e0e SP	165	MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
d26dd3eb	166	ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
a64d7e0e	167	MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
d26dd3eb	168	ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
a64d7e0e	169	MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
d26dd3eb	170	ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
a64d7e0e	171	MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
d26dd3eb SP	172	ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
	173	if ((y = s2 + (ss2 - U6 * s2)) == s2 + (ss2 + U6 * s2))
	174	return __signArctan (x, y);
	175
	176	return atanMp (x, pr);
	177	}
e4d82761	178	}
d26dd3eb SP	179	else
	180	{
	181	if (u < D)
	182	{ /* C <= u < D */
c2d94018	183	w = 1 / u;
d26dd3eb	184	EMULV (w, u, t1, t2, t3, t4, t5, t6, t7);
c2d94018	185	ww = w * ((1 - t1) - t2);
d26dd3eb SP	186	i = (TWO52 + TWO8 * w) - TWO52;
	187	i -= 16;
	188	z = (w - cij[i][0].d) + ww;
	189
	190	yy = cij[i][5].d + z * cij[i][6].d;
	191	yy = cij[i][4].d + z * yy;
	192	yy = cij[i][3].d + z * yy;
	193	yy = cij[i][2].d + z * yy;
c5d5d574	194	yy = HPI1 - z * yy;
f7eac6eb	195
d26dd3eb SP	196	t1 = HPI - cij[i][1].d;
d26dd3eb SP	197	if (i < 112)
c5d5d574	198	u3 = U31; /* w < 1/2 */
d26dd3eb	199	else
c5d5d574	200	u3 = U32; /* w >= 1/2 */
d26dd3eb SP	201	if ((y = t1 + (yy - u3)) == t1 + (yy + u3))
	202	return __signArctan (x, y);
	203
c5d5d574	204	DIV2 (1, 0, u, 0, w, ww, t1, t2, t3, t4, t5, t6, t7, t8, t9,
d26dd3eb SP	205	t10);
	206	t1 = w - hij[i][0].d;
	207	EADD (t1, ww, z, zz);
	208
	209	s1 = hij[i][14].d + z * hij[i][15].d;
	210	s1 = hij[i][13].d + z * s1;
	211	s1 = hij[i][12].d + z * s1;
	212	s1 = hij[i][11].d + z * s1;
	213	s1 *= z;
	214
a64d7e0e	215	ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
d26dd3eb SP	216	MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	217	ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
	218	MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	219	ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
	220	MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	221	ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
	222	MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	223	ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
	224	SUB2 (HPI, HPI1, s2, ss2, s1, ss1, t1, t2);
	225	if ((y = s1 + (ss1 - U7)) == s1 + (ss1 + U7))
	226	return __signArctan (x, y);
	227
	228	return atanMp (x, pr);
	229	}
	230	else
	231	{
	232	if (u < E)
c5d5d574	233	{ /* D <= u < E */
c2d94018	234	w = 1 / u;
d26dd3eb SP	235	v = w * w;
	236	EMULV (w, u, t1, t2, t3, t4, t5, t6, t7);
	237
	238	yy = d11.d + v * d13.d;
	239	yy = d9.d + v * yy;
	240	yy = d7.d + v * yy;
	241	yy = d5.d + v * yy;
	242	yy = d3.d + v * yy;
	243	yy = w v;
	244
c2d94018	245	ww = w * ((1 - t1) - t2);
d26dd3eb SP	246	ESUB (HPI, w, t3, cor);
	247	yy = ((HPI1 + cor) - ww) - yy;
	248	if ((y = t3 + (yy - U4)) == t3 + (yy + U4))
	249	return __signArctan (x, y);
	250
c5d5d574	251	DIV2 (1, 0, u, 0, w, ww, t1, t2, t3, t4, t5, t6, t7, t8,
d26dd3eb SP	252	t9, t10);
	253	MUL2 (w, ww, w, ww, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
	254
	255	s1 = f17.d + v * f19.d;
	256	s1 = f15.d + v * s1;
	257	s1 = f13.d + v * s1;
	258	s1 = f11.d + v * s1;
	259	s1 *= v;
	260
a64d7e0e	261	ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
d26dd3eb SP	262	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	263	ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
	264	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	265	ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
	266	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	267	ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
	268	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	269	MUL2 (w, ww, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
	270	ADD2 (w, ww, s2, ss2, s1, ss1, t1, t2);
	271	SUB2 (HPI, HPI1, s1, ss1, s2, ss2, t1, t2);
	272
	273	if ((y = s2 + (ss2 - U8)) == s2 + (ss2 + U8))
	274	return __signArctan (x, y);
	275
	276	return atanMp (x, pr);
	277	}
	278	else
	279	{
	280	/* u >= E */
	281	if (x > 0)
	282	return HPI;
	283	else
	284	return MHPI;
	285	}
	286	}
	287	}
f7eac6eb	288	}
e4d82761	289
e4d82761	290	/* Final stages. Compute atan(x) by multiple precision arithmetic */
d26dd3eb SP	291	static double
	292	atanMp (double x, const int pr[])
	293	{
	294	mp_no mpx, mpy, mpy2, mperr, mpt1, mpy1;
	295	double y1, y2;
	296	int i, p;
	297
	298	for (i = 0; i < M; i++)
	299	{
	300	p = pr[i];
	301	__dbl_mp (x, &mpx, p);
	302	__mpatan (&mpx, &mpy, p);
	303	__dbl_mp (u9[i].d, &mpt1, p);
	304	__mul (&mpy, &mpt1, &mperr, p);
	305	__add (&mpy, &mperr, &mpy1, p);
	306	__sub (&mpy, &mperr, &mpy2, p);
	307	__mp_dbl (&mpy1, &y1, p);
	308	__mp_dbl (&mpy2, &y2, p);
	309	if (y1 == y2)
10e1cf6b SP	310	{
	311	LIBC_PROBE (slowatan, 3, &p, &x, &y1);
	312	return y1;
	313	}
d26dd3eb	314	}
10e1cf6b	315	LIBC_PROBE (slowatan_inexact, 3, &p, &x, &y1);
d26dd3eb	316	return y1; /if impossible to do exact computing /
e4d82761 UD	317	}
e4d82761 UD	318
cccda09f	319	#ifdef NO_LONG_DOUBLE
e4d82761	320	weak_alias (atan, atanl)
cccda09f	321	#endif