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, 124 insertions, 0 deletions

diff --git a/sysdeps/ieee754/dbl-64/halfulp.c b/sysdeps/ieee754/dbl-64/halfulp.c
new file mode 100644
index 0000000000..7901ec4e36
--- /dev/null
+++ b/sysdeps/ieee754/dbl-64/halfulp.c
@@ -0,0 +1,124 @@
+
+/*
+ * IBM Accurate Mathematical Library
+ * Copyright (c) International Business Machines Corp., 2001
+ *
+ * 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:halfulp.c                                                */ 
+/*                                                                      */ 
+/*  FUNCTIONS:halfulp                                                   */
+/*  FILES NEEDED: mydefs.h dla.h endian.h                               */
+/*                uroot.c                                               */
+/*                                                                      */
+/*Routine halfulp(double x, double y) computes x^y where result does    */
+/*not need rounding. If the result is closer to 0 than can be           */
+/*represented it returns 0.                                             */
+/*     In the following cases the function does not compute anything    */
+/*and returns a negative number:                                        */
+/*1. if the result needs rounding,                                      */
+/*2. if y is outside the interval [0,  2^20-1],                         */
+/*3. if x can be represented by  x=2**n for some integer n.             */
+/************************************************************************/
+
+#include "endian.h"        
+#include "mydefs.h"
+#include "dla.h"
+
+double usqrt(double x);
+
+int4 tab54[32] = {
+   262143, 11585, 1782, 511, 210, 107, 63, 42,
+       30,    22,   17,  14,  12,  10,  9,  7,
+        7,     6,    5,   5,   5,   4,  4,  4,
+        3,     3,    3,   3,   3,   3,  3,  3 };
+
+
+double halfulp(double x, double y)
+{
+  mynumber v;
+  double z,u,uu,j1,j2,j3,j4,j5;
+  int4 k,l,m,n;
+  if (y <= 0) {               /*if power is negative or zero */
+    v.x = y;
+    if (v.i[LOW_HALF] != 0) return -10.0;  
+    v.x = x;
+    if (v.i[LOW_HALF] != 0) return -10.0;   
+    if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10;   /* if x =2 ^ n */ 
+    k = ((v.i[HIGH_HALF]&0x7fffffff)>>20)-1023;         /* find this n */
+    z = (double) k;
+    return (z*y == -1075.0)?0: -10.0;
+  }
+                              /* if y > 0  */
+  v.x = y;
+    if (v.i[LOW_HALF] != 0) return -10.0;
+  
+  v.x=x;
+                              /*  case where x = 2**n for some integer n */ 
+  if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) {
+    k=(v.i[HIGH_HALF]>>20)-1023;
+    return (((double) k)*y == -1075.0)?0:-10.0;
+  }  
+  
+  v.x = y;
+  k = v.i[HIGH_HALF];
+  m = k<<12;
+  l = 0;
+  while (m) 
+    {m = m<<1; l++; }
+  n = (k&0x000fffff)|0x00100000;
+  n = n>>(20-l);                       /*   n is the odd integer of y    */
+  k = ((k>>20) -1023)-l;               /*   y = n*2**k                   */
+  if (k>5) return -10.0;
+  if (k>0) for (;k>0;k--) n *= 2;
+  if (n > 34) return -10.0;
+  k = -k;
+  if (k>5) return -10.0;
+  
+                            /*   now treat x        */
+  while (k>0) {
+    z = usqrt(x);
+    EMULV(z,z,u,uu,j1,j2,j3,j4,j5);
+    if (((u-x)+uu) != 0) break;
+    x = z;
+    k--;
+ }
+  if (k) return -10.0; 
+  
+  /* it is impossible that n == 2,  so the mantissa of x must be short  */
+  
+  v.x = x;
+  if (v.i[LOW_HALF]) return -10.0;
+  k = v.i[HIGH_HALF];
+  m = k<<12;
+  l = 0;
+  while (m) {m = m<<1; l++; }
+  m = (k&0x000fffff)|0x00100000;
+  m = m>>(20-l);                       /*   m is the odd integer of x    */
+  
+            /*   now check whether the length of m**n is at most 54 bits */
+  
+  if  (m > tab54[n-3]) return -10.0;
+  
+             /* yes, it is - now compute x**n by simple multiplications  */
+  
+  u = x;
+  for (k=1;k<n;k++) u = u*x;
+  return u;
+}
+
+