Partial backports of:
=====================
commit c5d5d574cbfa96d0f6c1db24d1e072c472627e41
Author: Ondřej Bílka <neleai@seznam.cz>
Date: Thu Oct 17 16:03:24 2013 +0200
Format floating routines.
commit da08f647d58d674db08cdb3e61c8826c89470e2e
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Fri Dec 21 09:15:10 2012 +0530
Move more constants into static variables
Code cleanup.
commit f93a8d15699ee699282465dc1e03e033f3fabb52
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Wed Jan 16 16:06:48 2013 +0530
Consolidate constant defines into mpa.h
commit caa99d06e7f1403887294442af520b0f8c6f3de0
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Fri Jan 18 11:18:13 2013 +0530
Simplify calculation of 2^-m in __mpexp
commit 107a5bf085f5c4ef8c28266a34d476724cfc3475
Author: Joseph Myers <joseph@codesourcery.com>
Date: Tue Nov 18 15:40:56 2014 +0000
Fix libm mpone, mptwo namespace (bug 17616).
To provided __mptwo for __inv.
Full backports of the following:
================================
commit 44e0d4c20ce5bf3825897e5d4b7caae94016214d
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Wed Jan 2 11:44:13 2013 +0530
Split mantissa calculation loop and add branch prediction
commit f8af25d218202ff2f5d167b8e44e4b79f91d147f
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Fri Jan 4 15:09:33 2013 +0530
Remove commented declarations
commit a9e48ab40e230c7fe34e4892bec8af4f3f975a20
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Fri Jan 4 15:42:09 2013 +0530
Clean up comment for MP_NO
commit fffb407f4668b40b3df1eb8ee3ae3bc64ee79e20
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Fri Jan 4 22:52:12 2013 +0530
Remove unused __cr and __cpymn
commit 950c99ca9094e7dc6394e90395f51e12093393aa
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Wed Jan 9 19:07:15 2013 +0530
Update comments in mpa.c
Fixed comment style and clearer wording in some cases.
commit 1066a53440d2744566e97c59bcd0d422186b3e90
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Mon Jan 14 21:31:25 2013 +0530
Fix code formatting in mpa.c
This includes the overridden mpa.c in power4.
commit 2a91b5735ac1bc65ce5c2a3646d75ba7208e26e9
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Mon Jan 14 21:36:58 2013 +0530
Minor tweak to mp multiplication
Add a local variable to remove extra copies to/from memory in the Z
array.
ommit 45f058844c33f670475bd02f266942746bcb332b
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Tue Feb 26 21:28:16 2013 +0530
Another tweak to the multiplication algorithm
Reduce the formula to calculate mantissa so that we reduce the net
number of multiplications performed.
commit bab8a695ee79a5a6e9b2b699938952b006fcbbec
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Thu Feb 21 14:29:18 2013 +0530
Fix whitespace differences between generic and powerpc mpa.c
commit 2d0e0f29f85036d1189231cb7c1f19f27ba14a89
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Fri Feb 15 23:56:20 2013 +0530
Fix determination of lower precision in __mul
commit 909279a5cfa938c989c9b01c8f48a0276291ec45
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Wed Feb 13 14:16:23 2013 +0530
Optimized mp multiplication
Don't bother multiplying zeroes since that only wastes cycles.
commit bdf028142eb77d6ae59500db875068fa5d7b059d
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Wed Feb 13 13:55:29 2013 +0530
Clean up add_magnitudes and sub_magnitudes
commit d6752ccd696c71d23cd3df8fb9cc60b61c32e65a
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Thu Feb 14 10:31:09 2013 +0530
New __sqr function as a faster special case of __mul
commit 4709fe7602b56e9f6ee1ab6afb4067409a784f29
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Sat Feb 16 00:09:29 2013 +0530
Use intermediate variable to compute exponent in __mul
commit 20cd7fb3ae63795ae7c9a464abf5ed19b364ade0
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Wed Feb 20 18:56:20 2013 +0530
Copy comment about inner loop from powerpc mpa.c to the default one
commit e69804d14e43f14c3c65dc570afdbfb822c9838b
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Mon Feb 25 16:43:02 2013 +0530
Use long wherever possible in mpa.c
Using long throughout like powerpc does is beneficial since it reduces
the need to switch to 32-bit instructions. It gives a very minor
performance improvement.
commit 82a9811d29c00161c7c8ea7f70e2cc30988e192e
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Thu Mar 7 12:23:29 2013 +0530
Use generic mpa.c code for everything except __mul and __sqr
commit 6f2e90e78f151bab153c2b38492505ae2012db06
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Tue Mar 26 19:28:50 2013 +0530
Make mantissa type of mp_no configurable
The mantissa of mp_no is intended to take only integral values. This
is a relatively good choice for powerpc due to its 4 fpus, but not for
other architectures, which suffer due to this choice. This change
makes the default mantissa a long integer and allows powerpc to
override it. Additionally, some operations have been optimized for
integer manipulation, resulting in a significant improvement in
performance.
commit 5739f705eed5cf58e7b439e5983542e06d7fc2da
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Tue Mar 26 20:24:04 2013 +0530
Use integral constants
The compiler is smart enough to convert those into double for powerpc,
but if we put them as doubles, it adds overhead by performing those
operations in floating point mode.
commit 89f3b6e18c6e7833438789746fcfc2e7189f7cac
Author: Joseph Myers <joseph@codesourcery.com>
Date: Thu May 21 23:05:45 2015 +0000
Fix sysdeps/ieee754/dbl-64/mpa.c for -Wuninitialized.
If you remove the "override CFLAGS += -Wno-uninitialized" in
math/Makefile, one of the errors you get is:
../sysdeps/ieee754/dbl-64/mpa.c: In function '__mp_dbl.part.0':
../sysdeps/ieee754/dbl-64/mpa.c:183:5: error: 'c' may be used uninitialized in this function [-Werror=maybe-uninitialized]
c *= X[0];
The problem is that the p < 5 case initializes c if p is 1, 2, 3 or 4
but not otherwise, and in fact p is positive for all calls to this
function so the uninitialized case can't actually occur. This patch
replaces the "if (p == 4)" last case with a comment so the compiler
can see that all paths do initialize c.
Tested for x86_64.
* sysdeps/ieee754/dbl-64/mpa.c (norm): Remove if condition on
(p == 4) case.
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpa.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation
+ * Copyright (C) 2001-2017 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
@@ -22,9 +22,7 @@
/* FUNCTIONS: */
/* mcr */
/* acr */
-/* cr */
/* cpy */
-/* cpymn */
/* norm */
/* denorm */
/* mp_dbl */
@@ -44,479 +42,868 @@
#include "endian.h"
#include "mpa.h"
-#include "mpa2.h"
-#include <sys/param.h> /* For MIN() */
+#include <sys/param.h>
+#include <alloca.h>
#ifndef SECTION
# define SECTION
#endif
+#ifndef NO__CONST
+/* TODO: With only a partial backport of the constant cleanup from
+ upstream we limit the definition of these two constants to
+ this file. */
+static const mp_no __mpone = { 1, { 1.0, 1.0 } };
+static const mp_no __mptwo = { 1, { 1.0, 2.0 } };
+#endif
+
#ifndef NO___ACR
-/* 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. */
+/* Compare mantissa of two multiple precision numbers regardless of the sign
+ and exponent of the numbers. */
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; }
+mcr (const mp_no *x, const mp_no *y, int p)
+{
+ long i;
+ long p2 = p;
+ for (i = 1; i <= p2; 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 */
+/* Compare 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;
-}
-#endif
-
+__acr (const mp_no *x, const mp_no *y, int p)
+{
+ long i;
-#if 0
-/* cr() compares the values of two multiple precision numbers */
-static 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);
+ if (X[0] == 0)
+ {
+ if (Y[0] == 0)
+ i = 0;
+ else
+ i = -1;
+ }
+ else if (Y[0] == 0)
+ i = 1;
+ else
+ {
+ if (EX > EY)
+ i = 1;
+ else if (EX < EY)
+ i = -1;
+ else
+ i = mcr (x, y, p);
+ }
return i;
}
#endif
-
#ifndef NO___CPY
-/* Copy a multiple precision number. Set *y=*x. x=y is permissible. */
-void __cpy(const mp_no *x, mp_no *y, int p) {
+/* Copy multiple precision number X into Y. They could be the same
+ number. */
+void
+__cpy (const mp_no *x, mp_no *y, int p)
+{
+ long i;
+
EY = EX;
- for (int i=0; i <= p; i++) Y[i] = X[i];
+ for (i = 0; i <= p; i++)
+ Y[i] = X[i];
}
#endif
+#ifndef NO___MP_DBL
+/* Convert a multiple precision number *X into a double precision
+ number *Y, normalized case (|x| >= 2**(-1022))). X has precision
+ P, which is positive. */
+static void
+norm (const mp_no *x, double *y, int p)
+{
+# define R RADIXI
+ long i;
+ double c;
+ mantissa_t a, 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 /* p == 4. */
+ c = (X[1] + R * X[2]) + R * R * (X[3] + R * X[4]);
+ }
+ else
+ {
+ for (a = 1, z[1] = X[1]; z[1] < TWO23; )
+ {
+ a *= 2;
+ z[1] *= 2;
+ }
-#if 0
-/* 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. */
-
-static 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;
-}
-#endif
+ for (i = 2; i < 5; i++)
+ {
+ mantissa_store_t d, r;
+ d = X[i] * (mantissa_store_t) a;
+ DIV_RADIX (d, r);
+ z[i] = r;
+ z[i - 1] += d;
+ }
+ u = ALIGN_DOWN_TO (z[3], TWO19);
+ v = z[3] - u;
-#ifndef NO___MP_DBL
-/* 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;
-#if 0
- int k;
-#endif
- 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; }
+ if (v == TWO18)
+ {
+ if (z[4] == 0)
+ {
+ for (i = 5; i <= p; i++)
+ {
+ if (X[i] == 0)
+ continue;
+ else
+ {
+ z[3] += 1;
+ break;
+ }
+ }
+ }
+ else
+ z[3] += 1;
}
- }
- else z[3] += ONE;
- }
- c = (z[1] + R *(z[2] + R * z[3]))/a;
- }
+ 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;
+ for (i = 1; i < EX; i++)
+ c *= RADIX;
+ for (i = 1; i > EX; i--)
+ c *= RADIXI;
*y = c;
-#undef R
+# 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)
+/* Convert a multiple precision number *X into a double precision
+ number *Y, Denormal case (|x| < 2**(-1022))). */
+static void
+denorm (const mp_no *x, double *y, int p)
{
- int i,k;
- double c,u,z[5];
-#if 0
- double a,v;
-#endif
+ long i, k;
+ long p2 = p;
+ double c;
+ mantissa_t u, z[5];
+
+# define R RADIXI
+ if (EX < -44 || (EX == -44 && X[1] < TWO5))
+ {
+ *y = 0;
+ return;
+ }
-#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;
-#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) {
-#if 0
- int i,k;
- double a,c,u,v,z[5];
-#endif
+ if (p2 == 1)
+ {
+ if (EX == -42)
+ {
+ z[1] = X[1] + TWO10;
+ z[2] = 0;
+ z[3] = 0;
+ k = 3;
+ }
+ else if (EX == -43)
+ {
+ z[1] = TWO10;
+ z[2] = X[1];
+ z[3] = 0;
+ k = 2;
+ }
+ else
+ {
+ z[1] = TWO10;
+ z[2] = 0;
+ z[3] = X[1];
+ k = 1;
+ }
+ }
+ else if (p2 == 2)
+ {
+ if (EX == -42)
+ {
+ z[1] = X[1] + TWO10;
+ z[2] = X[2];
+ z[3] = 0;
+ 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] = 0;
+ 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] = 0;
+ k = 1;
+ }
+ z[3] = X[k];
+ }
- if (X[0] == ZERO) {*y = ZERO; return; }
+ u = ALIGN_DOWN_TO (z[3], TWO5);
- if (EX> -42) norm(x,y,p);
- else if (EX==-42 && X[1]>=TWO10) norm(x,y,p);
- else denorm(x,y,p);
+ if (u == z[3])
+ {
+ for (i = k + 1; i <= p2; i++)
+ {
+ if (X[i] == 0)
+ continue;
+ else
+ {
+ z[3] += 1;
+ break;
+ }
+ }
+ }
+
+ c = X[0] * ((z[1] + R * (z[2] + R * z[3])) - TWO10);
+
+ *y = c * TWOM1032;
+# undef R
}
-#endif
+/* Convert multiple precision number *X into double precision number *Y. The
+ result is correctly rounded to the nearest/even. */
+void
+__mp_dbl (const mp_no *x, double *y, int p)
+{
+ if (X[0] == 0)
+ {
+ *y = 0;
+ return;
+ }
-/* 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. */
+ if (__glibc_likely (EX > -42 || (EX == -42 && X[1] >= TWO10)))
+ norm (x, y, p);
+ else
+ denorm (x, y, p);
+}
+#endif
+/* Get the multiple precision equivalent of X into *Y. If the precision is too
+ small, the result is truncated. */
void
SECTION
-__dbl_mp(double x, mp_no *y, int p) {
+__dbl_mp (double x, mp_no *y, int p)
+{
+ long i, n;
+ long p2 = p;
- int i,n;
- double u;
+ /* Sign. */
+ if (x == 0)
+ {
+ Y[0] = 0;
+ return;
+ }
+ else if (x > 0)
+ Y[0] = 1;
+ else
+ {
+ Y[0] = -1;
+ x = -x;
+ }
- /* 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;
+ /* Exponent. */
+ for (EY = 1; x >= RADIX; EY += 1)
+ x *= RADIXI;
+ for (; x < 1; EY -= 1)
+ x *= RADIX;
+
+ /* Digits. */
+ n = MIN (p2, 4);
+ for (i = 1; i <= n; i++)
+ {
+ INTEGER_OF (x, Y[i]);
+ x *= RADIX;
+ }
+ for (; i <= p2; i++)
+ Y[i] = 0;
}
-
-/* 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. */
-
+/* Add magnitudes of *X and *Y assuming that abs (*X) >= abs (*Y) > 0. The
+ sign of the sum *Z is not changed. X and Y may overlap but not X and Z or
+ Y and Z. No guard digit is used. The result equals the exact sum,
+ truncated. */
static void
SECTION
-add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
- int i,j,k;
+add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
+{
+ long i, j, k;
+ long p2 = p;
+ mantissa_t zk;
EZ = EX;
- i=p; j=p+ EY - EX; k=p+1;
+ i = p2;
+ j = p2 + EY - EX;
+ k = p2 + 1;
+
+ if (__glibc_unlikely (j < 1))
+ {
+ __cpy (x, z, p);
+ return;
+ }
- 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;
-}
+ zk = 0;
+ for (; j > 0; i--, j--)
+ {
+ zk += X[i] + Y[j];
+ if (zk >= RADIX)
+ {
+ Z[k--] = zk - RADIX;
+ zk = 1;
+ }
+ else
+ {
+ Z[k--] = zk;
+ zk = 0;
+ }
+ }
-/* 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. */
+ for (; i > 0; i--)
+ {
+ zk += X[i];
+ if (zk >= RADIX)
+ {
+ Z[k--] = zk - RADIX;
+ zk = 1;
+ }
+ else
+ {
+ Z[k--] = zk;
+ zk = 0;
+ }
+ }
+ if (zk == 0)
+ {
+ for (i = 1; i <= p2; i++)
+ Z[i] = Z[i + 1];
+ }
+ else
+ {
+ Z[1] = zk;
+ EZ += 1;
+ }
+}
+
+/* Subtract the magnitudes of *X and *Y assuming that abs (*x) > abs (*y) > 0.
+ The sign of the difference *Z is not changed. X and Y may overlap but not X
+ and Z or Y and Z. One guard digit is used. The error is less than one
+ ULP. */
static void
SECTION
-sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
- int i,j,k;
+sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
+{
+ long i, j, k;
+ long p2 = p;
+ mantissa_t zk;
EZ = EX;
+ i = p2;
+ j = p2 + EY - EX;
+ k = p2;
+
+ /* Y is too small compared to X, copy X over to the result. */
+ if (__glibc_unlikely (j < 1))
+ {
+ __cpy (x, z, p);
+ return;
+ }
+
+ /* The relevant least significant digit in Y is non-zero, so we factor it in
+ to enhance accuracy. */
+ if (j < p2 && Y[j + 1] > 0)
+ {
+ Z[k + 1] = RADIX - Y[j + 1];
+ zk = -1;
+ }
+ else
+ zk = Z[k + 1] = 0;
- 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;
- }
+ /* Subtract and borrow. */
+ for (; j > 0; i--, j--)
+ {
+ zk += (X[i] - Y[j]);
+ if (zk < 0)
+ {
+ Z[k--] = zk + RADIX;
+ zk = -1;
+ }
+ else
+ {
+ Z[k--] = zk;
+ zk = 0;
+ }
+ }
- for (i=1; Z[i] == ZERO; i++) ;
+ /* We're done with digits from Y, so it's just digits in X. */
+ for (; i > 0; i--)
+ {
+ zk += X[i];
+ if (zk < 0)
+ {
+ Z[k--] = zk + RADIX;
+ zk = -1;
+ }
+ else
+ {
+ Z[k--] = zk;
+ zk = 0;
+ }
+ }
+
+ /* Normalize. */
+ for (i = 1; Z[i] == 0; i++)
+ ;
EZ = EZ - i + 1;
- for (k=1; i <= p+1; )
+ for (k = 1; i <= p2 + 1; )
Z[k++] = Z[i++];
- for (; k <= p; )
- Z[k++] = ZERO;
+ for (; k <= p2; )
+ Z[k++] = 0;
}
-
-/* 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. */
-
+/* Add *X and *Y and store the result in *Z. X and Y may overlap, but not X
+ and Z or Y and Z. One guard digit is used. The error is less than one
+ ULP. */
void
SECTION
-__add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
+__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] == 0)
+ {
+ __cpy (y, z, p);
+ return;
+ }
+ else if (Y[0] == 0)
+ {
+ __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;
- }
+ 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] = 0;
+ }
}
+/* Subtract *Y from *X and return the result in *Z. X and Y may overlap but
+ not X and Z or Y and Z. One guard digit is used. The error is less than
+ one ULP. */
+void
+SECTION
+__sub (const mp_no *x, const mp_no *y, mp_no *z, int p)
+{
+ int n;
+
+ if (X[0] == 0)
+ {
+ __cpy (y, z, p);
+ Z[0] = -Z[0];
+ return;
+ }
+ else if (Y[0] == 0)
+ {
+ __cpy (x, z, p);
+ 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. */
+ 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] = 0;
+ }
+}
+#ifndef NO__MUL
+/* Multiply *X and *Y and store result in *Z. X and Y may overlap but not X
+ and Z or Y and Z. For P in [1, 2, 3], the exact result is truncated to P
+ digits. In case P > 3 the error is bounded by 1.001 ULP. */
void
SECTION
-__sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+__mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
+{
+ long i, j, k, ip, ip2;
+ long p2 = p;
+ mantissa_store_t zk;
+ const mp_no *a;
+ mantissa_store_t *diag;
+
+ /* Is z=0? */
+ if (__glibc_unlikely (X[0] * Y[0] == 0))
+ {
+ Z[0] = 0;
+ return;
+ }
- int n;
+ /* We need not iterate through all X's and Y's since it's pointless to
+ multiply zeroes. Here, both are zero... */
+ for (ip2 = p2; ip2 > 0; ip2--)
+ if (X[ip2] != 0 || Y[ip2] != 0)
+ break;
+
+ a = X[ip2] != 0 ? y : x;
+
+ /* ... and here, at least one of them is still zero. */
+ for (ip = ip2; ip > 0; ip--)
+ if (a->d[ip] != 0)
+ break;
+
+ /* The product looks like this for p = 3 (as an example):
+
+
+ a1 a2 a3
+ x b1 b2 b3
+ -----------------------------
+ a1*b3 a2*b3 a3*b3
+ a1*b2 a2*b2 a3*b2
+ a1*b1 a2*b1 a3*b1
+
+ So our K needs to ideally be P*2, but we're limiting ourselves to P + 3
+ for P >= 3. We compute the above digits in two parts; the last P-1
+ digits and then the first P digits. The last P-1 digits are a sum of
+ products of the input digits from P to P-k where K is 0 for the least
+ significant digit and increases as we go towards the left. The product
+ term is of the form X[k]*X[P-k] as can be seen in the above example.
+
+ The first P digits are also a sum of products with the same product term,
+ except that the sum is from 1 to k. This is also evident from the above
+ example.
+
+ Another thing that becomes evident is that only the most significant
+ ip+ip2 digits of the result are non-zero, where ip and ip2 are the
+ 'internal precision' of the input numbers, i.e. digits after ip and ip2
+ are all 0. */
+
+ k = (__glibc_unlikely (p2 < 3)) ? p2 + p2 : p2 + 3;
+
+ while (k > ip + ip2 + 1)
+ Z[k--] = 0;
+
+ zk = 0;
+
+ /* Precompute sums of diagonal elements so that we can directly use them
+ later. See the next comment to know we why need them. */
+ diag = alloca (k * sizeof (mantissa_store_t));
+ mantissa_store_t d = 0;
+ for (i = 1; i <= ip; i++)
+ {
+ d += X[i] * (mantissa_store_t) Y[i];
+ diag[i] = d;
+ }
+ while (i < k)
+ diag[i++] = d;
- if (X[0] == ZERO) {__cpy(y,z,p); Z[0] = -Z[0]; return; }
- else if (Y[0] == ZERO) {__cpy(x,z,p); return; }
+ while (k > p2)
+ {
+ long lim = k / 2;
- 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;
- }
-}
+ if (k % 2 == 0)
+ /* We want to add this only once, but since we subtract it in the sum
+ of products above, we add twice. */
+ zk += 2 * X[lim] * (mantissa_store_t) Y[lim];
+ for (i = k - p2, j = p2; i < j; i++, j--)
+ zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]);
-/* 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. */
+ zk -= diag[k - 1];
-void
-SECTION
-__mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+ DIV_RADIX (zk, Z[k]);
+ k--;
+ }
- int i, i1, i2, j, k, k2;
- double u;
+ /* The real deal. Mantissa digit Z[k] is the sum of all X[i] * Y[j] where i
+ goes from 1 -> k - 1 and j goes the same range in reverse. To reduce the
+ number of multiplications, we halve the range and if k is an even number,
+ add the diagonal element X[k/2]Y[k/2]. Through the half range, we compute
+ X[i] * Y[j] as (X[i] + X[j]) * (Y[i] + Y[j]) - X[i] * Y[i] - X[j] * Y[j].
+
+ This reduction tells us that we're summing two things, the first term
+ through the half range and the negative of the sum of the product of all
+ terms of X and Y in the full range. i.e.
+
+ SUM(X[i] * Y[i]) for k terms. This is precalculated above for each k in
+ a single loop so that it completes in O(n) time and can hence be directly
+ used in the loop below. */
+ while (k > 1)
+ {
+ long lim = k / 2;
+
+ if (k % 2 == 0)
+ /* We want to add this only once, but since we subtract it in the sum
+ of products above, we add twice. */
+ zk += 2 * X[lim] * (mantissa_store_t) Y[lim];
- /* 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;
+ for (i = 1, j = k - 1; i < j; i++, j--)
+ zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]);
+
+ zk -= diag[k - 1];
+
+ DIV_RADIX (zk, Z[k]);
+ k--;
+ }
+ Z[k] = zk;
+
+ /* Get the exponent sum into an intermediate variable. This is a subtle
+ optimization, where given enough registers, all operations on the exponent
+ happen in registers and the result is written out only once into EZ. */
+ int e = EX + EY;
+
+ /* Is there a carry beyond the most significant digit? */
+ if (__glibc_unlikely (Z[1] == 0))
+ {
+ for (i = 1; i <= p2; i++)
+ Z[i] = Z[i + 1];
+ e--;
+ }
+ EZ = e;
Z[0] = X[0] * Y[0];
}
+#endif
+
+#ifndef NO__SQR
+/* Square *X and store result in *Y. X and Y may not overlap. For P in
+ [1, 2, 3], the exact result is truncated to P digits. In case P > 3 the
+ error is bounded by 1.001 ULP. This is a faster special case of
+ multiplication. */
+void
+SECTION
+__sqr (const mp_no *x, mp_no *y, int p)
+{
+ long i, j, k, ip;
+ mantissa_store_t yk;
+ /* Is z=0? */
+ if (__glibc_unlikely (X[0] == 0))
+ {
+ Y[0] = 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. */
+ /* We need not iterate through all X's since it's pointless to
+ multiply zeroes. */
+ for (ip = p; ip > 0; ip--)
+ if (X[ip] != 0)
+ break;
-static
-SECTION
-void __inv(const mp_no *x, mp_no *y, int p) {
- int i;
-#if 0
- int l;
+ k = (__glibc_unlikely (p < 3)) ? p + p : p + 3;
+
+ while (k > 2 * ip + 1)
+ Y[k--] = 0;
+
+ yk = 0;
+
+ while (k > p)
+ {
+ mantissa_store_t yk2 = 0;
+ long lim = k / 2;
+
+ if (k % 2 == 0)
+ yk += X[lim] * (mantissa_store_t) X[lim];
+
+ /* In __mul, this loop (and the one within the next while loop) run
+ between a range to calculate the mantissa as follows:
+
+ Z[k] = X[k] * Y[n] + X[k+1] * Y[n-1] ... + X[n-1] * Y[k+1]
+ + X[n] * Y[k]
+
+ For X == Y, we can get away with summing halfway and doubling the
+ result. For cases where the range size is even, the mid-point needs
+ to be added separately (above). */
+ for (i = k - p, j = p; i < j; i++, j--)
+ yk2 += X[i] * (mantissa_store_t) X[j];
+
+ yk += 2 * yk2;
+
+ DIV_RADIX (yk, Y[k]);
+ k--;
+ }
+
+ while (k > 1)
+ {
+ mantissa_store_t yk2 = 0;
+ long lim = k / 2;
+
+ if (k % 2 == 0)
+ yk += X[lim] * (mantissa_store_t) X[lim];
+
+ /* Likewise for this loop. */
+ for (i = 1, j = k - 1; i < j; i++, j--)
+ yk2 += X[i] * (mantissa_store_t) X[j];
+
+ yk += 2 * yk2;
+
+ DIV_RADIX (yk, Y[k]);
+ k--;
+ }
+ Y[k] = yk;
+
+ /* Squares are always positive. */
+ Y[0] = 1;
+
+ /* Get the exponent sum into an intermediate variable. This is a subtle
+ optimization, where given enough registers, all operations on the exponent
+ happen in registers and the result is written out only once into EZ. */
+ int e = EX * 2;
+
+ /* Is there a carry beyond the most significant digit? */
+ if (__glibc_unlikely (Y[1] == 0))
+ {
+ for (i = 1; i <= p; i++)
+ Y[i] = Y[i + 1];
+ e--;
+ }
+
+ EY = e;
+}
#endif
+
+/* Invert *X and store in *Y. Relative error bound:
+ - For P = 2: 1.001 * R ^ (1 - P)
+ - For P = 3: 1.063 * R ^ (1 - P)
+ - For P > 3: 2.001 * R ^ (1 - P)
+
+ *X = 0 is not permissible. */
+static void
+SECTION
+__inv (const mp_no *x, mp_no *y, int p)
+{
+ long i;
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);
- }
+ 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
+ };
+
+ __cpy (x, &z, p);
+ z.e = 0;
+ __mp_dbl (&z, &t, p);
+ t = 1 / 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);
+ }
}
+/* Divide *X by *Y and store result in *Z. X and Y may overlap but not X and Z
+ or Y and Z. Relative error bound:
+ - For P = 2: 2.001 * R ^ (1 - P)
+ - For P = 3: 2.063 * R ^ (1 - P)
+ - For P > 3: 3.001 * R ^ (1 - P)
-/* 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. */
-
+ *X = 0 is not permissible. */
void
SECTION
-__dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
+__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);}
+ if (X[0] == 0)
+ Z[0] = 0;
+ else
+ {
+ __inv (y, &w, p);
+ __mul (x, &w, z, p);
+ }
}
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpa.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa.h
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* Written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation, Inc.
+ * Copyright (C) 2001-2017 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
@@ -23,36 +23,58 @@
/* FUNCTIONS: */
/* mcr */
/* acr */
-/* cr */
/* cpy */
-/* cpymn */
/* mp_dbl */
/* dbl_mp */
/* add */
/* sub */
/* mul */
-/* inv */
/* dvd */
/* */
/* Arithmetic functions for multiple precision numbers. */
/* Common types and definition */
/************************************************************************/
+#include <mpa-arch.h>
-typedef struct {/* This structure holds the details of a multi-precision */
- int e; /* floating point number, x: d[0] holds its sign (-1,0 or 1) */
- double d[40]; /* e holds its exponent (...,-2,-1,0,1,2,...) and */
-} mp_no; /* d[1]...d[p] hold its mantissa digits. The value of x is, */
- /* x = d[1]*r**(e-1) + d[2]*r**(e-2) + ... + d[p]*r**(e-p). */
- /* Here r = 2**24, 0 <= d[i] < r and 1 <= p <= 32. */
- /* p is a global variable. A multi-precision number is */
- /* always normalized. Namely, d[1] > 0. An exception is */
- /* a zero which is characterized by d[0] = 0. The terms */
- /* d[p+1], d[p+2], ... of a none zero number have no */
- /* significance and so are the terms e, d[1],d[2],... */
- /* of a zero. */
+/* The mp_no structure holds the details of a multi-precision floating point
+ number.
-typedef union { int i[2]; double d; } number;
+ - The radix of the number (R) is 2 ^ 24.
+
+ - E: The exponent of the number.
+
+ - D[0]: The sign (-1, 1) or 0 if the value is 0. In the latter case, the
+ values of the remaining members of the structure are ignored.
+
+ - D[1] - D[p]: The mantissa of the number where:
+
+ 0 <= D[i] < R and
+ P is the precision of the number and 1 <= p <= 32
+
+ D[p+1] ... D[39] have no significance.
+
+ - The value of the number is:
+
+ D[1] * R ^ (E - 1) + D[2] * R ^ (E - 2) ... D[p] * R ^ (E - p)
+
+ */
+typedef struct
+{
+ int e;
+ mantissa_t d[40];
+} mp_no;
+
+typedef union
+{
+ int i[2];
+ double d;
+} number;
+
+/* TODO: With only a partial backport of the constant cleanup we don't
+ define __mpone or __mptwo here for other code to use. */
+/* extern const mp_no __mpone;
+extern const mp_no __mptwo; */
#define X x->d
#define Y y->d
@@ -63,21 +85,73 @@ typedef union { int i[2]; double d; } nu
#define ABS(x) ((x) < 0 ? -(x) : (x))
-int __acr(const mp_no *, const mp_no *, int);
-// int __cr(const mp_no *, const mp_no *, int);
-void __cpy(const mp_no *, mp_no *, int);
-// void __cpymn(const mp_no *, int, mp_no *, int);
-void __mp_dbl(const mp_no *, double *, int);
-void __dbl_mp(double, mp_no *, int);
-void __add(const mp_no *, const mp_no *, mp_no *, int);
-void __sub(const mp_no *, const mp_no *, mp_no *, int);
-void __mul(const mp_no *, const mp_no *, mp_no *, int);
-// void __inv(const mp_no *, mp_no *, int);
-void __dvd(const mp_no *, const mp_no *, mp_no *, int);
+#ifndef RADIXI
+# define RADIXI 0x1.0p-24 /* 2^-24 */
+#endif
+
+#ifndef TWO52
+# define TWO52 0x1.0p52 /* 2^52 */
+#endif
+
+#define TWO5 TWOPOW (5) /* 2^5 */
+#define TWO8 TWOPOW (8) /* 2^52 */
+#define TWO10 TWOPOW (10) /* 2^10 */
+#define TWO18 TWOPOW (18) /* 2^18 */
+#define TWO19 TWOPOW (19) /* 2^19 */
+#define TWO23 TWOPOW (23) /* 2^23 */
+
+#define TWO57 0x1.0p57 /* 2^57 */
+#define TWO71 0x1.0p71 /* 2^71 */
+#define TWOM1032 0x1.0p-1032 /* 2^-1032 */
+#define TWOM1022 0x1.0p-1022 /* 2^-1022 */
+
+#define HALF 0x1.0p-1 /* 1/2 */
+#define MHALF -0x1.0p-1 /* -1/2 */
+#define HALFRAD 0x1.0p23 /* 2^23 */
+
+int __acr (const mp_no *, const mp_no *, int);
+void __cpy (const mp_no *, mp_no *, int);
+void __mp_dbl (const mp_no *, double *, int);
+void __dbl_mp (double, mp_no *, int);
+void __add (const mp_no *, const mp_no *, mp_no *, int);
+void __sub (const mp_no *, const mp_no *, mp_no *, int);
+void __mul (const mp_no *, const mp_no *, mp_no *, int);
+void __sqr (const mp_no *, mp_no *, int);
+void __dvd (const mp_no *, const mp_no *, mp_no *, int);
extern void __mpatan (mp_no *, mp_no *, int);
extern void __mpatan2 (mp_no *, mp_no *, mp_no *, int);
extern void __mpsqrt (mp_no *, mp_no *, int);
-extern void __mpexp (mp_no *, mp_no *__y, int);
+extern void __mpexp (mp_no *, mp_no *, int);
extern void __c32 (mp_no *, mp_no *, mp_no *, int);
extern int __mpranred (double, mp_no *, int);
+
+/* Given a power POW, build a multiprecision number 2^POW. */
+static inline void
+__pow_mp (int pow, mp_no *y, int p)
+{
+ int i, rem;
+
+ /* The exponent is E such that E is a factor of 2^24. The remainder (of the
+ form 2^x) goes entirely into the first digit of the mantissa as it is
+ always less than 2^24. */
+ EY = pow / 24;
+ rem = pow - EY * 24;
+ EY++;
+
+ /* If the remainder is negative, it means that POW was negative since
+ |EY * 24| <= |pow|. Adjust so that REM is positive and still less than
+ 24 because of which, the mantissa digit is less than 2^24. */
+ if (rem < 0)
+ {
+ EY--;
+ rem += 24;
+ }
+ /* The sign of any 2^x is always positive. */
+ Y[0] = 1;
+ Y[1] = 1 << rem;
+
+ /* Everything else is 0. */
+ for (i = 2; i <= p; i++)
+ Y[i] = 0;
+}
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpexp.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpexp.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpexp.c
@@ -69,13 +69,13 @@ __mpexp(mp_no *x, mp_no *y, int p) {
for (i=0; i<EX; i++) a *= RADIXI;
for ( ; i>EX; i--) a *= RADIX;
b = X[1]*RADIXI; m2 = 24*EX;
- for (; b<HALF; m2--) { a *= TWO; b *= TWO; }
+ for (; b<HALF; m2--) { a *= 2; b *= 2; }
if (b == HALF) {
- for (i=2; i<=p; i++) { if (X[i]!=ZERO) break; }
- if (i==p+1) { m2--; a *= TWO; }
+ for (i=2; i<=p; i++) { if (X[i]!=0) break; }
+ if (i==p+1) { m2--; a *= 2; }
}
if ((m=m1+m2) <= 0) {
- m=0; a=ONE;
+ m=0; a=1;
for (i=n-1; i>0; i--,n--) { if (m1np[i][p]+m2>0) break; }
}
@@ -84,8 +84,8 @@ __mpexp(mp_no *x, mp_no *y, int p) {
__mul(x,&mpt1,&mps,p);
/* Evaluate the polynomial. Put result in mpt2 */
- mpone.e=1; mpone.d[0]=ONE; mpone.d[1]=ONE;
- mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=__mpexp_nn[n].d;
+ mpone.e=1; mpone.d[0]=1; mpone.d[1]=1;
+ mpk.e = 1; mpk.d[0] = 1; mpk.d[1]=__mpexp_nn[n].d;
__dvd(&mps,&mpk,&mpt1,p);
__add(&mpone,&mpt1,&mpak,p);
for (k=n-1; k>1; k--) {
@@ -99,9 +99,9 @@ __mpexp(mp_no *x, mp_no *y, int p) {
/* Raise polynomial value to the power of 2**m. Put result in y */
for (k=0,j=0; k<m; ) {
- __mul(&mpt2,&mpt2,&mpt1,p); k++;
+ __sqr (&mpt2, &mpt1, p); k++;
if (k==m) { j=1; break; }
- __mul(&mpt1,&mpt1,&mpt2,p); k++;
+ __sqr (&mpt1, &mpt2, p); k++;
}
if (j) __cpy(&mpt1,y,p);
else __cpy(&mpt2,y,p);
Index: glibc-2.17-c758a686/sysdeps/x86_64/fpu/multiarch/mpa-avx.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/x86_64/fpu/multiarch/mpa-avx.c
+++ glibc-2.17-c758a686/sysdeps/x86_64/fpu/multiarch/mpa-avx.c
@@ -1,5 +1,6 @@
#define __add __add_avx
#define __mul __mul_avx
+#define __sqr __sqr_avx
#define __sub __sub_avx
#define __dbl_mp __dbl_mp_avx
#define __dvd __dvd_avx
Index: glibc-2.17-c758a686/sysdeps/x86_64/fpu/multiarch/mpa-fma4.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/x86_64/fpu/multiarch/mpa-fma4.c
+++ glibc-2.17-c758a686/sysdeps/x86_64/fpu/multiarch/mpa-fma4.c
@@ -1,5 +1,6 @@
#define __add __add_fma4
#define __mul __mul_fma4
+#define __sqr __sqr_fma4
#define __sub __sub_fma4
#define __dbl_mp __dbl_mp_fma4
#define __dvd __dvd_fma4
Index: glibc-2.17-c758a686/sysdeps/powerpc/power4/fpu/mpa.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/powerpc/power4/fpu/mpa.c
+++ glibc-2.17-c758a686/sysdeps/powerpc/power4/fpu/mpa.c
@@ -2,7 +2,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001, 2006 Free Software Foundation
+ * Copyright (C) 2001-2017 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
@@ -17,532 +17,198 @@
* 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: 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"
-#include <sys/param.h> /* For MIN() */
-/* 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) {
- long i;
- long p2 = p;
- for (i=1; i<=p2; 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) {
- long 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) {
- long 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) {
-
- long i,k;
- long n2 = n;
- long m2 = m;
-
- EY = EX; k=MIN(m2,n2);
- for (i=0; i <= k; i++) Y[i] = X[i];
- for ( ; i <= n2; 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
- long i;
-#if 0
- int k;
-#endif
- 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)
+/* Define __mul and __sqr and use the rest from generic code. */
+#define NO__MUL
+#define NO__SQR
+
+#include <sysdeps/ieee754/dbl-64/mpa.c>
+
+/* Multiply *X and *Y and store result in *Z. X and Y may overlap but not X
+ and Z or Y and Z. For P in [1, 2, 3], the exact result is truncated to P
+ digits. In case P > 3 the error is bounded by 1.001 ULP. */
+void
+__mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
{
- long i,k;
+ long i, i1, i2, j, k, k2;
long p2 = p;
- double c,u,z[5];
-#if 0
- double a,v;
-#endif
+ double u, zk, zk2;
-#define R radixi.d
- if (EX<-44 || (EX==-44 && X[1]<TWO5))
- { *y=ZERO; return; }
-
- if (p2==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 (p2==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 <= p2; i++) {
- if (X[i] == ZERO) continue;
- else {z[3] += ONE; break; }
+ /* Is z=0? */
+ if (__glibc_unlikely (X[0] * Y[0] == 0))
+ {
+ Z[0] = 0;
+ return;
}
- }
-
- 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) {
-#if 0
- int i,k;
- double a,c,u,v,z[5];
+ /* Multiply, add and carry */
+ k2 = (p2 < 3) ? p2 + p2 : p2 + 3;
+ zk = Z[k2] = 0;
+ for (k = k2; k > 1;)
+ {
+ if (k > p2)
+ {
+ i1 = k - p2;
+ i2 = p2 + 1;
+ }
+ else
+ {
+ i1 = 1;
+ i2 = k;
+ }
+#if 1
+ /* Rearrange this inner loop to allow the fmadd instructions to be
+ independent and execute in parallel on processors that have
+ dual symmetrical FP pipelines. */
+ if (i1 < (i2 - 1))
+ {
+ /* Make sure we have at least 2 iterations. */
+ if (((i2 - i1) & 1L) == 1L)
+ {
+ /* Handle the odd iterations case. */
+ zk2 = x->d[i2 - 1] * y->d[i1];
+ }
+ else
+ zk2 = 0.0;
+ /* Do two multiply/adds per loop iteration, using independent
+ accumulators; zk and zk2. */
+ for (i = i1, j = i2 - 1; i < i2 - 1; i += 2, j -= 2)
+ {
+ zk += x->d[i] * y->d[j];
+ zk2 += x->d[i + 1] * y->d[j - 1];
+ }
+ zk += zk2; /* Final sum. */
+ }
+ else
+ {
+ /* Special case when iterations is 1. */
+ zk += x->d[i1] * y->d[i1];
+ }
+#else
+ /* The original code. */
+ for (i = i1, j = i2 - 1; i < i2; i++, j--)
+ zk += X[i] * Y[j];
#endif
- 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) {
+ u = (zk + CUTTER) - CUTTER;
+ if (u > zk)
+ u -= RADIX;
+ Z[k] = zk - u;
+ zk = u * RADIXI;
+ --k;
+ }
+ Z[k] = zk;
- long i,n;
- long p2 = p;
- double u;
+ int e = EX + EY;
+ /* Is there a carry beyond the most significant digit? */
+ if (Z[1] == 0)
+ {
+ for (i = 1; i <= p2; i++)
+ Z[i] = Z[i + 1];
+ e--;
+ }
- /* 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(p2,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<=p2; i++) Y[i] = ZERO;
- return;
+ EZ = e;
+ Z[0] = X[0] * Y[0];
}
+/* Square *X and store result in *Y. X and Y may not overlap. For P in
+ [1, 2, 3], the exact result is truncated to P digits. In case P > 3 the
+ error is bounded by 1.001 ULP. This is a faster special case of
+ multiplication. */
+void
+__sqr (const mp_no *x, mp_no *y, int p)
+{
+ long i, j, k, ip;
+ double u, yk;
-/* 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. */
+ /* Is z=0? */
+ if (__glibc_unlikely (X[0] == 0))
+ {
+ Y[0] = 0;
+ return;
+ }
-static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+ /* We need not iterate through all X's since it's pointless to
+ multiply zeroes. */
+ for (ip = p; ip > 0; ip--)
+ if (X[ip] != 0)
+ break;
- long i,j,k;
- long p2 = p;
+ k = (__glibc_unlikely (p < 3)) ? p + p : p + 3;
- EZ = EX;
+ while (k > 2 * ip + 1)
+ Y[k--] = 0;
- i=p2; j=p2+ EY - EX; k=p2+1;
+ yk = 0;
- 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<=p2; i++) Z[i] = Z[i+1]; }
- else EZ += ONE;
-}
+ while (k > p)
+ {
+ double yk2 = 0.0;
+ long lim = k / 2;
+ if (k % 2 == 0)
+ {
+ yk += X[lim] * X[lim];
+ lim--;
+ }
-/* 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. */
+ /* In __mul, this loop (and the one within the next while loop) run
+ between a range to calculate the mantissa as follows:
-static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+ Z[k] = X[k] * Y[n] + X[k+1] * Y[n-1] ... + X[n-1] * Y[k+1]
+ + X[n] * Y[k]
- long i,j,k;
- long p2 = p;
+ For X == Y, we can get away with summing halfway and doubling the
+ result. For cases where the range size is even, the mid-point needs
+ to be added separately (above). */
+ for (i = k - p, j = p; i <= lim; i++, j--)
+ yk2 += X[i] * X[j];
- EZ = EX;
+ yk += 2.0 * yk2;
- if (EX == EY) {
- i=j=k=p2;
- Z[k] = Z[k+1] = ZERO; }
- else {
- j= EX - EY;
- if (j > p2) {__cpy(x,z,p); return; }
- else {
- i=p2; j=p2+1-j; k=p2;
- if (Y[j] > ZERO) {
- Z[k+1] = RADIX - Y[j--];
- Z[k] = MONE; }
- else {
- Z[k+1] = ZERO;
- Z[k] = ZERO; j--;}
+ u = (yk + CUTTER) - CUTTER;
+ if (u > yk)
+ u -= RADIX;
+ Y[k--] = yk - u;
+ yk = u * RADIXI;
}
- }
-
- 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 <= p2+1; )
- Z[k++] = Z[i++];
- for (; k <= p2; )
- 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) {
-
- long i, i1, i2, j, k, k2;
- long p2 = p;
- double u, zk, zk2;
-
- /* Is z=0? */
- if (X[0]*Y[0]==ZERO)
- { Z[0]=ZERO; return; }
-
- /* Multiply, add and carry */
- k2 = (p2<3) ? p2+p2 : p2+3;
- zk = Z[k2]=ZERO;
- for (k=k2; k>1; ) {
- if (k > p2) {i1=k-p2; i2=p2+1; }
- else {i1=1; i2=k; }
-#if 1
- /* rearange this inner loop to allow the fmadd instructions to be
- independent and execute in parallel on processors that have
- dual symetrical FP pipelines. */
- if (i1 < (i2-1))
+ while (k > 1)
{
- /* make sure we have at least 2 iterations */
- if (((i2 - i1) & 1L) == 1L)
- {
- /* Handle the odd iterations case. */
- zk2 = x->d[i2-1]*y->d[i1];
- }
- else
- zk2 = zero.d;
- /* Do two multiply/adds per loop iteration, using independent
- accumulators; zk and zk2. */
- for (i=i1,j=i2-1; i<i2-1; i+=2,j-=2)
- {
- zk += x->d[i]*y->d[j];
- zk2 += x->d[i+1]*y->d[j-1];
+ double yk2 = 0.0;
+ long lim = k / 2;
+
+ if (k % 2 == 0)
+ {
+ yk += X[lim] * X[lim];
+ lim--;
}
- zk += zk2; /* final sum. */
- }
- else
+
+ /* Likewise for this loop. */
+ for (i = 1, j = k - 1; i <= lim; i++, j--)
+ yk2 += X[i] * X[j];
+
+ yk += 2.0 * yk2;
+
+ u = (yk + CUTTER) - CUTTER;
+ if (u > yk)
+ u -= RADIX;
+ Y[k--] = yk - u;
+ yk = u * RADIXI;
+ }
+ Y[k] = yk;
+
+ /* Squares are always positive. */
+ Y[0] = 1.0;
+
+ int e = EX * 2;
+ /* Is there a carry beyond the most significant digit? */
+ if (__glibc_unlikely (Y[1] == 0))
{
- /* Special case when iterations is 1. */
- zk += x->d[i1]*y->d[i1];
+ for (i = 1; i <= p; i++)
+ Y[i] = Y[i + 1];
+ e--;
}
-#else
- /* The orginal code. */
- for (i=i1,j=i2-1; i<i2; i++,j--) zk += X[i]*Y[j];
-#endif
-
- u = (zk + CUTTER)-CUTTER;
- if (u > zk) u -= RADIX;
- Z[k] = zk - u;
- zk = u*RADIXI;
- --k;
- }
- Z[k] = zk;
-
- /* Is there a carry beyond the most significant digit? */
- if (Z[1] == ZERO) {
- for (i=1; i<=p2; 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) {
- long i;
-#if 0
- int l;
-#endif
- 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;
+ EY = e;
}
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/atnat.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/atnat.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/atnat.h
@@ -138,8 +138,6 @@
#endif
#endif
-#define ZERO zero.d
-#define ONE one.d
#define A a.d
#define B b.d
#define C c.d
@@ -160,7 +158,5 @@
#define U6 u6.d
#define U7 u7.d
#define U8 u8.d
-#define TWO8 two8.d
-#define TWO52 two52.d
#endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/atnat2.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/atnat2.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/atnat2.h
@@ -174,11 +174,4 @@
#endif
#endif
-#define ZERO zero.d
-#define MZERO mzero.d
-#define ONE one.d
-#define TWO8 two8.d
-#define TWO52 two52.d
-#define TWOM1022 twom1022.d
-
#endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa2.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpa2.h
+++ /dev/null
@@ -1,94 +0,0 @@
-
-/*
- * IBM Accurate Mathematical Library
- * Written by International Business Machines Corp.
- * Copyright (C) 2001 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:mpa2.h */
-/* */
-/* */
-/* variables prototype and definition according to type of processor */
-/* types definition */
-/**************************************************************************/
-
-#ifndef MPA2_H
-#define MPA2_H
-
-
-#ifdef BIG_ENDI
-static const number
-/**/ radix = {{0x41700000, 0x00000000} }, /* 2**24 */
-/**/ radixi = {{0x3e700000, 0x00000000} }, /* 2**-24 */
-/**/ cutter = {{0x44b00000, 0x00000000} }, /* 2**76 */
-/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */
-/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */
-/**/ mone = {{0xbff00000, 0x00000000} }, /* -1 */
-/**/ two = {{0x40000000, 0x00000000} }, /* 2 */
-/**/ two5 = {{0x40400000, 0x00000000} }, /* 2**5 */
-/**/ two10 = {{0x40900000, 0x00000000} }, /* 2**10 */
-/**/ two18 = {{0x41100000, 0x00000000} }, /* 2**18 */
-/**/ two19 = {{0x41200000, 0x00000000} }, /* 2**19 */
-/**/ two23 = {{0x41600000, 0x00000000} }, /* 2**23 */
-/**/ two52 = {{0x43300000, 0x00000000} }, /* 2**52 */
-/**/ two57 = {{0x43800000, 0x00000000} }, /* 2**57 */
-/**/ two71 = {{0x44600000, 0x00000000} }, /* 2**71 */
-/**/ twom1032 = {{0x00000400, 0x00000000} }; /* 2**-1032 */
-
-#else
-#ifdef LITTLE_ENDI
-static const number
-/**/ radix = {{0x00000000, 0x41700000} }, /* 2**24 */
-/**/ radixi = {{0x00000000, 0x3e700000} }, /* 2**-24 */
-/**/ cutter = {{0x00000000, 0x44b00000} }, /* 2**76 */
-/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */
-/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */
-/**/ mone = {{0x00000000, 0xbff00000} }, /* -1 */
-/**/ two = {{0x00000000, 0x40000000} }, /* 2 */
-/**/ two5 = {{0x00000000, 0x40400000} }, /* 2**5 */
-/**/ two10 = {{0x00000000, 0x40900000} }, /* 2**10 */
-/**/ two18 = {{0x00000000, 0x41100000} }, /* 2**18 */
-/**/ two19 = {{0x00000000, 0x41200000} }, /* 2**19 */
-/**/ two23 = {{0x00000000, 0x41600000} }, /* 2**23 */
-/**/ two52 = {{0x00000000, 0x43300000} }, /* 2**52 */
-/**/ two57 = {{0x00000000, 0x43800000} }, /* 2**57 */
-/**/ two71 = {{0x00000000, 0x44600000} }, /* 2**71 */
-/**/ twom1032 = {{0x00000000, 0x00000400} }; /* 2**-1032 */
-
-#endif
-#endif
-
-#define RADIX radix.d
-#define RADIXI radixi.d
-#define CUTTER cutter.d
-#define ZERO zero.d
-#define ONE one.d
-#define MONE mone.d
-#define TWO two.d
-#define TWO5 two5.d
-#define TWO10 two10.d
-#define TWO18 two18.d
-#define TWO19 two19.d
-#define TWO23 two23.d
-#define TWO52 two52.d
-#define TWO57 two57.d
-#define TWO71 two71.d
-#define TWOM1032 twom1032.d
-
-
-#endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpatan.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan.h
@@ -177,6 +177,3 @@ __atan_twonm1[33] = {
#endif
#endif
-
-#define ONE __atan_one.d
-#define TWO __atan_two.d
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan2.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpatan2.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan2.c
@@ -49,18 +49,16 @@ void
SECTION
__mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) {
- static const double ZERO = 0.0, ONE = 1.0;
-
mp_no mpone = {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,0.0,0.0}};
mp_no mpt1,mpt2,mpt3;
- if (X[0] <= ZERO) {
- mpone.e = 1; mpone.d[0] = mpone.d[1] = ONE;
+ if (X[0] <= 0) {
+ mpone.e = 1; mpone.d[0] = mpone.d[1] = 1;
__dvd(x,y,&mpt1,p); __mul(&mpt1,&mpt1,&mpt2,p);
- if (mpt1.d[0] != ZERO) mpt1.d[0] = ONE;
+ if (mpt1.d[0] != 0) mpt1.d[0] = 1;
__add(&mpt2,&mpone,&mpt3,p); __mpsqrt(&mpt3,&mpt2,p);
__add(&mpt1,&mpt2,&mpt3,p); mpt3.d[0]=Y[0];
__mpatan(&mpt3,&mpt1,p); __add(&mpt1,&mpt1,z,p);
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpexp.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpexp.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpexp.h
@@ -159,11 +159,4 @@ extern const number __mpexp_half attribu
#endif
#endif
-#define RADIX __mpexp_radix.d
-#define RADIXI __mpexp_radixi.d
-#define ZERO __mpexp_zero.d
-#define ONE __mpexp_one.d
-#define TWO __mpexp_two.d
-#define HALF __mpexp_half.d
-
#endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpsqrt.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpsqrt.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpsqrt.h
@@ -51,7 +51,4 @@ extern const int __mpsqrt_mp[33] attribu
4,4,4,4,4,4,4,4,4};
#endif
-#define ONE __mpsqrt_one.d
-#define HALFRAD __mpsqrt_halfrad.d
-
#endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mptan.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mptan.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mptan.c
@@ -47,8 +47,6 @@ void
SECTION
__mptan(double x, mp_no *mpy, int p) {
- static const double MONE = -1.0;
-
int n;
mp_no mpw, mpc, mps;
@@ -56,7 +54,7 @@ __mptan(double x, mp_no *mpy, int p) {
__c32(&mpw, &mpc, &mps, p); /* computing sin(x) and cos(x) */
if (n) /* second or fourth quarter of unit circle */
{ __dvd(&mpc,&mps,mpy,p);
- mpy->d[0] *= MONE;
+ mpy->d[0] *= -1;
} /* tan is negative in this area */
else __dvd(&mps,&mpc,mpy,p);
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/ulog.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/ulog.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/ulog.h
@@ -181,10 +181,6 @@
#endif
#endif
-#define ZERO zero.d
-#define ONE one.d
-#define HALF half.d
-#define MHALF mhalf.d
#define SQRT_2 sqrt_2.d
#define DEL_U delu.d
#define DEL_V delv.d
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/utan.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/utan.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/utan.h
@@ -270,10 +270,4 @@
#endif
#endif
-
-#define ZERO zero.d
-#define ONE one.d
-#define MONE mone.d
-#define TWO8 two8.d
-
#endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa-arch.h
===================================================================
--- /dev/null
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa-arch.h
@@ -0,0 +1,47 @@
+/* Overridable constants and operations.
+ Copyright (C) 2013 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/>. */
+
+#include <stdint.h>
+
+typedef long mantissa_t;
+typedef int64_t mantissa_store_t;
+
+#define TWOPOW(i) (1L << i)
+
+#define RADIX_EXP 24
+#define RADIX TWOPOW (RADIX_EXP) /* 2^24 */
+
+/* Divide D by RADIX and put the remainder in R. D must be a non-negative
+ integral value. */
+#define DIV_RADIX(d, r) \
+ ({ \
+ r = d & (RADIX - 1); \
+ d >>= RADIX_EXP; \
+ })
+
+/* Put the integer component of a double X in R and retain the fraction in
+ X. This is used in extracting mantissa digits for MP_NO by using the
+ integer portion of the current value of the number as the current mantissa
+ digit and then scaling by RADIX to get the next mantissa digit in the same
+ manner. */
+#define INTEGER_OF(x, i) \
+ ({ \
+ i = (mantissa_t) x; \
+ x -= i; \
+ })
+
+/* Align IN down to F. The code assumes that F is a power of two. */
+#define ALIGN_DOWN_TO(in, f) ((in) & -(f))
Index: glibc-2.17-c758a686/sysdeps/powerpc/power4/fpu/mpa-arch.h
===================================================================
--- /dev/null
+++ glibc-2.17-c758a686/sysdeps/powerpc/power4/fpu/mpa-arch.h
@@ -0,0 +1,56 @@
+/* Overridable constants and operations.
+ Copyright (C) 2013-2017 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/>. */
+
+typedef double mantissa_t;
+typedef double mantissa_store_t;
+
+#define TWOPOW(i) (0x1.0p##i)
+
+#define RADIX TWOPOW (24) /* 2^24 */
+#define CUTTER TWOPOW (76) /* 2^76 */
+#define RADIXI 0x1.0p-24 /* 2^-24 */
+#define TWO52 TWOPOW (52) /* 2^52 */
+
+/* Divide D by RADIX and put the remainder in R. */
+#define DIV_RADIX(d,r) \
+ ({ \
+ double u = ((d) + CUTTER) - CUTTER; \
+ if (u > (d)) \
+ u -= RADIX; \
+ r = (d) - u; \
+ (d) = u * RADIXI; \
+ })
+
+/* Put the integer component of a double X in R and retain the fraction in
+ X. */
+#define INTEGER_OF(x, r) \
+ ({ \
+ double u = ((x) + TWO52) - TWO52; \
+ if (u > (x)) \
+ u -= 1; \
+ (r) = u; \
+ (x) -= u; \
+ })
+
+/* Align IN down to a multiple of F, where F is a power of two. */
+#define ALIGN_DOWN_TO(in, f) \
+ ({ \
+ double factor = f * TWO52; \
+ double u = (in + factor) - factor; \
+ if (u > in) \
+ u -= f; \
+ u; \
+ })
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/e_atan2.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/e_atan2.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/e_atan2.c
@@ -98,15 +98,15 @@ __ieee754_atan2(double y,double x) {
/* y=+-0 */
if (uy==0x00000000) {
if (dy==0x00000000) {
- if ((ux&0x80000000)==0x00000000) return ZERO;
+ if ((ux&0x80000000)==0x00000000) return 0;
else return opi.d; } }
else if (uy==0x80000000) {
if (dy==0x00000000) {
- if ((ux&0x80000000)==0x00000000) return MZERO;
+ if ((ux&0x80000000)==0x00000000) return -0.0;
else return mopi.d;} }
/* x=+-0 */
- if (x==ZERO) {
+ if (x==0) {
if ((uy&0x80000000)==0x00000000) return hpi.d;
else return mhpi.d; }
@@ -118,8 +118,8 @@ __ieee754_atan2(double y,double x) {
else if (uy==0xfff00000) {
if (dy==0x00000000) return mqpi.d; }
else {
- if ((uy&0x80000000)==0x00000000) return ZERO;
- else return MZERO; }
+ if ((uy&0x80000000)==0x00000000) return 0;
+ else return -0.0; }
}
}
else if (ux==0xfff00000) {
@@ -141,14 +141,14 @@ __ieee754_atan2(double y,double x) {
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;
+ ax = (x<0) ? -x : x; ay = (y<0) ? -y : y;
de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
- if (de>=ep) { return ((y>ZERO) ? hpi.d : mhpi.d); }
+ if (de>=ep) { return ((y>0) ? hpi.d : mhpi.d); }
else if (de<=em) {
- if (x>ZERO) {
+ if (x>0) {
if ((z=ay/ax)<TWOM1022) return normalized(ax,ay,y,z);
else return signArctan2(y,z); }
- else { return ((y>ZERO) ? opi.d : mopi.d); } }
+ 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; }
@@ -170,7 +170,7 @@ __ieee754_atan2(double y,double x) {
EMULV(ay,u,v,vv,t1,t2,t3,t4,t5)
du=((ax-v)-vv)/ay; }
- if (x>ZERO) {
+ if (x>0) {
/* (i) x>0, abs(y)< abs(x): atan(ay/ax) */
if (ay<ax) {
@@ -180,7 +180,7 @@ __ieee754_atan2(double y,double x) {
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)
+ 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)
@@ -212,7 +212,7 @@ __ieee754_atan2(double y,double x) {
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)
+ 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)
@@ -237,7 +237,7 @@ __ieee754_atan2(double y,double x) {
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)
+ 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)
@@ -265,7 +265,7 @@ __ieee754_atan2(double y,double x) {
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)
+ 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)
@@ -293,7 +293,7 @@ __ieee754_atan2(double y,double x) {
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)
+ 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)
@@ -321,7 +321,7 @@ __ieee754_atan2(double y,double x) {
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)
+ 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)
@@ -347,7 +347,7 @@ __ieee754_atan2(double y,double x) {
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)
+ 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)
@@ -375,7 +375,7 @@ __ieee754_atan2(double y,double x) {
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)
+ 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)
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/e_log.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/e_log.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/e_log.c
@@ -78,9 +78,9 @@ __ieee754_log(double x) {
n=0;
if (__builtin_expect(ux < 0x00100000, 0)) {
if (__builtin_expect(((ux & 0x7fffffff) | dx) == 0, 0))
- return MHALF/ZERO; /* return -INF */
+ return MHALF/0; /* return -INF */
if (__builtin_expect(ux < 0, 0))
- return (x-x)/ZERO; /* return NaN */
+ return (x-x)/0; /* return NaN */
n -= 54; x *= two54.d; /* scale x */
num.d = x;
}
@@ -89,7 +89,7 @@ __ieee754_log(double x) {
/* Regular values of x */
- w = x-ONE;
+ w = x-1;
if (__builtin_expect(ABS(w) > U03, 1)) { goto case_03; }
@@ -113,25 +113,25 @@ __ieee754_log(double x) {
w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d))))))));
EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5)
ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2)
- MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(d9.d,dd9.d,s2,ss2,s3,ss3,t1,t2)
- MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(d8.d,dd8.d,s2,ss2,s3,ss3,t1,t2)
- MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(d7.d,dd7.d,s2,ss2,s3,ss3,t1,t2)
- MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(d6.d,dd6.d,s2,ss2,s3,ss3,t1,t2)
- MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(d5.d,dd5.d,s2,ss2,s3,ss3,t1,t2)
- MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(d4.d,dd4.d,s2,ss2,s3,ss3,t1,t2)
- MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(d3.d,dd3.d,s2,ss2,s3,ss3,t1,t2)
- MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(d2.d,dd2.d,s2,ss2,s3,ss3,t1,t2)
- MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
- MUL2(w,ZERO,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
- ADD2(w,ZERO, s3,ss3, b, bb,t1,t2)
+ MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(w,0,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(w,0, s3,ss3, b, bb,t1,t2)
/* End stage II, case abs(x-1) < 0.03 */
if ((y=b+(bb+b*E4)) == b+(bb-b*E4)) return y;
@@ -155,7 +155,7 @@ __ieee754_log(double x) {
j = (num.i[HIGH_HALF] & 0x000fffff) >> 4;
/* Compute w=(u-ui*vj)/(ui*vj) */
- p0=(ONE+(i-75)*DEL_U)*(ONE+(j-180)*DEL_V);
+ p0=(1+(i-75)*DEL_U)*(1+(j-180)*DEL_V);
q=u-p0; r0=Iu[i].d*Iv[j].d; w=q*r0;
/* Evaluate polynomial I */
@@ -178,11 +178,11 @@ __ieee754_log(double x) {
/* Improve the accuracy of r0 */
EMULV(p0,r0,sa,sb,t1,t2,t3,t4,t5)
- t=r0*((ONE-sa)-sb);
+ t=r0*((1-sa)-sb);
EADD(r0,t,ra,rb)
/* Compute w */
- MUL2(q,ZERO,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(q,0,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8)
EADD(A,B0,a0,aa0)
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/s_atan.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/s_atan.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/s_atan.c
@@ -82,7 +82,7 @@ double atan(double x) {
return x+x;
/* Regular values of x, including denormals +-0 and +-INF */
- u = (x<ZERO) ? -x : x;
+ u = (x<0) ? -x : x;
if (u<C) {
if (u<B) {
if (u<A) { /* u < A */
@@ -93,7 +93,7 @@ double atan(double x) {
EMULV(x,x,v,vv,t1,t2,t3,t4,t5) /* v+vv=x^2 */
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)
+ 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)
@@ -101,8 +101,8 @@ double atan(double x) {
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(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
- ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2)
+ MUL2(x,0,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(x,0,s2,ss2,s1,ss1,t1,t2)
if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1)) return y;
return atanMp(x,pr);
@@ -124,14 +124,14 @@ double atan(double x) {
z=u-hij[i][0].d;
s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
z*(hij[i][14].d+z* hij[i][15].d))));
- ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
- MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(hij[i][9].d,hij[i][10].d,s1,0,s2,ss2,t1,t2)
+ MUL2(z,0,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(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(z,0,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(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(z,0,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(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(z,0,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 ((y=s2+(ss2-U6*s2)) == s2+(ss2+U6*s2)) return __signArctan(x,y);
@@ -140,9 +140,9 @@ double atan(double x) {
}
else {
if (u<D) { /* C <= u < D */
- w=ONE/u;
+ w=1/u;
EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
- ww=w*((ONE-t1)-t2);
+ ww=w*((1-t1)-t2);
i=(TWO52+TWO8*w)-TWO52; i-=16;
z=(w-cij[i][0].d)+ww;
yy=HPI1-z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
@@ -152,12 +152,12 @@ double atan(double x) {
else u3=U32; /* w >= 1/2 */
if ((y=t1+(yy-u3)) == t1+(yy+u3)) return __signArctan(x,y);
- DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+ DIV2(1,0,u,0,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
t1=w-hij[i][0].d;
EADD(t1,ww,z,zz)
s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
z*(hij[i][14].d+z* hij[i][15].d))));
- ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+ ADD2(hij[i][9].d,hij[i][10].d,s1,0,s2,ss2,t1,t2)
MUL2(z,zz,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(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -173,18 +173,18 @@ double atan(double x) {
}
else {
if (u<E) { /* D <= u < E */
- w=ONE/u; v=w*w;
+ w=1/u; v=w*w;
EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
- ww=w*((ONE-t1)-t2);
+ ww=w*((1-t1)-t2);
ESUB(HPI,w,t3,cor)
yy=((HPI1+cor)-ww)-yy;
if ((y=t3+(yy-U4)) == t3+(yy+U4)) return __signArctan(x,y);
- DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+ DIV2(1,0,u,0,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
MUL2(w,ww,w,ww,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)
+ 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)
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/s_tan.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/s_tan.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/s_tan.c
@@ -84,7 +84,7 @@ tan(double x) {
goto ret;
}
- w=(x<ZERO) ? -x : x;
+ w=(x<0) ? -x : x;
/* (I) The case abs(x) <= 1.259e-8 */
if (w<=g1.d) { retval = x; goto ret; }
@@ -125,11 +125,11 @@ tan(double x) {
/* First stage */
i = ((int) (mfftnhf.d+TWO8*w));
- z = w-xfg[i][0].d; z2 = z*z; s = (x<ZERO) ? MONE : ONE;
+ z = w-xfg[i][0].d; z2 = z*z; s = (x<0) ? -1 : 1;
pz = z+z*z2*(e0.d+z2*e1.d);
fi = xfg[i][1].d; gi = xfg[i][2].d; t2 = pz*(gi+fi)/(gi-pz);
if ((y=fi+(t2-fi*u3.d))==fi+(t2+fi*u3.d)) { retval = (s*y); goto ret; }
- t3 = (t2<ZERO) ? -t2 : t2;
+ t3 = (t2<0) ? -t2 : t2;
t4 = fi*ua3.d+t3*ub3.d;
if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (s*y); goto ret; }
@@ -165,8 +165,8 @@ tan(double x) {
da = xn*mp3.d;
a=t1-da;
da = (t1-a)-da;
- if (a<ZERO) {ya=-a; yya=-da; sy=MONE;}
- else {ya= a; yya= da; sy= ONE;}
+ if (a<0) {ya=-a; yya=-da; sy=-1;}
+ else {ya= a; yya= da; sy= 1;}
/* (IV),(V) The case 0.787 < abs(x) <= 25, abs(y) <= 1e-7 */
if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
@@ -240,14 +240,14 @@ tan(double x) {
/* -cot */
t2 = pz*(fi+gi)/(fi+pz);
if ((y=gi-(t2-gi*u10.d))==gi-(t2+gi*u10.d)) { retval = (-sy*y); goto ret; }
- t3 = (t2<ZERO) ? -t2 : t2;
+ t3 = (t2<0) ? -t2 : t2;
t4 = gi*ua10.d+t3*ub10.d;
if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
else {
/* tan */
t2 = pz*(gi+fi)/(gi-pz);
if ((y=fi+(t2-fi*u9.d))==fi+(t2+fi*u9.d)) { retval = (sy*y); goto ret; }
- t3 = (t2<ZERO) ? -t2 : t2;
+ t3 = (t2<0) ? -t2 : t2;
t4 = fi*ua9.d+t3*ub9.d;
if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
@@ -295,8 +295,8 @@ tan(double x) {
a = t - t1;
da = ((t-a)-t1)+da;
EADD(a,da,t1,t2) a=t1; da=t2;
- if (a<ZERO) {ya=-a; yya=-da; sy=MONE;}
- else {ya= a; yya= da; sy= ONE;}
+ if (a<0) {ya=-a; yya=-da; sy=-1;}
+ else {ya= a; yya= da; sy= 1;}
/* (+++) The case 25 < abs(x) <= 1e8, abs(y) <= 1e-7 */
if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
@@ -355,14 +355,14 @@ tan(double x) {
/* -cot */
t2 = pz*(fi+gi)/(fi+pz);
if ((y=gi-(t2-gi*u18.d))==gi-(t2+gi*u18.d)) { retval = (-sy*y); goto ret; }
- t3 = (t2<ZERO) ? -t2 : t2;
+ t3 = (t2<0) ? -t2 : t2;
t4 = gi*ua18.d+t3*ub18.d;
if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
else {
/* tan */
t2 = pz*(gi+fi)/(gi-pz);
if ((y=fi+(t2-fi*u17.d))==fi+(t2+fi*u17.d)) { retval = (sy*y); goto ret; }
- t3 = (t2<ZERO) ? -t2 : t2;
+ t3 = (t2<0) ? -t2 : t2;
t4 = fi*ua17.d+t3*ub17.d;
if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
@@ -398,8 +398,8 @@ tan(double x) {
/* Range reduction by algorithm iii */
n = (__branred(x,&a,&da)) & 0x00000001;
EADD(a,da,t1,t2) a=t1; da=t2;
- if (a<ZERO) {ya=-a; yya=-da; sy=MONE;}
- else {ya= a; yya= da; sy= ONE;}
+ if (a<0) {ya=-a; yya=-da; sy=-1;}
+ else {ya= a; yya= da; sy= 1;}
/* (+++) The case 1e8 < abs(x) < 2**1024, abs(y) <= 1e-7 */
if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
@@ -463,14 +463,14 @@ tan(double x) {
/* -cot */
t2 = pz*(fi+gi)/(fi+pz);
if ((y=gi-(t2-gi*u26.d))==gi-(t2+gi*u26.d)) { retval = (-sy*y); goto ret; }
- t3 = (t2<ZERO) ? -t2 : t2;
+ t3 = (t2<0) ? -t2 : t2;
t4 = gi*ua26.d+t3*ub26.d;
if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
else {
/* tan */
t2 = pz*(gi+fi)/(gi-pz);
if ((y=fi+(t2-fi*u25.d))==fi+(t2+fi*u25.d)) { retval = (sy*y); goto ret; }
- t3 = (t2<ZERO) ? -t2 : t2;
+ t3 = (t2<0) ? -t2 : t2;
t4 = fi*ua25.d+t3*ub25.d;
if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpatan.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan.c
@@ -69,8 +69,8 @@ __mpatan(mp_no *x, mp_no *y, int p) {
{if (dx>__atan_xm[m].d) break;}
}
mpone.e = mptwo.e = mptwoim1.e = 1;
- mpone.d[0] = mpone.d[1] = mptwo.d[0] = mptwoim1.d[0] = ONE;
- mptwo.d[1] = TWO;
+ mpone.d[0] = mpone.d[1] = mptwo.d[0] = mptwoim1.d[0] = 1;
+ mptwo.d[1] = 2;
/* Reduce x m times */
__mul(x,x,&mpsm,p);
@@ -92,7 +92,7 @@ __mpatan(mp_no *x, mp_no *y, int p) {
n=__atan_np[p]; mptwoim1.d[1] = __atan_twonm1[p].d;
__dvd(&mpsm,&mptwoim1,&mpt,p);
for (i=n-1; i>1; i--) {
- mptwoim1.d[1] -= TWO;
+ mptwoim1.d[1] -= 2;
__dvd(&mpsm,&mptwoim1,&mpt1,p);
__mul(&mpsm,&mpt,&mpt2,p);
__sub(&mpt1,&mpt2,&mpt,p);
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpsqrt.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpsqrt.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpsqrt.c
@@ -62,8 +62,8 @@ __mpsqrt(mp_no *x, mp_no *y, int p) {
mp_no mpxn,mpz,mpu,mpt1,mpt2;
/* Prepare multi-precision 1/2 and 3/2 */
- mphalf.e =0; mphalf.d[0] =ONE; mphalf.d[1] =HALFRAD;
- mp3halfs.e=1; mp3halfs.d[0]=ONE; mp3halfs.d[1]=ONE; mp3halfs.d[2]=HALFRAD;
+ mphalf.e =0; mphalf.d[0] =1; mphalf.d[1] =HALFRAD;
+ mp3halfs.e=1; mp3halfs.d[0]=1; mp3halfs.d[1]=1; mp3halfs.d[2]=HALFRAD;
ey=EX/2; __cpy(x,&mpxn,p); mpxn.e -= (ey+ey);
__mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p);