+2011-10-25 Ulrich Drepper <drepper@gmail.com>
+
+ * sysdeps/ieee754/dbl-64/branred.c: Move FMA4 code into separate
+ .text section. Avoid duplicate constants.
+ * sysdeps/ieee754/dbl-64/doasin.c: Likewise.
+ * sysdeps/ieee754/dbl-64/dosincos.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_asin.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_atan2.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_exp.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_log.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_pow.c: Likewise.
+ * sysdeps/ieee754/dbl-64/halfulp.c: Likewise.
+ * sysdeps/ieee754/dbl-64/mpa.c: Likewise.
+ * sysdeps/ieee754/dbl-64/mpa.h: Likewise.
+ * sysdeps/ieee754/dbl-64/mpatan.c: Likewise.
+ * sysdeps/ieee754/dbl-64/mpatan.h: Likewise.
+ * sysdeps/ieee754/dbl-64/mpatan2.c: Likewise.
+ * sysdeps/ieee754/dbl-64/mpexp.c: Likewise.
+ * sysdeps/ieee754/dbl-64/mpexp.h: Likewise.
+ * sysdeps/ieee754/dbl-64/mpsqrt.c: Likewise.
+ * sysdeps/ieee754/dbl-64/mpsqrt.h: Likewise.
+ * sysdeps/ieee754/dbl-64/mptan.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_sin.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_tan.c: Likewise.
+ * sysdeps/ieee754/dbl-64/sincos32.c: Likewise.
+ * sysdeps/ieee754/dbl-64/slowexp.c: Likewise.
+ * sysdeps/ieee754/dbl-64/slowpow.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/brandred-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/doasin-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/dosincos-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/e_asin-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/e_atan2-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/e_log-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/e_pow-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/halfulp-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/mpa-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/mpatan-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/mpatan2-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/mpexp-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/mplog-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/mpsqrt-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/mptan-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/s_atan-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/s_sin-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/s_tan-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/sincos32-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/slowexp-fma4.c: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/slowpow-fma4.c: Likewise.
+
2011-10-24 Ulrich Drepper <drepper@gmail.com>
* sysdeps/x86_64/dla.h: Move to ...
/*
* 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
#include "branred.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
/*******************************************************************/
/* Routine branred() performs range reduction of a double number */
/* 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
#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"
#include "dosincos.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
extern const union
{
int4 i[880];
/*(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
/*(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
/* 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;}
#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[]);
/* 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;
/* */
/*******************************************************************/
-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
#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). */
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
#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;
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;
#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
/*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
#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;
#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);
/* 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;
/**************************************************************************/
/* 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;
/* + 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;
/* 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;
/* 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
#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,
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;
#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;
}
-
/* 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) {
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;
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
for (i=1; i>EX; i--) c *= RADIXI;
*y = c;
- return;
#undef R
}
c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10);
*y = c*TWOM1032;
- return;
-
#undef R
}
/* 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;
if (u>x) u -= ONE;
Y[i] = u; x -= u; x *= RADIX; }
for ( ; i<=p; i++) Y[i] = ZERO;
- return;
}
/* 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;
/* 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;
Z[k++] = Z[i++];
for (; k <= p; )
Z[k++] = ZERO;
-
- return;
}
/* 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;
else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; }
else Z[0] = ZERO;
}
- return;
}
/* 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;
else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = -Y[0]; }
else Z[0] = ZERO;
}
- return;
}
/* 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;
EZ = EX + EY;
Z[0] = X[0] * Y[0];
- return;
}
/* 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;
__sub(&mptwo,y,&z,p);
__mul(&w,&z,y,p);
}
- return;
}
/* 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;
}
#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);
-
/*
* 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
#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
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 {
__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;
__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;
-
/*
* 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
#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 */
/**/ {{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 */
/**/ {{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 */
/**/ {{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 */
/**/ {{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 */
/**/ {{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 */
/**/ {{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
-
/*
* 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
#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;
-
/*
* 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
#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;
/* 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);
}
/*
* 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
#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 */
/**/ {{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 */
/**/ {{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 */
/**/ {{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 */
/**/ {{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
#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. */
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
__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);
/* 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;
/*
* 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
#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
-
/*
* 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
#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;
#include "MathLib.h"
#include "math_private.h"
+#ifndef SECTION
+# define SECTION
+#endif
+
extern const union
{
int4 i[880];
/* 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];
/* 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;
/* 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;
/* 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;
/* 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;
/* 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;
/* 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;
/* 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;
/* 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
/* 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;
/* 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;
/* 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;
/***************************************************************************/
-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;
/* 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;
/* 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;
#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"
/* 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;
/*
* 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
#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
/* 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
/***************************************************************************/
/* 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;
/*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;
/*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;
}
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;
}
/*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;
/* 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;
/* 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;
/* 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;
/* based on the routines mpranred() and c32(). */
/*****************************************************************/
-double __mpcos1(double x)
+double
+SECTION
+__mpcos1(double x)
{
int p;
int n;
/*
* 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
#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;
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);
/*
* 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
#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}};
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);
#define __branred __branred_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/branred.c>
#define __doasin __doasin_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/doasin.c>
#define __docos __docos_fma4
#define __dubcos __dubcos_fma4
#define __dubsin __dubsin_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/dosincos.c>
#define __dubcos __dubcos_fma4
#define __dubsin __dubsin_fma4
#define __sin32 __sin32_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/e_asin.c>
#define __mpatan2 __mpatan2_fma4
#define __mul __mul_fma4
#define __sub __sub_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/e_atan2.c>
#define __ieee754_exp __ieee754_exp_fma4
#define __exp1 __exp1_fma4
#define __slowexp __slowexp_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/e_exp.c>
#define __add __add_fma4
#define __dbl_mp __dbl_mp_fma4
#define __sub __sub_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/e_log.c>
#define __ieee754_pow __ieee754_pow_fma4
#define __exp1 __exp1_fma4
#define __slowpow __slowpow_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/e_pow.c>
#define __halfulp __halfulp_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/halfulp.c>
#define NO___CPY 1
#define NO___MP_DBL 1
+#define NO___ACR 1
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/mpa.c>
#define __mpsqrt __mpsqrt_fma4
#define __mul __mul_fma4
#define __sub __sub_fma4
+#define AVOID_MPATAN_H 1
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/mpatan.c>
#define __mpatan __mpatan_fma4
#define __mpsqrt __mpsqrt_fma4
#define __mul __mul_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/mpatan2.c>
#define __dbl_mp __dbl_mp_fma4
#define __dvd __dvd_fma4
#define __mul __mul_fma4
+#define AVOID_MPEXP_H 1
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/mpexp.c>
#define __mpexp __mpexp_fma4
#define __mul __mul_fma4
#define __sub __sub_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/mplog.c>
#define __dbl_mp __dbl_mp_fma4
#define __mul __mul_fma4
#define __sub __sub_fma4
+#define AVOID_MPSQRT_H 1
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/mpsqrt.c>
#define __c32 __c32_fma4
#define __dvd __dvd_fma4
#define __mpranred __mpranred_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/mptan.c>
#define __mpatan __mpatan_fma4
#define __mul __mul_fma4
#define __sub __sub_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/s_atan.c>
#define __mpcos1 __mpcos1_fma4
#define __mpsin __mpsin_fma4
#define __mpsin1 __mpsin1_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/s_sin.c>
#define __mpranred __mpranred_fma4
#define __mptan __mptan_fma4
#define __sub __sub_fma4
-
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/s_tan.c>
#define __dbl_mp __dbl_mp_fma4
#define __mul __mul_fma4
#define __sub __sub_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/sincos32.c>
#define __mpexp __mpexp_fma4
#define __mul __mul_fma4
#define __sub __sub_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/slowexp.c>
#define __mul __mul_fma4
#define __sub __sub_fma4
#define __halfulp __halfulp_fma4
+#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/slowpow.c>