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

/*
 * IBM Accurate Mathematical Library
 * written by International Business Machines Corp.
 * Copyright (C) 2001-2014 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: atnat2.c                                               */
/*                                                                      */
/*  FUNCTIONS: uatan2                                                   */
/*             atan2Mp                                                  */
/*             signArctan2                                              */
/*             normalized                                               */
/*                                                                      */
/*  FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat2.h                */
/*                mpatan.c mpatan2.c mpsqrt.c                           */
/*                uatan.tbl                                             */
/*                                                                      */
/* An ultimate atan2() routine. Given two IEEE double machine numbers y,*/
/* x it computes the correctly rounded (to nearest) value of atan2(y,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 "atnat2.h"
#include <math_private.h>
#include <stap-probe.h>

#ifndef SECTION
# define SECTION
#endif

/************************************************************************/
/* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */
/* it computes the correctly rounded (to nearest) value of atan2(y,x).  */
/* Assumption: Machine arithmetic operations are performed in           */
/* round to nearest mode of IEEE 754 standard.                          */
/************************************************************************/
static double atan2Mp (double, double, const int[]);
  /* Fix the sign and return after stage 1 or stage 2 */
static double
signArctan2 (double y, double z)
{
  return __copysign (z, y);
}

static double normalized (double, double, double, double);
void __mpatan2 (mp_no *, mp_no *, mp_no *, int);

double
SECTION
__ieee754_atan2 (double y, double x)
{
  int i, de, ux, dx, uy, dy;
  static const int pr[MM] = { 6, 8, 10, 20, 32 };
  double ax, ay, u, du, u9, ua, v, vv, dv, t1, t2, t3, t7, t8,
	 z, zz, cor, s1, ss1, s2, ss2;
#ifndef DLA_FMS
  double t4, t5, t6;
#endif
  number num;

  static const int ep = 59768832,      /*  57*16**5   */
		   em = -59768832;      /* -57*16**5   */

  /* x=NaN or y=NaN */
  num.d = x;
  ux = num.i[HIGH_HALF];
  dx = num.i[LOW_HALF];
  if ((ux & 0x7ff00000) == 0x7ff00000)
    {
      if (((ux & 0x000fffff) | dx) != 0x00000000)
	return x + x;
    }
  num.d = y;
  uy = num.i[HIGH_HALF];
  dy = num.i[LOW_HALF];
  if ((uy & 0x7ff00000) == 0x7ff00000)
    {
      if (((uy & 0x000fffff) | dy) != 0x00000000)
	return y + y;
    }

  /* y=+-0 */
  if (uy == 0x00000000)
    {
      if (dy == 0x00000000)
	{
	  if ((ux & 0x80000000) == 0x00000000)
	    return 0;
	  else
	    return opi.d;
	}
    }
  else if (uy == 0x80000000)
    {
      if (dy == 0x00000000)
	{
	  if ((ux & 0x80000000) == 0x00000000)
	    return -0.0;
	  else
	    return mopi.d;
	}
    }

  /* x=+-0 */
  if (x == 0)
    {
      if ((uy & 0x80000000) == 0x00000000)
	return hpi.d;
      else
	return mhpi.d;
    }

  /* x=+-INF */
  if (ux == 0x7ff00000)
    {
      if (dx == 0x00000000)
	{
	  if (uy == 0x7ff00000)
	    {
	      if (dy == 0x00000000)
		return qpi.d;
	    }
	  else if (uy == 0xfff00000)
	    {
	      if (dy == 0x00000000)
		return mqpi.d;
	    }
	  else
	    {
	      if ((uy & 0x80000000) == 0x00000000)
		return 0;
	      else
		return -0.0;
	    }
	}
    }
  else if (ux == 0xfff00000)
    {
      if (dx == 0x00000000)
	{
	  if (uy == 0x7ff00000)
	    {
	      if (dy == 0x00000000)
		return tqpi.d;
	    }
	  else if (uy == 0xfff00000)
	    {
	      if (dy == 0x00000000)
		return mtqpi.d;
	    }
	  else
	    {
	      if ((uy & 0x80000000) == 0x00000000)
		return opi.d;
	      else
		return mopi.d;
	    }
	}
    }

  /* y=+-INF */
  if (uy == 0x7ff00000)
    {
      if (dy == 0x00000000)
	return hpi.d;
    }
  else if (uy == 0xfff00000)
    {
      if (dy == 0x00000000)
	return mhpi.d;
    }

  /* either x/y or y/x is very close to zero */
  ax = (x < 0) ? -x : x;
  ay = (y < 0) ? -y : y;
  de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
  if (de >= ep)
    {
      return ((y > 0) ? hpi.d : mhpi.d);
    }
  else if (de <= em)
    {
      if (x > 0)
	{
	  if ((z = ay / ax) < TWOM1022)
	    return normalized (ax, ay, y, z);
	  else
	    return signArctan2 (y, z);
	}
      else
	{
	  return ((y > 0) ? opi.d : mopi.d);
	}
    }

  /* if either x or y is extremely close to zero, scale abs(x), abs(y). */
  if (ax < twom500.d || ay < twom500.d)
    {
      ax *= two500.d;
      ay *= two500.d;
    }

  /* Likewise for large x and y.  */
  if (ax > two500.d || ay > two500.d)
    {
      ax *= twom500.d;
      ay *= twom500.d;
    }

  /* x,y which are neither special nor extreme */
  if (ay < ax)
    {
      u = ay / ax;
      EMULV (ax, u, v, vv, t1, t2, t3, t4, t5);
      du = ((ay - v) - vv) / ax;
    }
  else
    {
      u = ax / ay;
      EMULV (ay, u, v, vv, t1, t2, t3, t4, t5);
      du = ((ax - v) - vv) / ay;
    }

  if (x > 0)
    {
      /* (i)   x>0, abs(y)< abs(x):  atan(ay/ax) */
      if (ay < ax)
	{
	  if (u < inv16.d)
	    {
	      v = u * u;

	      zz = du + u * v * (d3.d
				 + v * (d5.d
					+ v * (d7.d
					       + v * (d9.d
						      + v * (d11.d
							     + v * d13.d)))));

	      if ((z = u + (zz - u1.d * u)) == u + (zz + u1.d * u))
		return signArctan2 (y, z);

	      MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
	      s1 = v * (f11.d + v * (f13.d
				     + v * (f15.d + v * (f17.d + v * f19.d))));
	      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 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
	      ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);

	      if ((z = s1 + (ss1 - u5.d * s1)) == s1 + (ss1 + u5.d * s1))
		return signArctan2 (y, z);

	      return atan2Mp (x, y, pr);
	    }

	  i = (TWO52 + TWO8 * u) - TWO52;
	  i -= 16;
	  t3 = u - cij[i][0].d;
	  EADD (t3, du, v, dv);
	  t1 = cij[i][1].d;
	  t2 = cij[i][2].d;
	  zz = v * t2 + (dv * t2
			 + v * v * (cij[i][3].d
				    + v * (cij[i][4].d
					   + v * (cij[i][5].d
						  + v * cij[i][6].d))));
	  if (i < 112)
	    {
	      if (i < 48)
		u9 = u91.d;     /* u < 1/4	*/
	      else
		u9 = u92.d;
	    }           /* 1/4 <= u < 1/2 */
	  else
	    {
	      if (i < 176)
		u9 = u93.d;     /* 1/2 <= u < 3/4 */
	      else
		u9 = u94.d;
	    }           /* 3/4 <= u <= 1  */
	  if ((z = t1 + (zz - u9 * t1)) == t1 + (zz + u9 * t1))
	    return signArctan2 (y, z);

	  t1 = u - hij[i][0].d;
	  EADD (t1, du, v, vv);
	  s1 = v * (hij[i][11].d
		    + v * (hij[i][12].d
			   + v * (hij[i][13].d
				  + v * (hij[i][14].d
					 + v * hij[i][15].d))));
	  ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
	  MUL2 (v, vv, 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 (v, vv, 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 (v, vv, 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 (v, vv, 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 ((z = s2 + (ss2 - ub.d * s2)) == s2 + (ss2 + ub.d * s2))
	    return signArctan2 (y, z);
	  return atan2Mp (x, y, pr);
	}

      /* (ii)  x>0, abs(x)<=abs(y):  pi/2-atan(ax/ay) */
      if (u < inv16.d)
	{
	  v = u * u;
	  zz = u * v * (d3.d
			+ v * (d5.d
			       + v * (d7.d
				      + v * (d9.d
					     + v * (d11.d
						    + v * d13.d)))));
	  ESUB (hpi.d, u, t2, cor);
	  t3 = ((hpi1.d + cor) - du) - zz;
	  if ((z = t2 + (t3 - u2.d)) == t2 + (t3 + u2.d))
	    return signArctan2 (y, z);

	  MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
	  s1 = v * (f11.d
		    + v * (f13.d
			   + v * (f15.d + v * (f17.d + v * f19.d))));
	  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 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
	  ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);
	  SUB2 (hpi.d, hpi1.d, s1, ss1, s2, ss2, t1, t2);

	  if ((z = s2 + (ss2 - u6.d)) == s2 + (ss2 + u6.d))
	    return signArctan2 (y, z);
	  return atan2Mp (x, y, pr);
	}

      i = (TWO52 + TWO8 * u) - TWO52;
      i -= 16;
      v = (u - cij[i][0].d) + du;

      zz = hpi1.d - v * (cij[i][2].d
			 + v * (cij[i][3].d
				+ v * (cij[i][4].d
				       + v * (cij[i][5].d
					      + v * cij[i][6].d))));
      t1 = hpi.d - cij[i][1].d;
      if (i < 112)
	ua = ua1.d;	/* w <  1/2 */
      else
	ua = ua2.d;	/* w >= 1/2 */
      if ((z = t1 + (zz - ua)) == t1 + (zz + ua))
	return signArctan2 (y, z);

      t1 = u - hij[i][0].d;
      EADD (t1, du, v, vv);

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

      ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
      MUL2 (v, vv, 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 (v, vv, 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 (v, vv, 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 (v, vv, 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.d, hpi1.d, s2, ss2, s1, ss1, t1, t2);

      if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d))
	return signArctan2 (y, z);
      return atan2Mp (x, y, pr);
    }

  /* (iii) x<0, abs(x)< abs(y):  pi/2+atan(ax/ay) */
  if (ax < ay)
    {
      if (u < inv16.d)
	{
	  v = u * u;
	  zz = u * v * (d3.d
			+ v * (d5.d
			       + v * (d7.d
				      + v * (d9.d
					     + v * (d11.d + v * d13.d)))));
	  EADD (hpi.d, u, t2, cor);
	  t3 = ((hpi1.d + cor) + du) + zz;
	  if ((z = t2 + (t3 - u3.d)) == t2 + (t3 + u3.d))
	    return signArctan2 (y, z);

	  MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
	  s1 = v * (f11.d
		    + v * (f13.d + v * (f15.d + v * (f17.d + v * f19.d))));
	  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 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
	  ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);
	  ADD2 (hpi.d, hpi1.d, s1, ss1, s2, ss2, t1, t2);

	  if ((z = s2 + (ss2 - u7.d)) == s2 + (ss2 + u7.d))
	    return signArctan2 (y, z);
	  return atan2Mp (x, y, pr);
	}

      i = (TWO52 + TWO8 * u) - TWO52;
      i -= 16;
      v = (u - cij[i][0].d) + du;
      zz = hpi1.d + v * (cij[i][2].d
			 + v * (cij[i][3].d
				+ v * (cij[i][4].d
				       + v * (cij[i][5].d
					      + v * cij[i][6].d))));
      t1 = hpi.d + cij[i][1].d;
      if (i < 112)
	ua = ua1.d;	/* w <  1/2 */
      else
	ua = ua2.d;	/* w >= 1/2 */
      if ((z = t1 + (zz - ua)) == t1 + (zz + ua))
	return signArctan2 (y, z);

      t1 = u - hij[i][0].d;
      EADD (t1, du, v, vv);
      s1 = v * (hij[i][11].d
		+ v * (hij[i][12].d
		       + v * (hij[i][13].d
			      + v * (hij[i][14].d
				     + v * hij[i][15].d))));
      ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
      MUL2 (v, vv, 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 (v, vv, 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 (v, vv, 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 (v, vv, 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);
      ADD2 (hpi.d, hpi1.d, s2, ss2, s1, ss1, t1, t2);

      if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d))
	return signArctan2 (y, z);
      return atan2Mp (x, y, pr);
    }

  /* (iv)  x<0, abs(y)<=abs(x):  pi-atan(ax/ay) */
  if (u < inv16.d)
    {
      v = u * u;
      zz = u * v * (d3.d
		    + v * (d5.d
			   + v * (d7.d
				  + v * (d9.d + v * (d11.d + v * d13.d)))));
      ESUB (opi.d, u, t2, cor);
      t3 = ((opi1.d + cor) - du) - zz;
      if ((z = t2 + (t3 - u4.d)) == t2 + (t3 + u4.d))
	return signArctan2 (y, z);

      MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
      s1 = v * (f11.d + v * (f13.d + v * (f15.d + v * (f17.d + v * f19.d))));
      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 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
      ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);
      SUB2 (opi.d, opi1.d, s1, ss1, s2, ss2, t1, t2);

      if ((z = s2 + (ss2 - u8.d)) == s2 + (ss2 + u8.d))
	return signArctan2 (y, z);
      return atan2Mp (x, y, pr);
    }

  i = (TWO52 + TWO8 * u) - TWO52;
  i -= 16;
  v = (u - cij[i][0].d) + du;
  zz = opi1.d - v * (cij[i][2].d
		     + v * (cij[i][3].d
			    + v * (cij[i][4].d
				   + v * (cij[i][5].d + v * cij[i][6].d))));
  t1 = opi.d - cij[i][1].d;
  if (i < 112)
    ua = ua1.d;	/* w <  1/2 */
  else
    ua = ua2.d;	/* w >= 1/2 */
  if ((z = t1 + (zz - ua)) == t1 + (zz + ua))
    return signArctan2 (y, z);

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

  EADD (t1, du, v, vv);

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

  ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
  MUL2 (v, vv, 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 (v, vv, 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 (v, vv, 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 (v, vv, 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 (opi.d, opi1.d, s2, ss2, s1, ss1, t1, t2);

  if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d))
    return signArctan2 (y, z);
  return atan2Mp (x, y, pr);
}

#ifndef __ieee754_atan2
strong_alias (__ieee754_atan2, __atan2_finite)
#endif

/* Treat the Denormalized case */
static double
SECTION
normalized (double ax, double ay, double y, double z)
{
  int p;
  mp_no mpx, mpy, mpz, mperr, mpz2, mpt1;
  p = 6;
  __dbl_mp (ax, &mpx, p);
  __dbl_mp (ay, &mpy, p);
  __dvd (&mpy, &mpx, &mpz, p);
  __dbl_mp (ue.d, &mpt1, p);
  __mul (&mpz, &mpt1, &mperr, p);
  __sub (&mpz, &mperr, &mpz2, p);
  __mp_dbl (&mpz2, &z, p);
  return signArctan2 (y, z);
}

/* Stage 3: Perform a multi-Precision computation */
static double
SECTION
atan2Mp (double x, double y, const int pr[])
{
  double z1, z2;
  int i, p;
  mp_no mpx, mpy, mpz, mpz1, mpz2, mperr, mpt1;
  for (i = 0; i < MM; i++)
    {
      p = pr[i];
      __dbl_mp (x, &mpx, p);
      __dbl_mp (y, &mpy, p);
      __mpatan2 (&mpy, &mpx, &mpz, p);
      __dbl_mp (ud[i].d, &mpt1, p);
      __mul (&mpz, &mpt1, &mperr, p);
      __add (&mpz, &mperr, &mpz1, p);
      __sub (&mpz, &mperr, &mpz2, p);
      __mp_dbl (&mpz1, &z1, p);
      __mp_dbl (&mpz2, &z2, p);
      if (z1 == z2)
	{
	  LIBC_PROBE (slowatan2, 4, &p, &x, &y, &z1);
	  return z1;
	}
    }
  LIBC_PROBE (slowatan2_inexact, 4, &p, &x, &y, &z1);
  return z1;			/*if impossible to do exact computing */
}
Commit	Line	Data
f7eac6eb	1	/*
e4d82761	2	* IBM Accurate Mathematical Library
aeb25823	3	* written by International Business Machines Corp.
d4697bc9	4	* Copyright (C) 2001-2014 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: atnat2.c */
	21	/* */
	22	/* FUNCTIONS: uatan2 */
	23	/* atan2Mp */
	24	/* signArctan2 */
	25	/* normalized */
	26	/* */
	27	/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat2.h */
	28	/* mpatan.c mpatan2.c mpsqrt.c */
	29	/* uatan.tbl */
	30	/* */
	31	/* An ultimate atan2() routine. Given two IEEE double machine numbers y,*/
	32	/* x it computes the correctly rounded (to nearest) value of atan2(y,x).*/
	33	/* */
	34	/* Assumption: Machine arithmetic operations are performed in */
	35	/* round to nearest mode of IEEE 754 standard. */
	36	/* */
	37	/************************************************************************/
	38
c8b3296b	39	#include <dla.h>
e4d82761 UD	40	#include "mpa.h"
	41	#include "MathLib.h"
	42	#include "uatan.tbl"
	43	#include "atnat2.h"
1ed0291c	44	#include <math_private.h>
10e1cf6b	45	#include <stap-probe.h>
e859d1d9	46
31d3cc00 UD	47	#ifndef SECTION
	48	# define SECTION
	49	#endif
	50
e4d82761 UD	51	/************************************************************************/
	52	/* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */
	53	/* it computes the correctly rounded (to nearest) value of atan2(y,x). */
	54	/* Assumption: Machine arithmetic operations are performed in */
	55	/* round to nearest mode of IEEE 754 standard. */
	56	/************************************************************************/
1728ab37	57	static double atan2Mp (double, double, const int[]);
af968f62	58	/* Fix the sign and return after stage 1 or stage 2 */
1728ab37 SP	59	static double
1728ab37 SP	60	signArctan2 (double y, double z)
af968f62	61	{
1728ab37	62	return __copysign (z, y);
af968f62	63	}
1728ab37 SP	64
	65	static double normalized (double, double, double, double);
	66	void __mpatan2 (mp_no , mp_no , mp_no *, int);
e4d82761	67
31d3cc00 UD	68	double
31d3cc00 UD	69	SECTION
1728ab37 SP	70	__ieee754_atan2 (double y, double x)
	71	{
	72	int i, de, ux, dx, uy, dy;
	73	static const int pr[MM] = { 6, 8, 10, 20, 32 };
	74	double ax, ay, u, du, u9, ua, v, vv, dv, t1, t2, t3, t7, t8,
c5d5d574	75	z, zz, cor, s1, ss1, s2, ss2;
58985aa9	76	#ifndef DLA_FMS
1728ab37	77	double t4, t5, t6;
50944bca	78	#endif
e4d82761	79	number num;
e4d82761	80
c5d5d574 OB	81	static const int ep = 59768832, /* 57165 /
c5d5d574 OB	82	em = -59768832; /* -57165 /
e4d82761 UD	83
e4d82761 UD	84	/* x=NaN or y=NaN */
1728ab37 SP	85	num.d = x;
	86	ux = num.i[HIGH_HALF];
	87	dx = num.i[LOW_HALF];
	88	if ((ux & 0x7ff00000) == 0x7ff00000)
	89	{
	90	if (((ux & 0x000fffff) \| dx) != 0x00000000)
	91	return x + x;
	92	}
	93	num.d = y;
	94	uy = num.i[HIGH_HALF];
	95	dy = num.i[LOW_HALF];
	96	if ((uy & 0x7ff00000) == 0x7ff00000)
	97	{
	98	if (((uy & 0x000fffff) \| dy) != 0x00000000)
	99	return y + y;
	100	}
e4d82761 UD	101
e4d82761 UD	102	/* y=+-0 */
1728ab37 SP	103	if (uy == 0x00000000)
	104	{
	105	if (dy == 0x00000000)
	106	{
	107	if ((ux & 0x80000000) == 0x00000000)
a64d7e0e	108	return 0;
1728ab37 SP	109	else
	110	return opi.d;
	111	}
	112	}
	113	else if (uy == 0x80000000)
	114	{
	115	if (dy == 0x00000000)
	116	{
	117	if ((ux & 0x80000000) == 0x00000000)
a64d7e0e	118	return -0.0;
1728ab37 SP	119	else
	120	return mopi.d;
	121	}
	122	}
e4d82761 UD	123
e4d82761 UD	124	/* x=+-0 */
a64d7e0e	125	if (x == 0)
1728ab37 SP	126	{
	127	if ((uy & 0x80000000) == 0x00000000)
	128	return hpi.d;
	129	else
	130	return mhpi.d;
	131	}
e4d82761 UD	132
e4d82761 UD	133	/* x=+-INF */
1728ab37 SP	134	if (ux == 0x7ff00000)
	135	{
	136	if (dx == 0x00000000)
	137	{
	138	if (uy == 0x7ff00000)
	139	{
	140	if (dy == 0x00000000)
	141	return qpi.d;
	142	}
	143	else if (uy == 0xfff00000)
	144	{
	145	if (dy == 0x00000000)
	146	return mqpi.d;
	147	}
	148	else
	149	{
	150	if ((uy & 0x80000000) == 0x00000000)
a64d7e0e	151	return 0;
1728ab37	152	else
a64d7e0e	153	return -0.0;
1728ab37 SP	154	}
1728ab37 SP	155	}
e4d82761	156	}
1728ab37 SP	157	else if (ux == 0xfff00000)
	158	{
	159	if (dx == 0x00000000)
	160	{
	161	if (uy == 0x7ff00000)
	162	{
	163	if (dy == 0x00000000)
	164	return tqpi.d;
	165	}
	166	else if (uy == 0xfff00000)
	167	{
	168	if (dy == 0x00000000)
	169	return mtqpi.d;
	170	}
	171	else
	172	{
	173	if ((uy & 0x80000000) == 0x00000000)
	174	return opi.d;
	175	else
	176	return mopi.d;
	177	}
	178	}
e4d82761	179	}
e4d82761 UD	180
e4d82761 UD	181	/* y=+-INF */
1728ab37 SP	182	if (uy == 0x7ff00000)
	183	{
	184	if (dy == 0x00000000)
	185	return hpi.d;
	186	}
	187	else if (uy == 0xfff00000)
	188	{
	189	if (dy == 0x00000000)
	190	return mhpi.d;
	191	}
e4d82761 UD	192
e4d82761 UD	193	/* either x/y or y/x is very close to zero */
a64d7e0e SP	194	ax = (x < 0) ? -x : x;
a64d7e0e SP	195	ay = (y < 0) ? -y : y;
e4d82761	196	de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
1728ab37 SP	197	if (de >= ep)
1728ab37 SP	198	{
a64d7e0e	199	return ((y > 0) ? hpi.d : mhpi.d);
1728ab37 SP	200	}
	201	else if (de <= em)
	202	{
a64d7e0e	203	if (x > 0)
1728ab37 SP	204	{
	205	if ((z = ay / ax) < TWOM1022)
	206	return normalized (ax, ay, y, z);
	207	else
	208	return signArctan2 (y, z);
	209	}
	210	else
	211	{
a64d7e0e	212	return ((y > 0) ? opi.d : mopi.d);
1728ab37 SP	213	}
1728ab37 SP	214	}
e4d82761 UD	215
e4d82761 UD	216	/* if either x or y is extremely close to zero, scale abs(x), abs(y). */
1728ab37 SP	217	if (ax < twom500.d \|\| ay < twom500.d)
	218	{
	219	ax *= two500.d;
	220	ay *= two500.d;
	221	}
e4d82761	222
7726d6a9 JM	223	/* Likewise for large x and y. */
	224	if (ax > two500.d \|\| ay > two500.d)
	225	{
	226	ax *= twom500.d;
	227	ay *= twom500.d;
	228	}
	229
e4d82761	230	/* x,y which are neither special nor extreme */
1728ab37 SP	231	if (ay < ax)
	232	{
	233	u = ay / ax;
	234	EMULV (ax, u, v, vv, t1, t2, t3, t4, t5);
	235	du = ((ay - v) - vv) / ax;
	236	}
	237	else
	238	{
	239	u = ax / ay;
	240	EMULV (ay, u, v, vv, t1, t2, t3, t4, t5);
	241	du = ((ax - v) - vv) / ay;
e4d82761 UD	242	}
e4d82761 UD	243
a64d7e0e	244	if (x > 0)
1728ab37 SP	245	{
	246	/* (i) x>0, abs(y)< abs(x): atan(ay/ax) */
	247	if (ay < ax)
	248	{
	249	if (u < inv16.d)
	250	{
	251	v = u * u;
	252
	253	zz = du + u * v * (d3.d
	254	+ v * (d5.d
	255	+ v * (d7.d
	256	+ v * (d9.d
	257	+ v * (d11.d
	258	+ v * d13.d)))));
	259
	260	if ((z = u + (zz - u1.d * u)) == u + (zz + u1.d * u))
	261	return signArctan2 (y, z);
	262
	263	MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
	264	s1 = v * (f11.d + v * (f13.d
	265	+ v * (f15.d + v * (f17.d + v * f19.d))));
a64d7e0e	266	ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
1728ab37 SP	267	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	268	ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
	269	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	270	ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
	271	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	272	ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
	273	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	274	MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
	275	ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);
	276
	277	if ((z = s1 + (ss1 - u5.d * s1)) == s1 + (ss1 + u5.d * s1))
	278	return signArctan2 (y, z);
	279
	280	return atan2Mp (x, y, pr);
	281	}
	282
	283	i = (TWO52 + TWO8 * u) - TWO52;
	284	i -= 16;
	285	t3 = u - cij[i][0].d;
	286	EADD (t3, du, v, dv);
	287	t1 = cij[i][1].d;
	288	t2 = cij[i][2].d;
	289	zz = v * t2 + (dv * t2
	290	+ v * v * (cij[i][3].d
	291	+ v * (cij[i][4].d
	292	+ v * (cij[i][5].d
	293	+ v * cij[i][6].d))));
	294	if (i < 112)
	295	{
	296	if (i < 48)
c5d5d574	297	u9 = u91.d; /* u < 1/4 */
1728ab37 SP	298	else
1728ab37 SP	299	u9 = u92.d;
c5d5d574	300	} /* 1/4 <= u < 1/2 */
1728ab37 SP	301	else
	302	{
	303	if (i < 176)
c5d5d574	304	u9 = u93.d; /* 1/2 <= u < 3/4 */
1728ab37 SP	305	else
1728ab37 SP	306	u9 = u94.d;
c5d5d574	307	} /* 3/4 <= u <= 1 */
1728ab37 SP	308	if ((z = t1 + (zz - u9 * t1)) == t1 + (zz + u9 * t1))
	309	return signArctan2 (y, z);
	310
	311	t1 = u - hij[i][0].d;
	312	EADD (t1, du, v, vv);
	313	s1 = v * (hij[i][11].d
	314	+ v * (hij[i][12].d
c5d5d574 OB	315	+ v * (hij[i][13].d
	316	+ v * (hij[i][14].d
	317	+ v * hij[i][15].d))));
a64d7e0e	318	ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
1728ab37 SP	319	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	320	ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
	321	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	322	ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
	323	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	324	ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
	325	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	326	ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
	327
	328	if ((z = s2 + (ss2 - ub.d * s2)) == s2 + (ss2 + ub.d * s2))
	329	return signArctan2 (y, z);
	330	return atan2Mp (x, y, pr);
	331	}
	332
	333	/* (ii) x>0, abs(x)<=abs(y): pi/2-atan(ax/ay) */
	334	if (u < inv16.d)
	335	{
	336	v = u * u;
	337	zz = u * v * (d3.d
	338	+ v * (d5.d
	339	+ v * (d7.d
	340	+ v * (d9.d
	341	+ v * (d11.d
	342	+ v * d13.d)))));
	343	ESUB (hpi.d, u, t2, cor);
	344	t3 = ((hpi1.d + cor) - du) - zz;
	345	if ((z = t2 + (t3 - u2.d)) == t2 + (t3 + u2.d))
	346	return signArctan2 (y, z);
	347
	348	MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
	349	s1 = v * (f11.d
	350	+ v * (f13.d
	351	+ v * (f15.d + v * (f17.d + v * f19.d))));
a64d7e0e	352	ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
1728ab37 SP	353	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	354	ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
	355	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	356	ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
	357	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	358	ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
	359	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	360	MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
	361	ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);
	362	SUB2 (hpi.d, hpi1.d, s1, ss1, s2, ss2, t1, t2);
	363
	364	if ((z = s2 + (ss2 - u6.d)) == s2 + (ss2 + u6.d))
	365	return signArctan2 (y, z);
	366	return atan2Mp (x, y, pr);
	367	}
	368
	369	i = (TWO52 + TWO8 * u) - TWO52;
	370	i -= 16;
	371	v = (u - cij[i][0].d) + du;
	372
	373	zz = hpi1.d - v * (cij[i][2].d
	374	+ v * (cij[i][3].d
	375	+ v * (cij[i][4].d
	376	+ v * (cij[i][5].d
	377	+ v * cij[i][6].d))));
	378	t1 = hpi.d - cij[i][1].d;
	379	if (i < 112)
	380	ua = ua1.d; /* w < 1/2 */
	381	else
	382	ua = ua2.d; /* w >= 1/2 */
	383	if ((z = t1 + (zz - ua)) == t1 + (zz + ua))
	384	return signArctan2 (y, z);
	385
	386	t1 = u - hij[i][0].d;
	387	EADD (t1, du, v, vv);
	388
	389	s1 = v * (hij[i][11].d
	390	+ v * (hij[i][12].d
	391	+ v * (hij[i][13].d
	392	+ v * (hij[i][14].d
	393	+ v * hij[i][15].d))));
	394
a64d7e0e	395	ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
1728ab37 SP	396	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	397	ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
	398	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	399	ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
	400	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	401	ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
	402	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	403	ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
	404	SUB2 (hpi.d, hpi1.d, s2, ss2, s1, ss1, t1, t2);
	405
	406	if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d))
	407	return signArctan2 (y, z);
	408	return atan2Mp (x, y, pr);
e4d82761	409	}
1728ab37 SP	410
	411	/* (iii) x<0, abs(x)< abs(y): pi/2+atan(ax/ay) */
	412	if (ax < ay)
	413	{
	414	if (u < inv16.d)
	415	{
	416	v = u * u;
	417	zz = u * v * (d3.d
	418	+ v * (d5.d
	419	+ v * (d7.d
	420	+ v * (d9.d
	421	+ v * (d11.d + v * d13.d)))));
	422	EADD (hpi.d, u, t2, cor);
	423	t3 = ((hpi1.d + cor) + du) + zz;
	424	if ((z = t2 + (t3 - u3.d)) == t2 + (t3 + u3.d))
	425	return signArctan2 (y, z);
	426
	427	MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
	428	s1 = v * (f11.d
	429	+ v * (f13.d + v * (f15.d + v * (f17.d + v * f19.d))));
a64d7e0e	430	ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
1728ab37 SP	431	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	432	ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
	433	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	434	ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
	435	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	436	ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
	437	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	438	MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
	439	ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);
	440	ADD2 (hpi.d, hpi1.d, s1, ss1, s2, ss2, t1, t2);
	441
	442	if ((z = s2 + (ss2 - u7.d)) == s2 + (ss2 + u7.d))
	443	return signArctan2 (y, z);
	444	return atan2Mp (x, y, pr);
	445	}
	446
	447	i = (TWO52 + TWO8 * u) - TWO52;
	448	i -= 16;
	449	v = (u - cij[i][0].d) + du;
	450	zz = hpi1.d + v * (cij[i][2].d
	451	+ v * (cij[i][3].d
	452	+ v * (cij[i][4].d
	453	+ v * (cij[i][5].d
	454	+ v * cij[i][6].d))));
	455	t1 = hpi.d + cij[i][1].d;
	456	if (i < 112)
	457	ua = ua1.d; /* w < 1/2 */
	458	else
	459	ua = ua2.d; /* w >= 1/2 */
	460	if ((z = t1 + (zz - ua)) == t1 + (zz + ua))
	461	return signArctan2 (y, z);
	462
	463	t1 = u - hij[i][0].d;
	464	EADD (t1, du, v, vv);
	465	s1 = v * (hij[i][11].d
	466	+ v * (hij[i][12].d
	467	+ v * (hij[i][13].d
	468	+ v * (hij[i][14].d
	469	+ v * hij[i][15].d))));
a64d7e0e	470	ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
1728ab37 SP	471	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	472	ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
	473	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	474	ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
	475	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	476	ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
	477	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	478	ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
	479	ADD2 (hpi.d, hpi1.d, s2, ss2, s1, ss1, t1, t2);
	480
	481	if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d))
	482	return signArctan2 (y, z);
	483	return atan2Mp (x, y, pr);
e4d82761 UD	484	}
e4d82761 UD	485
1728ab37 SP	486	/* (iv) x<0, abs(y)<=abs(x): pi-atan(ax/ay) */
	487	if (u < inv16.d)
	488	{
	489	v = u * u;
	490	zz = u * v * (d3.d
	491	+ v * (d5.d
	492	+ v * (d7.d
	493	+ v * (d9.d + v * (d11.d + v * d13.d)))));
	494	ESUB (opi.d, u, t2, cor);
	495	t3 = ((opi1.d + cor) - du) - zz;
	496	if ((z = t2 + (t3 - u4.d)) == t2 + (t3 + u4.d))
	497	return signArctan2 (y, z);
	498
	499	MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
	500	s1 = v * (f11.d + v * (f13.d + v * (f15.d + v * (f17.d + v * f19.d))));
a64d7e0e	501	ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
1728ab37 SP	502	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	503	ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
	504	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	505	ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
	506	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	507	ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
	508	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	509	MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
	510	ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);
	511	SUB2 (opi.d, opi1.d, s1, ss1, s2, ss2, t1, t2);
	512
	513	if ((z = s2 + (ss2 - u8.d)) == s2 + (ss2 + u8.d))
	514	return signArctan2 (y, z);
	515	return atan2Mp (x, y, pr);
e4d82761	516	}
1728ab37 SP	517
	518	i = (TWO52 + TWO8 * u) - TWO52;
	519	i -= 16;
	520	v = (u - cij[i][0].d) + du;
	521	zz = opi1.d - v * (cij[i][2].d
	522	+ v * (cij[i][3].d
	523	+ v * (cij[i][4].d
	524	+ v * (cij[i][5].d + v * cij[i][6].d))));
	525	t1 = opi.d - cij[i][1].d;
	526	if (i < 112)
	527	ua = ua1.d; /* w < 1/2 */
	528	else
	529	ua = ua2.d; /* w >= 1/2 */
	530	if ((z = t1 + (zz - ua)) == t1 + (zz + ua))
	531	return signArctan2 (y, z);
	532
	533	t1 = u - hij[i][0].d;
	534
	535	EADD (t1, du, v, vv);
	536
	537	s1 = v * (hij[i][11].d
	538	+ v * (hij[i][12].d
	539	+ v * (hij[i][13].d
	540	+ v * (hij[i][14].d + v * hij[i][15].d))));
	541
a64d7e0e	542	ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
1728ab37 SP	543	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	544	ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
	545	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	546	ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
	547	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	548	ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
	549	MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
	550	ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
	551	SUB2 (opi.d, opi1.d, s2, ss2, s1, ss1, t1, t2);
	552
	553	if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d))
	554	return signArctan2 (y, z);
	555	return atan2Mp (x, y, pr);
e4d82761	556	}
1728ab37	557
af968f62	558	#ifndef __ieee754_atan2
0ac5ae23	559	strong_alias (__ieee754_atan2, __atan2_finite)
af968f62	560	#endif
0ac5ae23	561
1728ab37	562	/* Treat the Denormalized case */
31d3cc00 UD	563	static double
31d3cc00 UD	564	SECTION
1728ab37 SP	565	normalized (double ax, double ay, double y, double z)
	566	{
	567	int p;
	568	mp_no mpx, mpy, mpz, mperr, mpz2, mpt1;
	569	p = 6;
	570	__dbl_mp (ax, &mpx, p);
	571	__dbl_mp (ay, &mpy, p);
	572	__dvd (&mpy, &mpx, &mpz, p);
	573	__dbl_mp (ue.d, &mpt1, p);
	574	__mul (&mpz, &mpt1, &mperr, p);
	575	__sub (&mpz, &mperr, &mpz2, p);
	576	__mp_dbl (&mpz2, &z, p);
	577	return signArctan2 (y, z);
e4d82761	578	}
1728ab37 SP	579
1728ab37 SP	580	/* Stage 3: Perform a multi-Precision computation */
31d3cc00 UD	581	static double
31d3cc00 UD	582	SECTION
1728ab37	583	atan2Mp (double x, double y, const int pr[])
e4d82761	584	{
1728ab37 SP	585	double z1, z2;
	586	int i, p;
	587	mp_no mpx, mpy, mpz, mpz1, mpz2, mperr, mpt1;
	588	for (i = 0; i < MM; i++)
	589	{
	590	p = pr[i];
	591	__dbl_mp (x, &mpx, p);
	592	__dbl_mp (y, &mpy, p);
	593	__mpatan2 (&mpy, &mpx, &mpz, p);
	594	__dbl_mp (ud[i].d, &mpt1, p);
	595	__mul (&mpz, &mpt1, &mperr, p);
	596	__add (&mpz, &mperr, &mpz1, p);
	597	__sub (&mpz, &mperr, &mpz2, p);
	598	__mp_dbl (&mpz1, &z1, p);
	599	__mp_dbl (&mpz2, &z2, p);
	600	if (z1 == z2)
10e1cf6b SP	601	{
	602	LIBC_PROBE (slowatan2, 4, &p, &x, &y, &z1);
	603	return z1;
	604	}
1728ab37	605	}
10e1cf6b	606	LIBC_PROBE (slowatan2_inexact, 4, &p, &x, &y, &z1);
1728ab37	607	return z1; /if impossible to do exact computing /
f7eac6eb	608	}