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

diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c
new file mode 100644
index 0000000000..cdd15a1602
--- /dev/null
+++ b/sysdeps/ieee754/dbl-64/mpa.c
@@ -0,0 +1,497 @@
+
+/*
+ * 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: mpa.c                                                  */
+/*                                                                      */
+/*  FUNCTIONS:                                                          */
+/*               mcr                                                    */
+/*               acr                                                    */
+/*               cr                                                     */
+/*               cpy                                                    */
+/*               cpymn                                                  */
+/*               norm                                                   */
+/*               denorm                                                 */
+/*               mp_dbl                                                 */
+/*               dbl_mp                                                 */
+/*               add_magnitudes                                         */
+/*               sub_magnitudes                                         */
+/*               add                                                    */
+/*               sub                                                    */
+/*               mul                                                    */
+/*               inv                                                    */
+/*               dvd                                                    */
+/*                                                                      */
+/* Arithmetic functions for multiple precision numbers.                 */
+/* Relative errors are bounded                                          */
+/************************************************************************/
+
+
+#include "endian.h"
+#include "mpa.h"
+#include "mpa2.h"
+/* mcr() compares the sizes of the mantissas of two multiple precision  */
+/* numbers. Mantissas are compared regardless of the signs of the       */
+/* numbers, even if x->d[0] or y->d[0] are zero. Exponents are also     */
+/* disregarded.                                                         */
+static int mcr(const mp_no *x, const mp_no *y, int p) {
+  int i;
+  for (i=1; i<=p; i++) {
+    if      (X[i] == Y[i])  continue;
+    else if (X[i] >  Y[i])  return  1;
+    else                    return -1; }
+  return 0;
+}
+
+
+
+/* acr() compares the absolute values of two multiple precision numbers */
+int acr(const mp_no *x, const mp_no *y, int p) {
+  int i;
+
+  if      (X[0] == ZERO) {
+    if    (Y[0] == ZERO) i= 0;
+    else                 i=-1;
+  }
+  else if (Y[0] == ZERO) i= 1;
+  else {
+    if      (EX >  EY)   i= 1;
+    else if (EX <  EY)   i=-1;
+    else                 i= mcr(x,y,p);
+  }
+
+  return i;
+}
+
+
+/* cr90 compares the values of two multiple precision numbers           */
+int  cr(const mp_no *x, const mp_no *y, int p) {
+  int i;
+
+  if      (X[0] > Y[0])  i= 1;
+  else if (X[0] < Y[0])  i=-1;
+  else if (X[0] < ZERO ) i= acr(y,x,p);
+  else                   i= acr(x,y,p);
+
+  return i;
+}
+
+
+/* Copy a multiple precision number. Set *y=*x. x=y is permissible.      */
+void cpy(const mp_no *x, mp_no *y, int p) {
+  int i;
+
+  EY = EX;
+  for (i=0; i <= p; i++)    Y[i] = X[i];
+
+  return;
+}
+
+
+/* Copy a multiple precision number x of precision m into a */
+/* multiple precision number y of precision n. In case n>m, */
+/* the digits of y beyond the m'th are set to zero. In case */
+/* n<m, the digits of x beyond the n'th are ignored.        */
+/* x=y is permissible.                                      */
+
+void cpymn(const mp_no *x, int m, mp_no *y, int n) {
+
+  int i,k;
+
+  EY = EX;     k=MIN(m,n);
+  for (i=0; i <= k; i++)    Y[i] = X[i];
+  for (   ; i <= n; i++)    Y[i] = ZERO;
+
+  return;
+}
+
+/* Convert a multiple precision number *x into a double precision */
+/* number *y, normalized case  (|x| >= 2**(-1022))) */
+static void norm(const mp_no *x, double *y, int p) 
+{  
+  #define R  radixi.d
+  int i,k;
+  double a,c,u,v,z[5];
+  if (p<5) {
+    if      (p==1) c = X[1];
+    else if (p==2) c = X[1] + R* X[2];
+    else if (p==3) c = X[1] + R*(X[2]  +   R* X[3]);
+    else if (p==4) c =(X[1] + R* X[2]) + R*R*(X[3] + R*X[4]);
+  }
+  else {
+    for (a=ONE, z[1]=X[1]; z[1] < TWO23; ) 
+        {a *= TWO;   z[1] *= TWO; }
+
+    for (i=2; i<5; i++) {
+      z[i] = X[i]*a;
+      u = (z[i] + CUTTER)-CUTTER;
+      if  (u > z[i])  u -= RADIX;
+      z[i] -= u;
+      z[i-1] += u*RADIXI;
+    }
+
+    u = (z[3] + TWO71) - TWO71;
+    if (u > z[3])   u -= TWO19;
+    v = z[3]-u;
+
+    if (v == TWO18) {
+      if (z[4] == ZERO) {
+        for (i=5; i <= p; i++) {
+          if (X[i] == ZERO)   continue;
+          else                {z[3] += ONE;   break; }
+        }
+      }
+      else              z[3] += ONE;
+    }
+
+    c = (z[1] + R *(z[2] + R * z[3]))/a;
+  }
+
+  c *= X[0];
+
+  for (i=1; i<EX; i++)   c *= RADIX;
+  for (i=1; i>EX; i--)   c *= RADIXI;
+
+  *y = c;
+  return;
+#undef R
+}
+
+/* Convert a multiple precision number *x into a double precision */
+/* number *y, denormalized case  (|x| < 2**(-1022))) */
+static void denorm(const mp_no *x, double *y, int p) 
+{ 
+  int i,k;
+  double a,c,u,v,z[5];
+
+#define R  radixi.d
+  if (EX<-44 || (EX==-44 && X[1]<TWO5))
+     { *y=ZERO; return; }
+
+  if      (p==1) {
+    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=ZERO;  z[3]=ZERO;  k=3;}
+    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=ZERO;  k=2;}
+    else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
+  }
+  else if (p==2) {
+    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  z[3]=ZERO;  k=3;}
+    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=X[2];  k=2;}
+    else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
+  }
+  else {
+    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  k=3;}
+    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  k=2;}
+    else              {z[1]=     TWO10;  z[2]=ZERO;  k=1;}
+    z[3] = X[k];
+  }
+
+  u = (z[3] + TWO57) - TWO57;
+  if  (u > z[3])   u -= TWO5;
+
+  if (u==z[3]) {
+    for (i=k+1; i <= p; i++) {
+      if (X[i] == ZERO)   continue;
+      else {z[3] += ONE;   break; }
+    }
+  }
+
+  c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10);
+
+  *y = c*TWOM1032;
+  return;
+
+#undef R
+}
+
+/* Convert a multiple precision number *x into a double precision number *y. */
+/* The result is correctly rounded to the nearest/even. *x is left unchanged */
+
+void mp_dbl(const mp_no *x, double *y, int p) {
+
+  int i,k;
+  double a,c,u,v,z[5];
+
+  if (X[0] == ZERO)  {*y = ZERO;  return; }
+
+  if      (EX> -42)                 norm(x,y,p);
+  else if (EX==-42 && X[1]>=TWO10)  norm(x,y,p);
+  else                              denorm(x,y,p);
+}
+
+
+/* dbl_mp() converts a double precision number x into a multiple precision  */
+/* number *y. If the precision p is too small the result is truncated. x is */
+/* left unchanged.                                                          */
+
+void dbl_mp(double x, mp_no *y, int p) {
+
+  int i,n;
+  double u;
+
+  /* Sign */
+  if      (x == ZERO)  {Y[0] = ZERO;  return; }
+  else if (x >  ZERO)   Y[0] = ONE;
+  else                 {Y[0] = MONE;  x=-x;   }
+
+  /* Exponent */
+  for (EY=ONE; x >= RADIX; EY += ONE)   x *= RADIXI;
+  for (      ; x <  ONE;   EY -= ONE)   x *= RADIX;
+
+  /* Digits */
+  n=MIN(p,4);
+  for (i=1; i<=n; i++) {
+    u = (x + TWO52) - TWO52;
+    if (u>x)   u -= ONE;
+    Y[i] = u;     x -= u;    x *= RADIX; }
+  for (   ; i<=p; i++)     Y[i] = ZERO;
+  return;
+}
+
+
+/*  add_magnitudes() adds the magnitudes of *x & *y assuming that           */
+/*  abs(*x) >= abs(*y) > 0.                                                 */
+/* The sign of the sum *z is undefined. x&y may overlap but not x&z or y&z. */
+/* No guard digit is used. The result equals the exact sum, truncated.      */
+/* *x & *y are left unchanged.                                              */
+
+static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+
+  int i,j,k;
+
+  EZ = EX;
+
+  i=p;    j=p+ EY - EX;    k=p+1;
+
+  if (j<1)
+     {cpy(x,z,p);  return; }
+  else   Z[k] = ZERO;
+
+  for (; j>0; i--,j--) {
+    Z[k] += X[i] + Y[j];
+    if (Z[k] >= RADIX) {
+      Z[k]  -= RADIX;
+      Z[--k] = ONE; }
+    else
+      Z[--k] = ZERO;
+  }
+
+  for (; i>0; i--) {
+    Z[k] += X[i];
+    if (Z[k] >= RADIX) {
+      Z[k]  -= RADIX;
+      Z[--k] = ONE; }
+    else
+      Z[--k] = ZERO;
+  }
+
+  if (Z[1] == ZERO) {
+    for (i=1; i<=p; i++)    Z[i] = Z[i+1]; }
+  else   EZ += ONE;
+}
+
+
+/*  sub_magnitudes() subtracts the magnitudes of *x & *y assuming that      */
+/*  abs(*x) > abs(*y) > 0.                                                  */
+/* The sign of the difference *z is undefined. x&y may overlap but not x&z  */
+/* or y&z. One guard digit is used. The error is less than one ulp.         */
+/* *x & *y are left unchanged.                                              */
+
+static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+
+  int i,j,k;
+
+  EZ = EX;
+
+  if (EX == EY) {
+    i=j=k=p;
+    Z[k] = Z[k+1] = ZERO; }
+  else {
+    j= EX - EY;
+    if (j > p)  {cpy(x,z,p);  return; }
+    else {
+      i=p;   j=p+1-j;   k=p;
+      if (Y[j] > ZERO) {
+        Z[k+1] = RADIX - Y[j--];
+        Z[k]   = MONE; }
+      else {
+        Z[k+1] = ZERO;
+        Z[k]   = ZERO;   j--;}
+    }
+  }
+
+  for (; j>0; i--,j--) {
+    Z[k] += (X[i] - Y[j]);
+    if (Z[k] < ZERO) {
+      Z[k]  += RADIX;
+      Z[--k] = MONE; }
+    else
+      Z[--k] = ZERO;
+  }
+
+  for (; i>0; i--) {
+    Z[k] += X[i];
+    if (Z[k] < ZERO) {
+      Z[k]  += RADIX;
+      Z[--k] = MONE; }
+    else
+      Z[--k] = ZERO;
+  }
+
+  for (i=1; Z[i] == ZERO; i++) ;
+  EZ = EZ - i + 1;
+  for (k=1; i <= p+1; )
+    Z[k++] = Z[i++];
+  for (; k <= p; )
+    Z[k++] = ZERO;
+
+  return;
+}
+
+
+/* Add two multiple precision numbers. Set *z = *x + *y. x&y may overlap  */
+/* but not x&z or y&z. One guard digit is used. The error is less than    */
+/* one ulp. *x & *y are left unchanged.                                   */
+
+void add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+
+  int n;
+
+  if      (X[0] == ZERO)     {cpy(y,z,p);  return; }
+  else if (Y[0] == ZERO)     {cpy(x,z,p);  return; }
+
+  if (X[0] == Y[0])   {
+    if (acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] = X[0]; }
+    else                     {add_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
+  }
+  else                       {
+    if ((n=acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] = X[0]; }
+    else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
+    else                      Z[0] = ZERO;
+  }
+  return;
+}
+
+
+/* Subtract two multiple precision numbers. *z is set to *x - *y. x&y may */
+/* overlap but not x&z or y&z. One guard digit is used. The error is      */
+/* less than one ulp. *x & *y are left unchanged.                         */
+
+void sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+
+  int n;
+
+  if      (X[0] == ZERO)     {cpy(y,z,p);  Z[0] = -Z[0];  return; }
+  else if (Y[0] == ZERO)     {cpy(x,z,p);                 return; }
+
+  if (X[0] != Y[0])    {
+    if (acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
+    else                     {add_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
+  }
+  else                       {
+    if ((n=acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
+    else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
+    else                      Z[0] = ZERO;
+  }
+  return;
+}
+
+
+/* Multiply two multiple precision numbers. *z is set to *x * *y. x&y      */
+/* may overlap but not x&z or y&z. In case p=1,2,3 the exact result is     */
+/* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp.   */
+/* *x & *y are left unchanged.                                             */
+
+void mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+
+  int i, i1, i2, j, k, k2;
+  double u;
+
+                      /* Is z=0? */
+  if (X[0]*Y[0]==ZERO)
+     { Z[0]=ZERO;  return; }
+
+                       /* Multiply, add and carry */
+  k2 = (p<3) ? p+p : p+3;
+  Z[k2]=ZERO;
+  for (k=k2; k>1; ) {
+    if (k > p)  {i1=k-p; i2=p+1; }
+    else        {i1=1;   i2=k;   }
+    for (i=i1,j=i2-1; i<i2; i++,j--)  Z[k] += X[i]*Y[j];
+
+    u = (Z[k] + CUTTER)-CUTTER;
+    if  (u > Z[k])  u -= RADIX;
+    Z[k]  -= u;
+    Z[--k] = u*RADIXI;
+  }
+
+                 /* Is there a carry beyond the most significant digit? */
+  if (Z[1] == ZERO) {
+    for (i=1; i<=p; i++)  Z[i]=Z[i+1];
+    EZ = EX + EY - 1; }
+  else
+    EZ = EX + EY;
+
+  Z[0] = X[0] * Y[0];
+  return;
+}
+
+
+/* Invert a multiple precision number. Set *y = 1 / *x.                     */
+/* Relative error bound = 1.001*r**(1-p) for p=2, 1.063*r**(1-p) for p=3,   */
+/* 2.001*r**(1-p) for p>3.                                                  */
+/* *x=0 is not permissible. *x is left unchanged.                           */
+
+void inv(const mp_no *x, mp_no *y, int p) {
+  int i,l;
+  double t;
+  mp_no z,w;
+  static const int np1[] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,3,
+                            4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};
+  const mp_no mptwo = {1,1.0,2.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
+                         0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
+                         0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
+                         0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,};
+
+  cpy(x,&z,p);  z.e=0;  mp_dbl(&z,&t,p);
+  t=ONE/t;   dbl_mp(t,y,p);    EY -= EX;
+
+  for (i=0; i<np1[p]; i++) {
+    cpy(y,&w,p);
+    mul(x,&w,y,p);
+    sub(&mptwo,y,&z,p);
+    mul(&w,&z,y,p);
+  }
+  return;
+}
+
+
+/* Divide one multiple precision number by another.Set *z = *x / *y. *x & *y */
+/* are left unchanged. x&y may overlap but not x&z or y&z.                   */
+/* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3     */
+/* and 3.001*r**(1-p) for p>3. *y=0 is not permissible.                      */
+
+void dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+
+  mp_no w;
+
+  if (X[0] == ZERO)    Z[0] = ZERO;
+  else                {inv(y,&w,p);   mul(x,&w,z,p);}
+  return;
+}
+