Update.

2001-03-11  Ulrich Drepper  <drepper@redhat.com>

	Last-bit accurate math library implementation by IBM Haifa.
	Contributed by Abraham Ziv <ziv@il.ibm.com>, Moshe Olshansky
	<olshansk@il.ibm.com>, Ealan Henis <ealan@il.ibm.com>, and
	Anna Reitman <reitman@il.ibm.com>.
	* math/Makefile (dbl-only-routines): New variable.
	(libm-routines): Add $(dbl-only-routines).
	* sysdeps/ieee754/dbl-64/e_acos.c: Empty, definition is in e_asin.c.
	* sysdeps/ieee754/dbl-64/e_asin.c: Replaced with accurate asin
	implementation.
	* sysdeps/ieee754/dbl-64/e_atan2.c: Replaced with accurate atan2
	implementation.
	* sysdeps/ieee754/dbl-64/e_exp.c: Replaced with accurate exp
	implementation.
	* sysdeps/ieee754/dbl-64/e_lgamma_r.c: Don't use __kernel_sin and
	__kernel_cos.
	* sysdeps/ieee754/dbl-64/e_log.c: Replaced with accurate log
	implementation.
	* sysdeps/ieee754/dbl-64/e_remainder.c: Replaced with accurate
	remainder implementation.
	* sysdeps/ieee754/dbl-64/e_pow.c: Replaced with accurate pow
	implementation.
	* sysdeps/ieee754/dbl-64/e_sqrt.c: Replaced with accurate sqrt
	implementation.
	* sysdeps/ieee754/dbl-64/k_cos.c: Empty, definition is in s_sin.c.
	* sysdeps/ieee754/dbl-64/k_sin.c: Empty, definition is in s_sin.c.
	* sysdeps/ieee754/dbl-64/s_atan.c: Replaced with accurate atan
	implementation.
	* sysdeps/ieee754/dbl-64/s_cos.c: Empty, definition is in s_sin.c.
	* sysdeps/ieee754/dbl-64/s_sin.c: Replaced with accurate sin/cos
	implementation.
	* sysdeps/ieee754/dbl-64/s_sincos.c: Rewritten to not use __kernel_sin
	and __kernel_cos.
	* sysdeps/ieee754/dbl-64/s_tan.c: Replaced with accurate tan
	implementation.
	* sysdeps/ieee754/dbl-64/Dist: Add new non-code files.
	* sysdeps/ieee754/dbl-64/MathLib.h: New file.
	* sysdeps/ieee754/dbl-64/asincos.tbl: New file.
	* sysdeps/ieee754/dbl-64/atnat.h: New file.
	* sysdeps/ieee754/dbl-64/atnat2.h: New file.
	* sysdeps/ieee754/dbl-64/branred.c: New file.
	* sysdeps/ieee754/dbl-64/branred.h: New file.
	* sysdeps/ieee754/dbl-64/dla.h: New file.
	* sysdeps/ieee754/dbl-64/doasin.c: New file.
	* sysdeps/ieee754/dbl-64/doasin.h: New file.
	* sysdeps/ieee754/dbl-64/dosincos.c: New file.
	* sysdeps/ieee754/dbl-64/dosincos.h: New file.
	* sysdeps/ieee754/dbl-64/endian.h: New file.
	* sysdeps/ieee754/dbl-64/halfulp.c: New file.
	* sysdeps/ieee754/dbl-64/mpa.c: New file.
	* sysdeps/ieee754/dbl-64/mpa.h: New file.
	* sysdeps/ieee754/dbl-64/mpa2.h: New file.
	* sysdeps/ieee754/dbl-64/mpatan.c: New file.
	* sysdeps/ieee754/dbl-64/mpatan.h: New file.
	* sysdeps/ieee754/dbl-64/mpatan2.c: New file.
	* sysdeps/ieee754/dbl-64/mpexp.c: New file.
	* sysdeps/ieee754/dbl-64/mpexp.h: New file.
	* sysdeps/ieee754/dbl-64/mplog.c: New file.
	* sysdeps/ieee754/dbl-64/mplog.h: New file.
	* sysdeps/ieee754/dbl-64/mpsqrt.c: New file.
	* sysdeps/ieee754/dbl-64/mpsqrt.h: New file.
	* sysdeps/ieee754/dbl-64/mptan.c: New file.
	* sysdeps/ieee754/dbl-64/mydefs.h: New file.
	* sysdeps/ieee754/dbl-64/powtwo.tbl: New file.
	* sysdeps/ieee754/dbl-64/root.tbl: New file.
	* sysdeps/ieee754/dbl-64/sincos.tbl: New file.
	* sysdeps/ieee754/dbl-64/sincos32.c: New file.
	* sysdeps/ieee754/dbl-64/sincos32.h: New file.
	* sysdeps/ieee754/dbl-64/slowexp.c: New file.
	* sysdeps/ieee754/dbl-64/slowpow.c: New file.
	* sysdeps/ieee754/dbl-64/uasncs.h: New file.
	* sysdeps/ieee754/dbl-64/uatan.tbl: New file.
	* sysdeps/ieee754/dbl-64/uexp.h: New file.
	* sysdeps/ieee754/dbl-64/uexp.tbl: New file.
	* sysdeps/ieee754/dbl-64/ulog.h: New file.
	* sysdeps/ieee754/dbl-64/ulog.tbl: New file.
	* sysdeps/ieee754/dbl-64/upow.h: New file.
	* sysdeps/ieee754/dbl-64/upow.tbl: New file.
	* sysdeps/ieee754/dbl-64/urem.h: New file.
	* sysdeps/ieee754/dbl-64/uroot.h: New file.
	* sysdeps/ieee754/dbl-64/usncs.h: New file.
	* sysdeps/ieee754/dbl-64/utan.h: New file.
	* sysdeps/ieee754/dbl-64/utan.tbl: New file.
	* sysdeps/i386/fpu/branred.c: New file.
	* sysdeps/i386/fpu/doasin.c: New file.
	* sysdeps/i386/fpu/dosincos.c: New file.
	* sysdeps/i386/fpu/halfulp.c: New file.
	* sysdeps/i386/fpu/mpa.c: New file.
	* sysdeps/i386/fpu/mpatan.c: New file.
	* sysdeps/i386/fpu/mpatan2.c: New file.
	* sysdeps/i386/fpu/mpexp.c: New file.
	* sysdeps/i386/fpu/mplog.c: New file.
	* sysdeps/i386/fpu/mpsqrt.c: New file.
	* sysdeps/i386/fpu/mptan.c: New file.
	* sysdeps/i386/fpu/sincos32.c: New file.
	* sysdeps/i386/fpu/slowexp.c: New file.
	* sysdeps/i386/fpu/slowpow.c: New file.
	* sysdeps/ia64/fpu/branred.c: New file.
	* sysdeps/ia64/fpu/doasin.c: New file.
	* sysdeps/ia64/fpu/dosincos.c: New file.
	* sysdeps/ia64/fpu/halfulp.c: New file.
	* sysdeps/ia64/fpu/mpa.c: New file.
	* sysdeps/ia64/fpu/mpatan.c: New file.
	* sysdeps/ia64/fpu/mpatan2.c: New file.
	* sysdeps/ia64/fpu/mpexp.c: New file.
	* sysdeps/ia64/fpu/mplog.c: New file.
	* sysdeps/ia64/fpu/mpsqrt.c: New file.
	* sysdeps/ia64/fpu/mptan.c: New file.
	* sysdeps/ia64/fpu/sincos32.c: New file.
	* sysdeps/ia64/fpu/slowexp.c: New file.
	* sysdeps/ia64/fpu/slowpow.c: New file.
	* sysdeps/m68k/fpu/branred.c: New file.
	* sysdeps/m68k/fpu/doasin.c: New file.
	* sysdeps/m68k/fpu/dosincos.c: New file.
	* sysdeps/m68k/fpu/halfulp.c: New file.
	* sysdeps/m68k/fpu/mpa.c: New file.
	* sysdeps/m68k/fpu/mpatan.c: New file.
	* sysdeps/m68k/fpu/mpatan2.c: New file.
	* sysdeps/m68k/fpu/mpexp.c: New file.
	* sysdeps/m68k/fpu/mplog.c: New file.
	* sysdeps/m68k/fpu/mpsqrt.c: New file.
	* sysdeps/m68k/fpu/mptan.c: New file.
	* sysdeps/m68k/fpu/sincos32.c: New file.
	* sysdeps/m68k/fpu/slowexp.c: New file.
	* sysdeps/m68k/fpu/slowpow.c: New file.

	* iconvdata/gconv-modules: Add a number of alias, mostly for IBM
	codepages.

1 files changed, 388 insertions, 123 deletions

diff --git a/sysdeps/ieee754/dbl-64/e_atan2.c b/sysdeps/ieee754/dbl-64/e_atan2.c
index ae7d759a9f..11342d87d3 100644
--- a/sysdeps/ieee754/dbl-64/e_atan2.c
+++ b/sysdeps/ieee754/dbl-64/e_atan2.c
@@ -1,130 +1,395 @@
-/* @(#)e_atan2.c 5.1 93/09/24 */
 /*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * IBM Accurate Mathematical Library
+ * Copyright (c) International Business Machines Corp., 2001
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
- * is preserved.
- * ====================================================
- */
-
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_atan2.c,v 1.8 1995/05/10 20:44:51 jtc Exp $";
-#endif
-
-/* __ieee754_atan2(y,x)
- * Method :
- *	1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
- *	2. Reduce x to positive by (if x and y are unexceptional): 
- *		ARG (x+iy) = arctan(y/x)   	   ... if x > 0,
- *		ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
- *
- * Special cases:
+ * 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 of the License, or
+ * (at your option) any later version.
  *
- *	ATAN2((anything), NaN ) is NaN;
- *	ATAN2(NAN , (anything) ) is NaN;
- *	ATAN2(+-0, +(anything but NaN)) is +-0  ;
- *	ATAN2(+-0, -(anything but NaN)) is +-pi ;
- *	ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
- *	ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
- *	ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
- *	ATAN2(+-INF,+INF ) is +-pi/4 ;
- *	ATAN2(+-INF,-INF ) is +-3pi/4;
- *	ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
+ * 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 General Public License for more details.
  *
- * Constants:
- * The hexadecimal values are the intended ones for the following 
- * constants. The decimal values may be used, provided that the 
- * compiler will convert from decimal to binary accurately enough 
- * to produce the hexadecimal values shown.
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  */
+/************************************************************************/
+/*  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"
+/************************************************************************/
+/* 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[]);
+static double signArctan2(double ,double);
+static double normalized(double ,double,double ,double);
+void __mpatan2(mp_no *,mp_no *,mp_no *,int);
+
+double __ieee754_atan2(double y,double x) {
+
+  int i,de,ux,dx,uy,dy,p;
+  static const int pr[MM]={6,8,10,20,32};
+  double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t4,t5,t6,t7,t8,
+         z,zz,z1,z2,cor,s1,ss1,s2,ss2;
+  number num;
+  mp_no mperr,mpt1,mpx,mpy,mpz,mpz1,mpz2;
+
+  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 ZERO;
+      else                               return opi.d; } }
+  else if (uy==0x80000000) {
+    if    (dy==0x00000000) {
+      if  ((ux&0x80000000)==0x00000000)  return MZERO;
+      else                               return mopi.d;} }
+
+  /* x=+-0 */
+  if (x==ZERO) {
+    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 ZERO;
+        else                                 return MZERO; }
+    }
+  }
+  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<ZERO) ? -x : x;    ay = (y<ZERO) ? -y : y;
+  de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
+  if      (de>=ep)  { return ((y>ZERO) ? hpi.d : mhpi.d); }
+  else if (de<=em)  {
+    if    (x>ZERO)  {
+      if  ((z=ay/ax)<TWOM1022)  return normalized(ax,ay,y,z);
+      else                      return signArctan2(y,z); }
+    else            { return ((y>ZERO) ? 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; }
+
+  /* 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; }
 
-#include "math.h"
-#include "math_private.h"
-
-#ifdef __STDC__
-static const double 
-#else
-static double 
-#endif
-tiny  = 1.0e-300,
-zero  = 0.0,
-pi_o_4  = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
-pi_o_2  = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
-pi      = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
-pi_lo   = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
-
-#ifdef __STDC__
-	double __ieee754_atan2(double y, double x)
-#else
-	double __ieee754_atan2(y,x)
-	double  y,x;
-#endif
-{  
-	double z;
-	int32_t k,m,hx,hy,ix,iy;
-	u_int32_t lx,ly;
-
-	EXTRACT_WORDS(hx,lx,x);
-	ix = hx&0x7fffffff;
-	EXTRACT_WORDS(hy,ly,y);
-	iy = hy&0x7fffffff;
-	if(((ix|((lx|-lx)>>31))>0x7ff00000)||
-	   ((iy|((ly|-ly)>>31))>0x7ff00000))	/* x or y is NaN */
-	   return x+y;
-	if(((hx-0x3ff00000)|lx)==0) return __atan(y);   /* x=1.0 */
-	m = ((hy>>31)&1)|((hx>>30)&2);	/* 2*sign(x)+sign(y) */
-
-    /* when y = 0 */
-	if((iy|ly)==0) {
-	    switch(m) {
-		case 0: 
-		case 1: return y; 	/* atan(+-0,+anything)=+-0 */
-		case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
-		case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
-	    }
-	}
-    /* when x = 0 */
-	if((ix|lx)==0) return (hy<0)?  -pi_o_2-tiny: pi_o_2+tiny;
-	    
-    /* when x is INF */
-	if(ix==0x7ff00000) {
-	    if(iy==0x7ff00000) {
-		switch(m) {
-		    case 0: return  pi_o_4+tiny;/* atan(+INF,+INF) */
-		    case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
-		    case 2: return  3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
-		    case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
-		}
-	    } else {
-		switch(m) {
-		    case 0: return  zero  ;	/* atan(+...,+INF) */
-		    case 1: return -zero  ;	/* atan(-...,+INF) */
-		    case 2: return  pi+tiny  ;	/* atan(+...,-INF) */
-		    case 3: return -pi-tiny  ;	/* atan(-...,-INF) */
-		}
-	    }
-	}
-    /* when y is INF */
-	if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
-
-    /* compute y/x */
-	k = (iy-ix)>>20;
-	if(k > 60) z=pi_o_2+0.5*pi_lo; 	/* |y/x| >  2**60 */
-	else if(hx<0&&k<-60) z=0.0; 	/* |y|/x < -2**60 */
-	else z=__atan(fabs(y/x));	/* safe to do y/x */
-	switch (m) {
-	    case 0: return       z  ;	/* atan(+,+) */
-	    case 1: {
-	    	      u_int32_t zh;
-		      GET_HIGH_WORD(zh,z);
-		      SET_HIGH_WORD(z,zh ^ 0x80000000);
-		    }
-		    return       z  ;	/* atan(-,+) */
-	    case 2: return  pi-(z-pi_lo);/* atan(+,-) */
-	    default: /* case 3 */
-	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
-	}
+  if (x>ZERO) {
+
+    /* (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,ZERO,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);
+      }
+      else {
+        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,ZERO,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) */
+    else {
+      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,ZERO,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);
+      }
+      else {
+        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,ZERO,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);
+      }
+    }
+  }
+  else {
+
+    /* (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,ZERO,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);
+      }
+      else {
+        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,ZERO,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) */
+    else {
+      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,ZERO,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);
+      }
+      else {
+        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,ZERO,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);
+      }
+    }
+  }
+}
+  /* Treat the Denormalized case */
+static double  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);
+}
+  /* Fix the sign and return after stage 1 or stage 2 */
+static double signArctan2(double y,double z)
+{
+  return ((y<ZERO) ? -z : z);
+}
+  /* Stage 3: Perform a multi-Precision computation */
+static double  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)   return z1;
+  }
+  return z1; /*if unpossible to do exact computing */
 }