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, 103 insertions, 72 deletions

diff --git a/sysdeps/ieee754/dbl-64/e_remainder.c b/sysdeps/ieee754/dbl-64/e_remainder.c
index 6418081..631c2fc 100644
--- a/sysdeps/ieee754/dbl-64/e_remainder.c
+++ b/sysdeps/ieee754/dbl-64/e_remainder.c
@@ -1,80 +1,111 @@
-/* @(#)e_remainder.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_remainder.c,v 1.8 1995/05/10 20:46:05 jtc Exp $";
-#endif
-
-/* __ieee754_remainder(x,p)
- * Return :                  
- * 	returns  x REM p  =  x - [x/p]*p as if in infinite 
- * 	precise arithmetic, where [x/p] is the (infinite bit) 
- *	integer nearest x/p (in half way case choose the even one).
- * Method : 
- *	Based on fmod() return x-[x/p]chopped*p exactlp.
+ * 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.
+ *
+ * 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.
+ *
+ * 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 urem.c                                                    */
+/*                                                                        */
+/*  FUNCTION: uremainder                                                  */
+/*                                                                        */
+/* An ultimate remainder routine. Given two IEEE double machine numbers x */
+/* ,y   it computes the correctly rounded (to nearest) value of remainder */
+/* of dividing x by y.                                                    */
+/* Assumption: Machine arithmetic operations are performed in             */
+/* round to nearest mode of IEEE 754 standard.                            */
+/*                                                                        */
+/* ************************************************************************/
 
-#include "math.h"
-#include "math_private.h"
+#include "endian.h"
+#include "mydefs.h"
+#include "urem.h"
+#include "MathLib.h"
 
-#ifdef __STDC__
-static const double zero = 0.0;
-#else
-static double zero = 0.0;
-#endif
-
-
-#ifdef __STDC__
-	double __ieee754_remainder(double x, double p)
-#else
-	double __ieee754_remainder(x,p)
-	double x,p;
-#endif
+/**************************************************************************/
+/* An ultimate remainder routine. Given two IEEE double machine numbers x */
+/* ,y   it computes the correctly rounded (to nearest) value of remainder */
+/**************************************************************************/
+double __ieee754_remainder(double x, double y)
 {
-	int32_t hx,hp;
-	u_int32_t sx,lx,lp;
-	double p_half;
-
-	EXTRACT_WORDS(hx,lx,x);
-	EXTRACT_WORDS(hp,lp,p);
-	sx = hx&0x80000000;
-	hp &= 0x7fffffff;
-	hx &= 0x7fffffff;
-
-    /* purge off exception values */
-	if((hp|lp)==0) return (x*p)/(x*p); 	/* p = 0 */
-	if((hx>=0x7ff00000)||			/* x not finite */
-	  ((hp>=0x7ff00000)&&			/* p is NaN */
-	  (((hp-0x7ff00000)|lp)!=0)))
-	    return (x*p)/(x*p);
-
+  double z,d,xx,yy;
+  int4 kx,ky,m,n,nn,n1,m1,l;
+  mynumber u,t,w={0,0},v={0,0},ww={0,0},r;
+  u.x=x;
+  t.x=y;
+  kx=u.i[HIGH_HALF]&0x7fffffff; /* no sign  for x*/
+  t.i[HIGH_HALF]&=0x7fffffff;   /*no sign for y */
+  ky=t.i[HIGH_HALF];
+  /*------ |x| < 2^1024  and   2^-970 < |y| < 2^1024 ------------------*/
+  if (kx<0x7ff00000 && ky<0x7ff00000 && ky>=0x03500000) {
+    if (kx+0x00100000<ky) return x;
+    if ((kx-0x01500000)<ky) {
+      z=x/t.x;
+      v.i[HIGH_HALF]=t.i[HIGH_HALF];
+      d=(z+big.x)-big.x;
+      xx=(x-d*v.x)-d*(t.x-v.x);
+      if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x);
+      else {
+	if (ABS(xx)>0.5*t.x) return (z>d)?xx-t.x:xx+t.x;
+	else return xx;
+      }
+    }   /*    (kx<(ky+0x01500000))         */
+    else  {
+      r.x=1.0/t.x;
+      n=t.i[HIGH_HALF];
+      nn=(n&0x7ff00000)+0x01400000;
+      w.i[HIGH_HALF]=n;
+      ww.x=t.x-w.x;
+      l=(kx-nn)&0xfff00000;
+      n1=ww.i[HIGH_HALF];
+      m1=r.i[HIGH_HALF];
+      while (l>0) {
+	r.i[HIGH_HALF]=m1-l;
+	z=u.x*r.x;
+	w.i[HIGH_HALF]=n+l;
+	ww.i[HIGH_HALF]=(n1)?n1+l:n1;
+	d=(z+big.x)-big.x;
+	u.x=(u.x-d*w.x)-d*ww.x;
+	l=(u.i[HIGH_HALF]&0x7ff00000)-nn;
+      }
+      r.i[HIGH_HALF]=m1;
+      w.i[HIGH_HALF]=n;
+      ww.i[HIGH_HALF]=n1;
+      z=u.x*r.x;
+      d=(z+big.x)-big.x;
+      u.x=(u.x-d*w.x)-d*ww.x;
+      if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x);
+      else
+        if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
+        else
+        {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);}
+    }
 
-	if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);	/* now x < 2p */
-	if (((hx-hp)|(lx-lp))==0) return zero*x;
-	x  = fabs(x);
-	p  = fabs(p);
-	if (hp<0x00200000) {
-	    if(x+x>p) {
-		x-=p;
-		if(x+x>=p) x -= p;
-	    }
-	} else {
-	    p_half = 0.5*p;
-	    if(x>p_half) {
-		x-=p;
-		if(x>=p_half) x -= p;
-	    }
-	}
-	GET_HIGH_WORD(hx,x);
-	SET_HIGH_WORD(x,hx^sx);
-	return x;
+  }   /*   (kx<0x7ff00000&&ky<0x7ff00000&&ky>=0x03500000)     */
+  else {
+    if (kx<0x7ff00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
+      y=ABS(y)*t128.x;
+      z=uremainder(x,y)*t128.x;
+      z=uremainder(z,y)*tm128.x;
+      return z;
+    }
+    else { /* if x is too big */
+      if (kx>=0x7ff00000||(ky==0&&t.i[LOW_HALF]==0)||ky>0x7ff00000||
+	  (ky==0x7ff00000&&t.i[LOW_HALF]!=0))
+	return (u.i[HIGH_HALF]&0x80000000)?nNAN.x:NAN.x;
+      else return x;
+    }
+  }
 }