aboutsummaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/dbl-64
diff options
context:
space:
mode:
authorUlrich Drepper <drepper@gmail.com>2011-10-25 00:56:33 -0400
committerUlrich Drepper <drepper@gmail.com>2011-10-25 00:56:33 -0400
commit31d3cc00b0cc5205b4b4efd73c911cfddff444c6 (patch)
tree89e95646d591fee083d62642085da77a01127292 /sysdeps/ieee754/dbl-64
parent202c9deb15ee43bcbe70b36fa9bae050b8633c27 (diff)
downloadglibc-31d3cc00b0cc5205b4b4efd73c911cfddff444c6.tar
glibc-31d3cc00b0cc5205b4b4efd73c911cfddff444c6.tar.gz
glibc-31d3cc00b0cc5205b4b4efd73c911cfddff444c6.tar.bz2
glibc-31d3cc00b0cc5205b4b4efd73c911cfddff444c6.zip
Cleanup FMA4 patch
Move the FMA4 code into its own section. Avoid some of the duplication of data resulting from the double use of source files.
Diffstat (limited to 'sysdeps/ieee754/dbl-64')
-rw-r--r--sysdeps/ieee754/dbl-64/branred.c10
-rw-r--r--sysdeps/ieee754/dbl-64/doasin.c8
-rw-r--r--sysdeps/ieee754/dbl-64/dosincos.c16
-rw-r--r--sysdeps/ieee754/dbl-64/e_asin.c12
-rw-r--r--sysdeps/ieee754/dbl-64/e_atan2.c16
-rw-r--r--sysdeps/ieee754/dbl-64/e_exp.c12
-rw-r--r--sysdeps/ieee754/dbl-64/e_log.c8
-rw-r--r--sysdeps/ieee754/dbl-64/e_pow.c24
-rw-r--r--sysdeps/ieee754/dbl-64/halfulp.c8
-rw-r--r--sysdeps/ieee754/dbl-64/mpa.c61
-rw-r--r--sysdeps/ieee754/dbl-64/mpa.h2
-rw-r--r--sysdeps/ieee754/dbl-64/mpatan.c29
-rw-r--r--sysdeps/ieee754/dbl-64/mpatan.h73
-rw-r--r--sysdeps/ieee754/dbl-64/mpatan2.c15
-rw-r--r--sysdeps/ieee754/dbl-64/mpexp.c43
-rw-r--r--sysdeps/ieee754/dbl-64/mpexp.h70
-rw-r--r--sysdeps/ieee754/dbl-64/mpsqrt.c18
-rw-r--r--sysdeps/ieee754/dbl-64/mpsqrt.h29
-rw-r--r--sysdeps/ieee754/dbl-64/mptan.c11
-rw-r--r--sysdeps/ieee754/dbl-64/s_sin.c64
-rw-r--r--sysdeps/ieee754/dbl-64/s_tan.c12
-rw-r--r--sysdeps/ieee754/dbl-64/sincos32.c50
-rw-r--r--sysdeps/ieee754/dbl-64/slowexp.c12
-rw-r--r--sysdeps/ieee754/dbl-64/slowpow.c12
24 files changed, 418 insertions, 197 deletions
diff --git a/sysdeps/ieee754/dbl-64/branred.c b/sysdeps/ieee754/dbl-64/branred.c
index 76015f0c5c..c8483034af 100644
--- a/sysdeps/ieee754/dbl-64/branred.c
+++ b/sysdeps/ieee754/dbl-64/branred.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* Written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation, Inc.
+ * Copyright (C) 2001, 2011 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
@@ -38,6 +38,10 @@
#include "branred.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
/*******************************************************************/
/* Routine branred() performs range reduction of a double number */
@@ -45,7 +49,9 @@
/* x=n*pi/2+(a+aa), abs(a+aa)<pi/4, n=0,+-1,+-2,.... */
/* Routine return integer (n mod 4) */
/*******************************************************************/
-int __branred(double x, double *a, double *aa)
+int
+SECTION
+__branred(double x, double *a, double *aa)
{
int i,k;
#if 0
diff --git a/sysdeps/ieee754/dbl-64/doasin.c b/sysdeps/ieee754/dbl-64/doasin.c
index 64abc3cbb1..14958b5ca2 100644
--- a/sysdeps/ieee754/dbl-64/doasin.c
+++ b/sysdeps/ieee754/dbl-64/doasin.c
@@ -34,11 +34,17 @@
#include <dla.h>
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
/********************************************************************/
/* Compute arcsin(x,dx,v) of double-length number (x+dx) the result */
/* stored in v where v= v[0]+v[1] =arcsin(x+dx) */
/********************************************************************/
-void __doasin(double x, double dx, double v[]) {
+void
+SECTION
+__doasin(double x, double dx, double v[]) {
#include "doasin.h"
diff --git a/sysdeps/ieee754/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c
index 712d585b9e..e8890ff8de 100644
--- a/sysdeps/ieee754/dbl-64/dosincos.c
+++ b/sysdeps/ieee754/dbl-64/dosincos.c
@@ -39,6 +39,10 @@
#include "dosincos.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
extern const union
{
int4 i[880];
@@ -52,7 +56,9 @@ extern const union
/*(x+dx) between 0 and PI/4 */
/***********************************************************************/
-void __dubsin(double x, double dx, double v[]) {
+void
+SECTION
+__dubsin(double x, double dx, double v[]) {
double r,s,c,cc,d,dd,d2,dd2,e,ee,
sn,ssn,cs,ccs,ds,dss,dc,dcc;
#ifndef DLA_FMS
@@ -106,7 +112,9 @@ void __dubsin(double x, double dx, double v[]) {
/*(x+dx) between 0 and PI/4 */
/**********************************************************************/
-void __dubcos(double x, double dx, double v[]) {
+void
+SECTION
+__dubcos(double x, double dx, double v[]) {
double r,s,c,cc,d,dd,d2,dd2,e,ee,
sn,ssn,cs,ccs,ds,dss,dc,dcc;
#ifndef DLA_FMS
@@ -172,7 +180,9 @@ void __dubcos(double x, double dx, double v[]) {
/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
/* as Double-Length number and store it in array v */
/**********************************************************************/
-void __docos(double x, double dx, double v[]) {
+void
+SECTION
+__docos(double x, double dx, double v[]) {
double y,yy,p,w[2];
if (x>0) {y=x; yy=dx;}
else {y=-x; yy=-dx;}
diff --git a/sysdeps/ieee754/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c
index cd4cc2e2c2..65319c0b58 100644
--- a/sysdeps/ieee754/dbl-64/e_asin.c
+++ b/sysdeps/ieee754/dbl-64/e_asin.c
@@ -42,6 +42,10 @@
#include "uasncs.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
void __doasin(double x, double dx, double w[]);
void __dubsin(double x, double dx, double v[]);
void __dubcos(double x, double dx, double v[]);
@@ -53,7 +57,9 @@ double __cos32(double x, double res, double res1);
/* An ultimate asin routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of arcsin(x) */
/***************************************************************************/
-double __ieee754_asin(double x){
+double
+SECTION
+__ieee754_asin(double x){
double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2];
mynumber u,v;
int4 k,m,n;
@@ -334,7 +340,9 @@ strong_alias (__ieee754_asin, __asin_finite)
/* */
/*******************************************************************/
-double __ieee754_acos(double x)
+double
+SECTION
+__ieee754_acos(double x)
{
double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps;
#if 0
diff --git a/sysdeps/ieee754/dbl-64/e_atan2.c b/sysdeps/ieee754/dbl-64/e_atan2.c
index 9caacccf4c..64dae3e8d5 100644
--- a/sysdeps/ieee754/dbl-64/e_atan2.c
+++ b/sysdeps/ieee754/dbl-64/e_atan2.c
@@ -44,6 +44,10 @@
#include "atnat2.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
/************************************************************************/
/* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */
/* it computes the correctly rounded (to nearest) value of atan2(y,x). */
@@ -59,7 +63,9 @@ static double signArctan2(double y,double z)
static double normalized(double ,double,double ,double);
void __mpatan2(mp_no *,mp_no *,mp_no *,int);
-double __ieee754_atan2(double y,double x) {
+double
+SECTION
+__ieee754_atan2(double y,double x) {
int i,de,ux,dx,uy,dy;
#if 0
@@ -384,7 +390,9 @@ strong_alias (__ieee754_atan2, __atan2_finite)
#endif
/* Treat the Denormalized case */
-static double normalized(double ax,double ay,double y, double z)
+static double
+SECTION
+normalized(double ax,double ay,double y, double z)
{ int p;
mp_no mpx,mpy,mpz,mperr,mpz2,mpt1;
p=6;
@@ -394,7 +402,9 @@ static double normalized(double ax,double ay,double y, double z)
return signArctan2(y,z);
}
/* Stage 3: Perform a multi-Precision computation */
-static double atan2Mp(double x,double y,const int pr[])
+static double
+SECTION
+atan2Mp(double x,double y,const int pr[])
{
double z1,z2;
int i,p;
diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c
index 48bbb05ed8..e7a839d42e 100644
--- a/sysdeps/ieee754/dbl-64/e_exp.c
+++ b/sysdeps/ieee754/dbl-64/e_exp.c
@@ -40,13 +40,19 @@
#include "uexp.tbl"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
double __slowexp(double);
/***************************************************************************/
/* An ultimate exp routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of e^x */
/***************************************************************************/
-double __ieee754_exp(double x) {
+double
+SECTION
+__ieee754_exp(double x) {
double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
mynumber junk1, junk2, binexp = {{0,0}};
#if 0
@@ -156,7 +162,9 @@ strong_alias (__ieee754_exp, __exp_finite)
/*else return e^(x + xx) (always positive ) */
/************************************************************************/
-double __exp1(double x, double xx, double error) {
+double
+SECTION
+__exp1(double x, double xx, double error) {
double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
mynumber junk1, junk2, binexp = {{0,0}};
#if 0
diff --git a/sysdeps/ieee754/dbl-64/e_log.c b/sysdeps/ieee754/dbl-64/e_log.c
index 7a0a26f251..e45520eba8 100644
--- a/sysdeps/ieee754/dbl-64/e_log.c
+++ b/sysdeps/ieee754/dbl-64/e_log.c
@@ -41,13 +41,19 @@
#include "MathLib.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
void __mplog(mp_no *, mp_no *, int);
/*********************************************************************/
/* An ultimate log routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of log(x). */
/*********************************************************************/
-double __ieee754_log(double x) {
+double
+SECTION
+__ieee754_log(double x) {
#define M 4
static const int pr[M]={8,10,18,32};
int i,j,n,ux,dx,p;
diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index 94b1ab8961..350e93986d 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -43,6 +43,10 @@
#include "upow.tbl"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
double __exp1(double x, double xx, double error);
static double log1(double x, double *delta, double *error);
@@ -55,7 +59,9 @@ static int checkint(double x);
/* An ultimate power routine. Given two IEEE double machine numbers y,x */
/* it computes the correctly rounded (to nearest) value of X^y. */
/***************************************************************************/
-double __ieee754_pow(double x, double y) {
+double
+SECTION
+__ieee754_pow(double x, double y) {
double z,a,aa,error, t,a1,a2,y1,y2;
#if 0
double gor=1.0;
@@ -160,7 +166,9 @@ strong_alias (__ieee754_pow, __pow_finite)
/**************************************************************************/
/* Computing x^y using more accurate but more slow log routine */
/**************************************************************************/
-static double power1(double x, double y) {
+static double
+SECTION
+power1(double x, double y) {
double z,a,aa,error, t,a1,a2,y1,y2;
z = my_log2(x,&aa,&error);
t = y*134217729.0;
@@ -183,7 +191,9 @@ static double power1(double x, double y) {
/* + the parameter delta. */
/* The result is bounded by error (rightmost argument) */
/****************************************************************************/
-static double log1(double x, double *delta, double *error) {
+static double
+SECTION
+log1(double x, double *delta, double *error) {
int i,j,m;
#if 0
int n;
@@ -275,7 +285,9 @@ static double log1(double x, double *delta, double *error) {
/* Computing log(x)(x is left argument).The result is return double + delta.*/
/* The result is bounded by error (right argument) */
/****************************************************************************/
-static double my_log2(double x, double *delta, double *error) {
+static double
+SECTION
+my_log2(double x, double *delta, double *error) {
int i,j,m;
#if 0
int n;
@@ -369,7 +381,9 @@ static double my_log2(double x, double *delta, double *error) {
/* Routine receives a double x and checks if it is an integer. If not */
/* it returns 0, else it returns 1 if even or -1 if odd. */
/**********************************************************************/
-static int checkint(double x) {
+static int
+SECTION
+checkint(double x) {
union {int4 i[2]; double x;} u;
int k,m,n;
#if 0
diff --git a/sysdeps/ieee754/dbl-64/halfulp.c b/sysdeps/ieee754/dbl-64/halfulp.c
index 31bd2daf4d..6018309427 100644
--- a/sysdeps/ieee754/dbl-64/halfulp.c
+++ b/sysdeps/ieee754/dbl-64/halfulp.c
@@ -40,6 +40,10 @@
#include <dla.h>
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
static const int4 tab54[32] = {
262143, 11585, 1782, 511, 210, 107, 63, 42,
30, 22, 17, 14, 12, 10, 9, 7,
@@ -47,7 +51,9 @@ static const int4 tab54[32] = {
3, 3, 3, 3, 3, 3, 3, 3 };
-double __halfulp(double x, double y)
+double
+SECTION
+__halfulp(double x, double y)
{
mynumber v;
double z,u,uu;
diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c
index ad5a639c4b..39c640882b 100644
--- a/sysdeps/ieee754/dbl-64/mpa.c
+++ b/sysdeps/ieee754/dbl-64/mpa.c
@@ -47,11 +47,18 @@
#include "mpa.h"
#include "mpa2.h"
#include <sys/param.h> /* For MIN() */
+
+#ifndef SECTION
+# define SECTION
+#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. */
-static int mcr(const mp_no *x, const mp_no *y, int p) {
+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;
@@ -61,9 +68,9 @@ static int mcr(const mp_no *x, const mp_no *y, int p) {
}
-
/* acr() compares the absolute values of two multiple precision numbers */
-static int __acr(const mp_no *x, const mp_no *y, int p) {
+int
+__acr(const mp_no *x, const mp_no *y, int p) {
int i;
if (X[0] == ZERO) {
@@ -79,10 +86,11 @@ static int __acr(const mp_no *x, const mp_no *y, int p) {
return i;
}
+#endif
#if 0
-/* cr90 compares the values of two multiple precision numbers */
+/* cr() compares the values of two multiple precision numbers */
static int __cr(const mp_no *x, const mp_no *y, int p) {
int i;
@@ -119,8 +127,6 @@ static void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
EY = EX; k=MIN(m,n);
for (i=0; i <= k; i++) Y[i] = X[i];
for ( ; i <= n; i++) Y[i] = ZERO;
-
- return;
}
#endif
@@ -177,7 +183,6 @@ static void norm(const mp_no *x, double *y, int p)
for (i=1; i>EX; i--) c *= RADIXI;
*y = c;
- return;
#undef R
}
@@ -225,8 +230,6 @@ static void denorm(const mp_no *x, double *y, int p)
c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10);
*y = c*TWOM1032;
- return;
-
#undef R
}
@@ -252,7 +255,9 @@ void __mp_dbl(const mp_no *x, double *y, int p) {
/* 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) {
+void
+SECTION
+__dbl_mp(double x, mp_no *y, int p) {
int i,n;
double u;
@@ -273,7 +278,6 @@ void __dbl_mp(double x, mp_no *y, int p) {
if (u>x) u -= ONE;
Y[i] = u; x -= u; x *= RADIX; }
for ( ; i<=p; i++) Y[i] = ZERO;
- return;
}
@@ -283,7 +287,9 @@ void __dbl_mp(double x, mp_no *y, int p) {
/* 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) {
+static void
+SECTION
+add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int i,j,k;
@@ -325,7 +331,9 @@ static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* 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) {
+static void
+SECTION
+sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int i,j,k;
@@ -372,8 +380,6 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
Z[k++] = Z[i++];
for (; k <= p; )
Z[k++] = ZERO;
-
- return;
}
@@ -381,7 +387,9 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* 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) {
+void
+SECTION
+__add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int n;
@@ -397,7 +405,6 @@ void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; }
else Z[0] = ZERO;
}
- return;
}
@@ -405,7 +412,9 @@ void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* 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) {
+void
+SECTION
+__sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int n;
@@ -421,7 +430,6 @@ void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = -Y[0]; }
else Z[0] = ZERO;
}
- return;
}
@@ -430,7 +438,9 @@ void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* 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) {
+void
+SECTION
+__mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int i, i1, i2, j, k, k2;
double u;
@@ -461,7 +471,6 @@ void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
EZ = EX + EY;
Z[0] = X[0] * Y[0];
- return;
}
@@ -470,7 +479,9 @@ void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* 2.001*r**(1-p) for p>3. */
/* *x=0 is not permissible. *x is left unchanged. */
-static void __inv(const mp_no *x, mp_no *y, int p) {
+static
+SECTION
+void __inv(const mp_no *x, mp_no *y, int p) {
int i;
#if 0
int l;
@@ -493,7 +504,6 @@ static void __inv(const mp_no *x, mp_no *y, int p) {
__sub(&mptwo,y,&z,p);
__mul(&w,&z,y,p);
}
- return;
}
@@ -502,11 +512,12 @@ static void __inv(const mp_no *x, mp_no *y, int p) {
/* 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) {
+void
+SECTION
+__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;
}
diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h
index 3ca0ca5f00..5647ab7b4f 100644
--- a/sysdeps/ieee754/dbl-64/mpa.h
+++ b/sysdeps/ieee754/dbl-64/mpa.h
@@ -64,7 +64,7 @@ typedef union { int i[2]; double d; } number;
#define ABS(x) ((x) < 0 ? -(x) : (x))
-// int __acr(const mp_no *, const mp_no *, int);
+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);
diff --git a/sysdeps/ieee754/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c
index ee21c25138..f40873ea59 100644
--- a/sysdeps/ieee754/dbl-64/mpatan.c
+++ b/sysdeps/ieee754/dbl-64/mpatan.c
@@ -1,8 +1,7 @@
-
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
*
* 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
@@ -34,11 +33,19 @@
#include "endian.h"
#include "mpa.h"
-void __mpsqrt(mp_no *, mp_no *, int);
-void __mpatan(mp_no *x, mp_no *y, int p) {
+#ifndef SECTION
+# define SECTION
+#endif
+
#include "mpatan.h"
+void __mpsqrt(mp_no *, mp_no *, int);
+
+void
+SECTION
+__mpatan(mp_no *x, mp_no *y, int p) {
+
int i,m,n;
double dx;
mp_no
@@ -54,19 +61,19 @@ void __mpatan(mp_no *x, mp_no *y, int p) {
mp_no mps,mpsm,mpt,mpt1,mpt2,mpt3;
- /* Choose m and initiate mpone, mptwo & mptwoim1 */
+ /* Choose m and initiate mpone, mptwo & mptwoim1 */
if (EX>0) m=7;
else if (EX<0) m=0;
else {
__mp_dbl(x,&dx,p); dx=ABS(dx);
for (m=6; m>0; m--)
- {if (dx>xm[m].d) break;}
+ {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;
- /* Reduce x m times */
+ /* Reduce x m times */
__mul(x,x,&mpsm,p);
if (m==0) __cpy(x,&mps,p);
else {
@@ -82,8 +89,8 @@ void __mpatan(mp_no *x, mp_no *y, int p) {
__mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0];
}
- /* Evaluate a truncated power series for Atan(s) */
- n=np[p]; mptwoim1.d[1] = twonm1[p].d;
+ /* Evaluate a truncated power series for Atan(s) */
+ 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;
@@ -94,8 +101,8 @@ void __mpatan(mp_no *x, mp_no *y, int p) {
__mul(&mps,&mpt,&mpt1,p);
__sub(&mps,&mpt1,&mpt,p);
- /* Compute Atan(x) */
- mptwoim1.d[1] = twom[m].d;
+ /* Compute Atan(x) */
+ mptwoim1.d[1] = __atan_twom[m].d;
__mul(&mptwoim1,&mpt,y,p);
return;
diff --git a/sysdeps/ieee754/dbl-64/mpatan.h b/sysdeps/ieee754/dbl-64/mpatan.h
index d420ff3408..003b06c695 100644
--- a/sysdeps/ieee754/dbl-64/mpatan.h
+++ b/sysdeps/ieee754/dbl-64/mpatan.h
@@ -1,8 +1,7 @@
-
/*
* IBM Accurate Mathematical Library
* Written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation, Inc.
+ * Copyright (C) 2001, 2011 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
@@ -29,9 +28,18 @@
#ifndef MPATAN_H
#define MPATAN_H
+extern const number __atan_xm[8] attribute_hidden;
+extern const number __atan_twonm1[33] attribute_hidden;
+extern const number __atan_twom[8] attribute_hidden;
+extern const number __atan_one attribute_hidden;
+extern const number __atan_two attribute_hidden;
+extern const int __atan_np[33] attribute_hidden;
+
+
+#ifndef AVOID_MPATAN_H
#ifdef BIG_ENDI
- static const number
- xm[8] = { /* x[m] */
+ const number
+ __atan_xm[8] = { /* x[m] */
/**/ {{0x00000000, 0x00000000} }, /* 0.0 */
/**/ {{0x3f8930be, 0x00000000} }, /* 0.0123 */
/**/ {{0x3f991687, 0x00000000} }, /* 0.0245 */
@@ -40,9 +48,9 @@
/**/ {{0x3fc95810, 0x00000000} }, /* 0.198 */
/**/ {{0x3fda7ef9, 0x00000000} }, /* 0.414 */
/**/ {{0x3ff00000, 0x00000000} }, /* 1.0 */
- };
- static const number
- twonm1[33] = { /* 2n-1 */
+ };
+ const number
+ __atan_twonm1[33] = { /* 2n-1 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
@@ -76,10 +84,10 @@
/**/ {{0x405b4000, 0x00000000} }, /* 109 */
/**/ {{0x405c4000, 0x00000000} }, /* 113 */
/**/ {{0x405d4000, 0x00000000} }, /* 117 */
- };
+ };
- static const number
- twom[8] = { /* 2**m */
+ const number
+ __atan_twom[8] = { /* 2**m */
/**/ {{0x3ff00000, 0x00000000} }, /* 1.0 */
/**/ {{0x40000000, 0x00000000} }, /* 2.0 */
/**/ {{0x40100000, 0x00000000} }, /* 4.0 */
@@ -88,17 +96,17 @@
/**/ {{0x40400000, 0x00000000} }, /* 32.0 */
/**/ {{0x40500000, 0x00000000} }, /* 64.0 */
/**/ {{0x40600000, 0x00000000} }, /* 128.0 */
- };
+ };
- static const number
-/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */
-/**/ two = {{0x40000000, 0x00000000} }; /* 2 */
+ const number
+/**/ __atan_one = {{0x3ff00000, 0x00000000} }, /* 1 */
+/**/ __atan_two = {{0x40000000, 0x00000000} }; /* 2 */
#else
#ifdef LITTLE_ENDI
- static const number
- xm[8] = { /* x[m] */
+ const number
+ __atan_xm[8] = { /* x[m] */
/**/ {{0x00000000, 0x00000000} }, /* 0.0 */
/**/ {{0x00000000, 0x3f8930be} }, /* 0.0123 */
/**/ {{0x00000000, 0x3f991687} }, /* 0.0245 */
@@ -107,9 +115,9 @@
/**/ {{0x00000000, 0x3fc95810} }, /* 0.198 */
/**/ {{0x00000000, 0x3fda7ef9} }, /* 0.414 */
/**/ {{0x00000000, 0x3ff00000} }, /* 1.0 */
- };
- static const number
- twonm1[33] = { /* 2n-1 */
+ };
+ const number
+__atan_twonm1[33] = { /* 2n-1 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
@@ -143,10 +151,10 @@
/**/ {{0x00000000, 0x405b4000} }, /* 109 */
/**/ {{0x00000000, 0x405c4000} }, /* 113 */
/**/ {{0x00000000, 0x405d4000} }, /* 117 */
- };
+ };
- static const number
- twom[8] = { /* 2**m */
+ const number
+ __atan_twom[8] = { /* 2**m */
/**/ {{0x00000000, 0x3ff00000} }, /* 1.0 */
/**/ {{0x00000000, 0x40000000} }, /* 2.0 */
/**/ {{0x00000000, 0x40100000} }, /* 4.0 */
@@ -155,20 +163,21 @@
/**/ {{0x00000000, 0x40400000} }, /* 32.0 */
/**/ {{0x00000000, 0x40500000} }, /* 64.0 */
/**/ {{0x00000000, 0x40600000} }, /* 128.0 */
- };
+ };
- static const number
-/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */
-/**/ two = {{0x00000000, 0x40000000} }; /* 2 */
+ const number
+/**/ __atan_one = {{0x00000000, 0x3ff00000} }, /* 1 */
+/**/ __atan_two = {{0x00000000, 0x40000000} }; /* 2 */
#endif
#endif
-#define ONE one.d
-#define TWO two.d
-
- static const int
- np[33] = { 0, 0, 0, 0, 6, 8,10,11,13,15,17,19,21,23,25,27,28,
- 30,32,34,36,38,40,42,43,45,47,49,51,53,55,57,59};
+ const int
+ __atan_np[33] = { 0, 0, 0, 0, 6, 8,10,11,13,15,17,19,21,23,25,27,28,
+ 30,32,34,36,38,40,42,43,45,47,49,51,53,55,57,59};
#endif
+#endif
+
+#define ONE __atan_one.d
+#define TWO __atan_two.d
diff --git a/sysdeps/ieee754/dbl-64/mpatan2.c b/sysdeps/ieee754/dbl-64/mpatan2.c
index 8977ec9042..1deb056417 100644
--- a/sysdeps/ieee754/dbl-64/mpatan2.c
+++ b/sysdeps/ieee754/dbl-64/mpatan2.c
@@ -1,8 +1,7 @@
-
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
*
* 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
@@ -38,18 +37,24 @@
#include "mpa.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
void __mpsqrt(mp_no *, mp_no *, int);
void __mpatan(mp_no *, mp_no *, int);
/* Multi-Precision Atan2(y,x) function subroutine, for p >= 4. */
/* y=0 is not permitted if x<=0. No error messages are given. */
-void __mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) {
+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}};
+ 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;
diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c
index e2ab71b2cc..b0cffe2fe5 100644
--- a/sysdeps/ieee754/dbl-64/mpexp.c
+++ b/sysdeps/ieee754/dbl-64/mpexp.c
@@ -1,8 +1,7 @@
-
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
*
* 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
@@ -34,34 +33,40 @@
#include "mpa.h"
#include "mpexp.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
/* Multi-Precision exponential function subroutine (for p >= 4, */
/* 2**(-55) <= abs(x) <= 1024). */
-void __mpexp(mp_no *x, mp_no *y, int p) {
+void
+SECTION
+__mpexp(mp_no *x, mp_no *y, int p) {
int i,j,k,m,m1,m2,n;
double a,b;
static const int np[33] = {0,0,0,0,3,3,4,4,5,4,4,5,5,5,6,6,6,6,6,6,
- 6,6,6,6,7,7,7,7,8,8,8,8,8};
+ 6,6,6,6,7,7,7,7,8,8,8,8,8};
static const int m1p[33]= {0,0,0,0,17,23,23,28,27,38,42,39,43,47,43,47,50,54,
- 57,60,64,67,71,74,68,71,74,77,70,73,76,78,81};
+ 57,60,64,67,71,74,68,71,74,77,70,73,76,78,81};
static const int m1np[7][18] = {
- { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
- { 0, 0, 0, 0,36,48,60,72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
- { 0, 0, 0, 0,24,32,40,48,56,64,72, 0, 0, 0, 0, 0, 0, 0},
- { 0, 0, 0, 0,17,23,29,35,41,47,53,59,65, 0, 0, 0, 0, 0},
- { 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0},
- { 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63},
- { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}};
+ { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ { 0, 0, 0, 0,36,48,60,72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ { 0, 0, 0, 0,24,32,40,48,56,64,72, 0, 0, 0, 0, 0, 0, 0},
+ { 0, 0, 0, 0,17,23,29,35,41,47,53,59,65, 0, 0, 0, 0, 0},
+ { 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0},
+ { 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63},
+ { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}};
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}};
+ 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 mpk = {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}};
+ 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 mps,mpak,mpt1,mpt2;
/* Choose m,n and compute a=2**(-m) */
- n = np[p]; m1 = m1p[p]; a = twomm1[p].d;
+ n = np[p]; m1 = m1p[p]; a = __mpexp_twomm1[p].d;
for (i=0; i<EX; i++) a *= RADIXI;
for ( ; i>EX; i--) a *= RADIX;
b = X[1]*RADIXI; m2 = 24*EX;
@@ -81,12 +86,12 @@ void __mpexp(mp_no *x, mp_no *y, int 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]=nn[n].d;
+ mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=__mpexp_nn[n].d;
__dvd(&mps,&mpk,&mpt1,p);
__add(&mpone,&mpt1,&mpak,p);
for (k=n-1; k>1; k--) {
__mul(&mps,&mpak,&mpt1,p);
- mpk.d[1]=nn[k].d;
+ mpk.d[1]=__mpexp_nn[k].d;
__dvd(&mpt1,&mpk,&mpt2,p);
__add(&mpone,&mpt2,&mpak,p);
}
diff --git a/sysdeps/ieee754/dbl-64/mpexp.h b/sysdeps/ieee754/dbl-64/mpexp.h
index a0b08ccb41..7985060a8c 100644
--- a/sysdeps/ieee754/dbl-64/mpexp.h
+++ b/sysdeps/ieee754/dbl-64/mpexp.h
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* Written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation, Inc.
+ * Copyright (C) 2001, 2011 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
@@ -28,9 +28,20 @@
#ifndef MPEXP_H
#define MPEXP_H
+extern const number __mpexp_twomm1[33] attribute_hidden;
+extern const number __mpexp_nn[9] attribute_hidden;
+extern const number __mpexp_radix attribute_hidden;
+extern const number __mpexp_radixi attribute_hidden;
+extern const number __mpexp_zero attribute_hidden;
+extern const number __mpexp_one attribute_hidden;
+extern const number __mpexp_two attribute_hidden;
+extern const number __mpexp_half attribute_hidden;
+
+
+#ifndef AVOID_MPEXP_H
#ifdef BIG_ENDI
- static const number
- twomm1[33] = { /* 2**-m1 */
+ const number
+ __mpexp_twomm1[33] = { /* 2**-m1 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
@@ -65,8 +76,8 @@
/**/ {{0x3b100000, 0x00000000} }, /* 2**-78 */
/**/ {{0x3ae00000, 0x00000000} }, /* 2**-81 */
};
- static const number
- nn[9]={ /* n */
+ const number
+ __mpexp_nn[9]={ /* n */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
/**/ {{0x3ff00000, 0x00000000} }, /* 1 */
/**/ {{0x40000000, 0x00000000} }, /* 2 */
@@ -78,18 +89,18 @@
/**/ {{0x40200000, 0x00000000} }, /* 8 */
};
- static const number
-/**/ radix = {{0x41700000, 0x00000000} }, /* 2**24 */
-/**/ radixi = {{0x3e700000, 0x00000000} }, /* 2**-24 */
-/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */
-/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */
-/**/ two = {{0x40000000, 0x00000000} }, /* 2 */
-/**/ half = {{0x3fe00000, 0x00000000} }; /* 1/2 */
+ const number
+/**/ __mpexp_radix = {{0x41700000, 0x00000000} }, /* 2**24 */
+/**/ __mpexp_radixi = {{0x3e700000, 0x00000000} }, /* 2**-24 */
+/**/ __mpexp_zero = {{0x00000000, 0x00000000} }, /* 0 */
+/**/ __mpexp_one = {{0x3ff00000, 0x00000000} }, /* 1 */
+/**/ __mpexp_two = {{0x40000000, 0x00000000} }, /* 2 */
+/**/ __mpexp_half = {{0x3fe00000, 0x00000000} }; /* 1/2 */
#else
#ifdef LITTLE_ENDI
- static const number
- twomm1[33] = { /* 2**-m1 */
+ const number
+ __mpexp_twomm1[33] = { /* 2**-m1 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
@@ -124,8 +135,8 @@
/**/ {{0x00000000, 0x3b100000} }, /* 2**-78 */
/**/ {{0x00000000, 0x3ae00000} }, /* 2**-81 */
};
- static const number
- nn[9]={ /* n */
+ const number
+ __mpexp_nn[9]={ /* n */
/**/ {{0x00000000, 0x00000000} }, /* 0 */
/**/ {{0x00000000, 0x3ff00000} }, /* 1 */
/**/ {{0x00000000, 0x40000000} }, /* 2 */
@@ -137,22 +148,23 @@
/**/ {{0x00000000, 0x40200000} }, /* 8 */
};
- static const number
-/**/ radix = {{0x00000000, 0x41700000} }, /* 2**24 */
-/**/ radixi = {{0x00000000, 0x3e700000} }, /* 2**-24 */
-/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */
-/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */
-/**/ two = {{0x00000000, 0x40000000} }, /* 2 */
-/**/ half = {{0x00000000, 0x3fe00000} }; /* 1/2 */
+ const number
+/**/ __mpexp_radix = {{0x00000000, 0x41700000} }, /* 2**24 */
+/**/ __mpexp_radixi = {{0x00000000, 0x3e700000} }, /* 2**-24 */
+/**/ __mpexp_zero = {{0x00000000, 0x00000000} }, /* 0 */
+/**/ __mpexp_one = {{0x00000000, 0x3ff00000} }, /* 1 */
+/**/ __mpexp_two = {{0x00000000, 0x40000000} }, /* 2 */
+/**/ __mpexp_half = {{0x00000000, 0x3fe00000} }; /* 1/2 */
#endif
#endif
+#endif
-#define RADIX radix.d
-#define RADIXI radixi.d
-#define ZERO zero.d
-#define ONE one.d
-#define TWO two.d
-#define HALF half.d
+#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
diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c
index bea623296b..d1a80f9091 100644
--- a/sysdeps/ieee754/dbl-64/mpsqrt.c
+++ b/sysdeps/ieee754/dbl-64/mpsqrt.c
@@ -32,6 +32,12 @@
#include "endian.h"
#include "mpa.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
+#include "mpsqrt.h"
+
/****************************************************************************/
/* Multi-Precision square root function subroutine for precision p >= 4. */
/* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */
@@ -42,9 +48,9 @@
static double fastiroot(double);
-void __mpsqrt(mp_no *x, mp_no *y, int p) {
-#include "mpsqrt.h"
-
+void
+SECTION
+__mpsqrt(mp_no *x, mp_no *y, int p) {
int i,m,ex,ey;
double dx,dy;
mp_no
@@ -64,7 +70,7 @@ void __mpsqrt(mp_no *x, mp_no *y, int p) {
__mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p);
__mul(&mpxn,&mphalf,&mpz,p);
- m=mp[p];
+ m=__mpsqrt_mp[p];
for (i=0; i<m; i++) {
__mul(&mpu,&mpu,&mpt1,p);
__mul(&mpt1,&mpz,&mpt2,p);
@@ -81,7 +87,9 @@ void __mpsqrt(mp_no *x, mp_no *y, int p) {
/* Compute a double precision approximation for 1/sqrt(x) */
/* with the relative error bounded by 2**-51. */
/***********************************************************/
-static double fastiroot(double x) {
+static double
+SECTION
+fastiroot(double x) {
union {int i[2]; double d;} p,q;
double y,z, t;
int n;
diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.h b/sysdeps/ieee754/dbl-64/mpsqrt.h
index 729d57af2c..86fa397b22 100644
--- a/sysdeps/ieee754/dbl-64/mpsqrt.h
+++ b/sysdeps/ieee754/dbl-64/mpsqrt.h
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* Written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation, Inc.
+ * Copyright (C) 2001, 2011 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
@@ -28,24 +28,31 @@
#ifndef MPSQRT_H
#define MPSQRT_H
+extern const number __mpsqrt_one attribute_hidden;
+extern const number __mpsqrt_halfrad attribute_hidden;
+extern const int __mpsqrt_mp[33] attribute_hidden;
+
+
+#ifndef AVOID_MPSQRT_H
#ifdef BIG_ENDI
- static const number
-/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */
-/**/ halfrad = {{0x41600000, 0x00000000} }; /* 2**23 */
+ const number
+/**/ __mpsqrt_one = {{0x3ff00000, 0x00000000} }, /* 1 */
+/**/ __mpsqrt_halfrad = {{0x41600000, 0x00000000} }; /* 2**23 */
#else
#ifdef LITTLE_ENDI
- static const number
-/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */
-/**/ halfrad = {{0x00000000, 0x41600000} }; /* 2**23 */
+ const number
+/**/ __mpsqrt_one = {{0x00000000, 0x3ff00000} }, /* 1 */
+/**/ __mpsqrt_halfrad = {{0x00000000, 0x41600000} }; /* 2**23 */
#endif
#endif
-#define ONE one.d
-#define HALFRAD halfrad.d
+ const int __mpsqrt_mp[33] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,
+ 4,4,4,4,4,4,4,4,4};
+#endif
- static const int mp[33] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,
- 4,4,4,4,4,4,4,4,4};
+#define ONE __mpsqrt_one.d
+#define HALFRAD __mpsqrt_halfrad.d
#endif
diff --git a/sysdeps/ieee754/dbl-64/mptan.c b/sysdeps/ieee754/dbl-64/mptan.c
index 267445a19e..e1e5d9b925 100644
--- a/sysdeps/ieee754/dbl-64/mptan.c
+++ b/sysdeps/ieee754/dbl-64/mptan.c
@@ -1,8 +1,7 @@
-
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
*
* 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
@@ -38,10 +37,16 @@
#include "endian.h"
#include "mpa.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
int __mpranred(double, mp_no *, int);
void __c32(mp_no *, mp_no *, mp_no *, int);
-void __mptan(double x, mp_no *mpy, int p) {
+void
+SECTION
+__mptan(double x, mp_no *mpy, int p) {
static const double MONE = -1.0;
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index 02d428ca03..6f19f158f1 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -55,6 +55,10 @@
#include "MathLib.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
extern const union
{
int4 i[880];
@@ -92,7 +96,9 @@ static double csloww2(double x, double dx, double orig, int n);
/* An ultimate sin routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of sin(x) */
/*******************************************************************/
-double __sin(double x){
+double
+SECTION
+__sin(double x){
double xx,res,t,cor,y,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2;
#if 0
double w[2];
@@ -349,7 +355,9 @@ double __sin(double x){
/* it computes the correctly rounded (to nearest) value of cos(x) */
/*******************************************************************/
-double __cos(double x)
+double
+SECTION
+__cos(double x)
{
double y,xx,res,t,cor,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2;
mynumber u,v;
@@ -596,7 +604,9 @@ double __cos(double x)
/* precision and if still doesn't accurate enough by mpsin or dubsin */
/************************************************************************/
-static double slow(double x) {
+static double
+SECTION
+slow(double x) {
static const double th2_36 = 206158430208.0; /* 1.5*2**37 */
double y,x1,x2,xx,r,t,res,cor,w[2];
x1=(x+th2_36)-th2_36;
@@ -620,7 +630,9 @@ static const double th2_36 = 206158430208.0; /* 1.5*2**37 */
/* and if result still doesn't accurate enough by mpsin or dubsin */
/*******************************************************************************/
-static double slow1(double x) {
+static double
+SECTION
+slow1(double x) {
mynumber u;
double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res;
static const double t22 = 6291456.0;
@@ -656,7 +668,9 @@ static double slow1(double x) {
/* Routine compute sin(x) for 0.855469 <|x|<2.426265 by __sincostab.tbl */
/* and if result still doesn't accurate enough by mpsin or dubsin */
/**************************************************************************/
-static double slow2(double x) {
+static double
+SECTION
+slow2(double x) {
mynumber u;
double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res,del;
static const double t22 = 6291456.0;
@@ -708,7 +722,9 @@ static double slow2(double x) {
/* result.And if result not accurate enough routine calls mpsin1 or dubsin */
/***************************************************************************/
-static double sloww(double x,double dx, double orig) {
+static double
+SECTION
+sloww(double x,double dx, double orig) {
static const double th2_36 = 206158430208.0; /* 1.5*2**37 */
double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn;
union {int4 i[2]; double x;} v;
@@ -755,7 +771,9 @@ static double sloww(double x,double dx, double orig) {
/* accurate enough routine calls mpsin1 or dubsin */
/***************************************************************************/
-static double sloww1(double x, double dx, double orig) {
+static double
+SECTION
+sloww1(double x, double dx, double orig) {
mynumber u;
double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res;
static const double t22 = 6291456.0;
@@ -797,7 +815,9 @@ static double sloww1(double x, double dx, double orig) {
/* accurate enough routine calls mpsin1 or dubsin */
/***************************************************************************/
-static double sloww2(double x, double dx, double orig, int n) {
+static double
+SECTION
+sloww2(double x, double dx, double orig, int n) {
mynumber u;
double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res;
static const double t22 = 6291456.0;
@@ -841,7 +861,9 @@ static double sloww2(double x, double dx, double orig, int n) {
/* result.And if result not accurate enough routine calls other routines */
/***************************************************************************/
-static double bsloww(double x,double dx, double orig,int n) {
+static double
+SECTION
+bsloww(double x,double dx, double orig,int n) {
static const double th2_36 = 206158430208.0; /* 1.5*2**37 */
double y,x1,x2,xx,r,t,res,cor,w[2];
#if 0
@@ -874,7 +896,9 @@ static double bsloww(double x,double dx, double orig,int n) {
/* And if result not accurate enough routine calls other routines */
/***************************************************************************/
-static double bsloww1(double x, double dx, double orig,int n) {
+static double
+SECTION
+bsloww1(double x, double dx, double orig,int n) {
mynumber u;
double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res;
static const double t22 = 6291456.0;
@@ -917,7 +941,9 @@ mynumber u;
/* And if result not accurate enough routine calls other routines */
/***************************************************************************/
-static double bsloww2(double x, double dx, double orig, int n) {
+static double
+SECTION
+bsloww2(double x, double dx, double orig, int n) {
mynumber u;
double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res;
static const double t22 = 6291456.0;
@@ -959,7 +985,9 @@ mynumber u;
/* precision and if still doesn't accurate enough by mpcos or docos */
/************************************************************************/
-static double cslow2(double x) {
+static double
+SECTION
+cslow2(double x) {
mynumber u;
double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res;
static const double t22 = 6291456.0;
@@ -1002,7 +1030,9 @@ static double cslow2(double x) {
/***************************************************************************/
-static double csloww(double x,double dx, double orig) {
+static double
+SECTION
+csloww(double x,double dx, double orig) {
static const double th2_36 = 206158430208.0; /* 1.5*2**37 */
double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn;
union {int4 i[2]; double x;} v;
@@ -1051,7 +1081,9 @@ static double csloww(double x,double dx, double orig) {
/* accurate enough routine calls other routines */
/***************************************************************************/
-static double csloww1(double x, double dx, double orig) {
+static double
+SECTION
+csloww1(double x, double dx, double orig) {
mynumber u;
double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res;
static const double t22 = 6291456.0;
@@ -1095,7 +1127,9 @@ static double csloww1(double x, double dx, double orig) {
/* accurate enough routine calls other routines */
/***************************************************************************/
-static double csloww2(double x, double dx, double orig, int n) {
+static double
+SECTION
+csloww2(double x, double dx, double orig, int n) {
mynumber u;
double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res;
static const double t22 = 6291456.0;
diff --git a/sysdeps/ieee754/dbl-64/s_tan.c b/sysdeps/ieee754/dbl-64/s_tan.c
index 4089e0d65f..f0fcd677ae 100644
--- a/sysdeps/ieee754/dbl-64/s_tan.c
+++ b/sysdeps/ieee754/dbl-64/s_tan.c
@@ -41,10 +41,16 @@
#include "MathLib.h"
#include "math.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
static double tanMp(double);
void __mptan(double, mp_no *, int);
-double tan(double x) {
+double
+SECTION
+tan(double x) {
#include "utan.h"
#include "utan.tbl"
@@ -486,7 +492,9 @@ double tan(double x) {
/* multiple precision stage */
/* Convert x to multi precision number,compute tan(x) by mptan() routine */
/* and converts result back to double */
-static double tanMp(double x)
+static double
+SECTION
+tanMp(double x)
{
int p;
double y;
diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c
index a4f896a465..e39aaeea06 100644
--- a/sysdeps/ieee754/dbl-64/sincos32.c
+++ b/sysdeps/ieee754/dbl-64/sincos32.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
*
* 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
@@ -45,11 +45,17 @@
#include "sincos32.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
/****************************************************************/
/* Compute Multi-Precision sin() function for given p. Receive */
/* Multi Precision number x and result stored at y */
/****************************************************************/
-static void ss32(mp_no *x, mp_no *y, int p) {
+static void
+SECTION
+ss32(mp_no *x, mp_no *y, int p) {
int i;
double a;
#if 0
@@ -79,7 +85,9 @@ static void ss32(mp_no *x, mp_no *y, int p) {
/* Compute Multi-Precision cos() function for given p. Receive Multi */
/* Precision number x and result stored at y */
/**********************************************************************/
-static void cc32(mp_no *x, mp_no *y, int p) {
+static void
+SECTION
+cc32(mp_no *x, mp_no *y, int p) {
int i;
double a;
#if 0
@@ -109,7 +117,9 @@ static void cc32(mp_no *x, mp_no *y, int p) {
/***************************************************************************/
/* c32() computes both sin(x), cos(x) as Multi precision numbers */
/***************************************************************************/
-void __c32(mp_no *x, mp_no *y, mp_no *z, int p) {
+void
+SECTION
+__c32(mp_no *x, mp_no *y, mp_no *z, int p) {
static const mp_no mpt={1,{1.0,2.0}}, one={1,{1.0,1.0}};
mp_no u,t,t1,t2,c,s;
int i;
@@ -134,7 +144,9 @@ void __c32(mp_no *x, mp_no *y, mp_no *z, int p) {
/*result which is more accurate */
/*Computing sin(x) with multi precision routine c32 */
/************************************************************************/
-double __sin32(double x, double res, double res1) {
+double
+SECTION
+__sin32(double x, double res, double res1) {
int p;
mp_no a,b,c;
p=32;
@@ -158,7 +170,9 @@ double __sin32(double x, double res, double res1) {
/*result which is more accurate */
/*Computing cos(x) with multi precision routine c32 */
/************************************************************************/
-double __cos32(double x, double res, double res1) {
+double
+SECTION
+__cos32(double x, double res, double res1) {
int p;
mp_no a,b,c;
p=32;
@@ -172,12 +186,12 @@ double __cos32(double x, double res, double res1) {
}
else if (x>0.8)
{ __sub(&hp,&c,&a,p);
- __c32(&a,&c,&b,p);
+ __c32(&a,&c,&b,p);
}
else __c32(&c,&b,&a,p); /* b=cos(0.5*(res+res1)) */
__dbl_mp(x,&c,p); /* c = x */
__sub(&b,&c,&a,p);
- /* if a>0 return max(res,res1), otherwise return min(res,res1) */
+ /* if a>0 return max(res,res1), otherwise return min(res,res1) */
if (a.d[0]>0) return (res>res1)?res:res1;
else return (res<res1)?res:res1;
}
@@ -186,7 +200,9 @@ double __cos32(double x, double res, double res1) {
/*Compute sin(x+dx) as Multi Precision number and return result as */
/* double */
/*******************************************************************/
-double __mpsin(double x, double dx) {
+double
+SECTION
+__mpsin(double x, double dx) {
int p;
double y;
mp_no a,b,c;
@@ -204,7 +220,9 @@ double __mpsin(double x, double dx) {
/* Compute cos()of double-length number (x+dx) as Multi Precision */
/* number and return result as double */
/*******************************************************************/
-double __mpcos(double x, double dx) {
+double
+SECTION
+__mpcos(double x, double dx) {
int p;
double y;
mp_no a,b,c;
@@ -227,7 +245,9 @@ double __mpcos(double x, double dx) {
/* n=0,+-1,+-2,.... */
/* Return int which indicates in which quarter of circle x is */
/******************************************************************/
-int __mpranred(double x, mp_no *y, int p)
+int
+SECTION
+__mpranred(double x, mp_no *y, int p)
{
number v;
double t,xn;
@@ -275,7 +295,9 @@ int __mpranred(double x, mp_no *y, int p)
/* Multi-Precision sin() function subroutine, for p=32. It is */
/* based on the routines mpranred() and c32(). */
/*******************************************************************/
-double __mpsin1(double x)
+double
+SECTION
+__mpsin1(double x)
{
int p;
int n;
@@ -314,7 +336,9 @@ double __mpsin1(double x)
/* based on the routines mpranred() and c32(). */
/*****************************************************************/
-double __mpcos1(double x)
+double
+SECTION
+__mpcos1(double x)
{
int p;
int n;
diff --git a/sysdeps/ieee754/dbl-64/slowexp.c b/sysdeps/ieee754/dbl-64/slowexp.c
index 78c107f709..6a6bce31ba 100644
--- a/sysdeps/ieee754/dbl-64/slowexp.c
+++ b/sysdeps/ieee754/dbl-64/slowexp.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
*
* 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
@@ -31,10 +31,16 @@
#include "mpa.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
void __mpexp(mp_no *x, mp_no *y, int p);
/*Converting from double precision to Multi-precision and calculating e^x */
-double __slowexp(double x) {
+double
+SECTION
+__slowexp(double x) {
double w,z,res,eps=3.0e-26;
#if 0
double y;
@@ -47,7 +53,7 @@ double __slowexp(double x) {
p=6;
__dbl_mp(x,&mpx,p); /* Convert a double precision number x */
- /* into a multiple precision number mpx with prec. p. */
+ /* into a multiple precision number mpx with prec. p. */
__mpexp(&mpx, &mpy, p); /* Multi-Precision exponential function */
__dbl_mp(eps,&mpeps,p);
__mul(&mpeps,&mpy,&mpcor,p);
diff --git a/sysdeps/ieee754/dbl-64/slowpow.c b/sysdeps/ieee754/dbl-64/slowpow.c
index e11a532bf8..0c57e6d4f2 100644
--- a/sysdeps/ieee754/dbl-64/slowpow.c
+++ b/sysdeps/ieee754/dbl-64/slowpow.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
*
* 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
@@ -35,12 +35,18 @@
#include "mpa.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
void __mpexp(mp_no *x, mp_no *y, int p);
void __mplog(mp_no *x, mp_no *y, int p);
double ulog(double);
double __halfulp(double x,double y);
-double __slowpow(double x, double y, double z) {
+double
+SECTION
+__slowpow(double x, double y, double z) {
double res,res1;
mp_no mpx, mpy, mpz,mpw,mpp,mpr,mpr1;
static const mp_no eps = {-3,{1.0,4.0}};
@@ -48,7 +54,7 @@ double __slowpow(double x, double y, double z) {
res = __halfulp(x,y); /* halfulp() returns -10 or x^y */
if (res >= 0) return res; /* if result was really computed by halfulp */
- /* else, if result was not really computed by halfulp */
+ /* else, if result was not really computed by halfulp */
p = 10; /* p=precision */
__dbl_mp(x,&mpx,p);
__dbl_mp(y,&mpy,p);