/*
*this function computes log, correctly rounded,
using double-extended arithmetic
THIS IS EXPERIMENTAL SOFTWARE
In particular it changes rounding modes all the time without warning
nor restoring.
*
* Author : Florent de Dinechin
* Florent.de.Dinechin at ens-lyon.fr
*
This function compiles both on IA32 and IA64 architectures. On IA64,
it needs icc 8.1 or higher, with the following flags (which should be
set up by the autoconf).
icc -DHAVE_CONFIG_H -Qoption,cpp,--extended_float_types \
-IPF_fp_speculationsafe -c log-de.c;\
mv log-de.o log-td.o; make
*/
#include <stdio.h>
#include <stdlib.h>
#include "crlibm.h"
#include "crlibm_private.h"
#include "double-extended.h"
#include "log-de.h"
static void log_accurate(double_ext* prh, double_ext* prl, double_ext z, int E, int index) {
double_ext eh,el, t13, t12, t11, t10, t9, t8,
p7h,p7l, t7h,t7l, t6h,t6l, t5h,t5l, t4h,t4l,
t3h,t3l, t2h,t2l, t1h,t1l, t0h,t0l;
/* Many temporary because single assignment form is nicer for Gappa */
#if !(defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64))
double_ext c1h,c2h,c3h,c4h,c5h,c6h,c7h,c8h,c9h,c10h,c11h,c12h,c13h,c14h,c15h;
double_ext c1l,c2l,c3l,c4l,c5l,c6l,c7l,c8l;
#endif
#if EVAL_PERF
crlibm_second_step_taken++;
#endif
/* TODO check the conditions for the double-double ops */
PREFETCH_POLY_ACCURATE;
t13 = c13h + z*c14h;
t12 = c12h + z*t13;
t11 = c11h + z*t12;
t10 = c10h + z*t11;
t9 = c9h + z*t10;
t8 = c8h + z*t9;
#if 1 /* This is faster on PIII. Your mileage may vary */
Mul12_ext(&p7h, &p7l, z, t8);
Add22_ext(&t7h, &t7l, p7h,p7l, c7h,c7l);
#else
FMA22_ext(&t7h, &t7l, z,0, t8,0, c7h,c7l);
#endif
FMA22_ext(&t6h, &t6l, z,0, t7h,t7l, c6h,c6l);
FMA22_ext(&t5h, &t5l, z,0, t6h,t6l, c5h,c5l);
FMA22_ext(&t4h, &t4l, z,0, t5h,t5l, c4h,c4l);
FMA22_ext(&t3h, &t3l, z,0, t4h,t4l, c3h,c3l);
FMA22_ext(&t2h, &t2l, z,0, t3h,t3l, c2h,c2l);
FMA22_ext(&t1h, &t1l, z,0, t2h,t2l, c1h,c1l);
FMA22_ext(&t0h, &t0l, z,0, t1h,t1l, argredtable[index].logirh, argredtable[index].logirl);
Mul22_ext(&eh, &el, log2H,log2L, E, 0);
Add22_ext(prh, prl, eh,el, t0h,t0l);
}
double log_rn(double x) {
double_ext logirh, r, y, z, th, tl, logde;
#if defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64)
db_number xdb;
int E, index, index0, roundtestmask;
#else /* assuming Itanium here */
int64_t E, i;
uint64_t index, roundtestmask;
double c2,c3,c4,c5,c6,c7;
#endif
#if defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64)
xdb.d=x;
index0 = (xdb.i[HI] & 0x000fffff);
index = (index0 + (1<<(20-L-1))) >> (20-L);
E = (xdb.i[HI]>>20)-1023; /* extract the exponent */
/* Filter cases */
if (xdb.i[HI] < 0x00100000){ /* x < 2^(-1022) */
if (((xdb.i[HI] & 0x7fffffff)|xdb.i[LO])==0) return -1.0/0.0; /* log(+/-0) = -Inf */
if (xdb.i[HI] < 0) return (x-x)/0.0; /* log(-x) = Nan */
/* Else subnormal number */
xdb.d *= two64; /* make x a normal number */
E = -64 + (xdb.i[HI]>>20)-1023; /* extract the exponent */
index0 = (xdb.i[HI] & 0x000fffff);
index = (index0 + (1<<(20-L-1))) >> (20-L);
}
if (xdb.i[HI] >= 0x7ff00000) return x+x; /* Inf or Nan */
DOUBLE_EXTENDED_MODE; /* This one should be overlapped with following integer computation */
/* Extract exponent and mantissa */
xdb.i[HI] = index0 | 0x3ff00000; /* do exponent = 0 */
/* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
if (index >= MAXINDEX){ /* corresponds to y>sqrt(2)*/
xdb.i[HI] -= 0x00100000;
index = index & INDEXMASK;
E++;
}
y = xdb.d;
#else /* defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64) */
/* Here come the code specific to Itanium processor */
E=0;
PREFETCH_POLY_QUICK; /* defined in log-de.h */
y=x;
i = _Asm_getf(2/*_FR_D*/, y); /* Cast y to a 64-bit integer */
/* Filter special cases */
if (i<(int64_t)ULL(0010000000000000)){ /* equivalent to : x < 2^(-1022) */
if ((i & ULL(7fffffffffffffff))==0) return -1.0/0.0; /* log(+/-0) = -Inf */
if (i<0) return (x-x)/0.0; /* log(-x) = Nan */
/* Else subnormal number */
y *= two64; /* make x a normal number */
E = -64;
i = _Asm_getf(2/*_FR_D*/, y); /* and update i */
}
if (i >= ULL(7ff0000000000000)) return x+x; /* Inf or Nan */
/* Extract exponent and mantissa */
E += (i>>52)-1023;
i = i & ULL(000fffffffffffff); /* keep only mantissa */
index = (i + (ULL(1)<<(52-L-1))) >> (52-L);
/* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
if (index >= MAXINDEX){ /* corresponds to y>sqrt(2)*/
y = _Asm_setf(2/*_FR_D*/, (i | ULL(3ff0000000000000)) - ULL(0010000000000000) ); /* exponent = -1 */
index = index & INDEXMASK;
E++;
}
else
y = _Asm_setf(2/*_FR_D*/, i | ULL(3ff0000000000000) ); /* exponent = 0*/
#endif /* defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64) */
/* All the previous argument reduction was exact */
/* now y holds 1+f, and E is the exponent */
r = (double_ext) (argredtable[index].r); /* approx to 1/y.d */
logirh = argredtable[index].logirh;
z = y*r - 1. ; /* even without an FMA, all exact */
#if 0
if(E==0)
roundtestmask=ACCURATE_TO_61_BITS;
else
roundtestmask=ACCURATE_TO_61_BITS;
#else
roundtestmask=ACCURATE_TO_62_BITS;
#endif
#ifdef ESTRIN
/* Estrin polynomial evaluation */
double_ext z2,z4, p01, p23, p45, p67, p03, p47,p07;
z2 = z*z; p67 = c6 + z*c7; p45 = c4 + z*c5; p23 = c2 + z*c3; p01 = logirh + z;
z4 = z2*z2; p47 = p45 + z2*p67; p03 = p01 + z2*p23;
p07 = p03 + z4*p47;
logde = p07 + E*log2H;
#endif
#ifdef PATERSON
double_ext z4,z2,t0,t1,t2,t3,t4,t5,t6,t7,t8;
z2 = z * z; t1 = z + ps_alpha; t2 = z + ps_beta; t3 = c3 * z + c2; t4 = z + logirh;
z4 = z2 * z2; t5 = z2 + ps_c; t6 = t3 * z2 + t4;
t7 = t5 * t1 + t2; t0 = z4 * c7; t8 = t7 * t0 + t6;
logde = t8 + E*log2H;
#endif
#if 0 /* to time the first step only */
BACK_TO_DOUBLE_MODE; return (double)t;
#endif
/* To test the second step only, comment out the following line */
DE_TEST_AND_RETURN_RN(logde, roundtestmask);
log_accurate(&th, &tl, z, E, index);
BACK_TO_DOUBLE_MODE;
return (double) (th+tl); /* The exact sum of these double-extended is rounded to the nearest */
}
double log_rd(double x) {
double_ext logirh, r, y, z, th, tl, logde;
#if defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64)
db_number xdb;
int E, index, roundtestmask;
#else
int64_t E, i;
uint64_t index, roundtestmask;
double_ext c1,c2,c3,c4,c5,c6,c7;
#endif
E=0;
#if defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64)
xdb.d=x;
/* Filter cases */
if (xdb.i[HI] < 0x00100000){ /* x < 2^(-1022) */
if (((xdb.i[HI] & 0x7fffffff)|xdb.i[LO])==0) return -1.0/0.0; /* log(+/-0) = -Inf */
if (xdb.i[HI] < 0) return (x-x)/0.0; /* log(-x) = Nan */
/* Else subnormal number */
E = -64;
xdb.d *= two64; /* make x a normal number */
}
if (xdb.i[HI] >= 0x7ff00000) return x+x; /* Inf or Nan */
DOUBLE_EXTENDED_MODE; /* This one should be overlapped with following integer computation */
/* Extract exponent and mantissa */
E += (xdb.i[HI]>>20)-1023; /* extract the exponent */
index = (xdb.i[HI] & 0x000fffff);
xdb.i[HI] = index | 0x3ff00000; /* do exponent = 0 */
index = (index + (1<<(20-L-1))) >> (20-L);
/* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
if (index >= MAXINDEX){ /* corresponds to y>sqrt(2)*/
xdb.i[HI] -= 0x00100000;
E++;
}
y = xdb.d;
#else /* defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64) */
/* Here come the code specific to Itanium processor */
PREFETCH_POLY_QUICK; /* defined in log-de.h */
y=x;
i = _Asm_getf(2/*_FR_D*/, y); /* Cast y to a 64-bit integer */
/* Filter special cases */
if (i<(int64_t)ULL(0010000000000000)){ /* equivalent to : x < 2^(-1022) */
if ((i & ULL(7fffffffffffffff))==0) return -1.0/0.0; /* log(+/-0) = -Inf */
if (i<0) return (x-x)/0.0; /* log(-x) = Nan */
/* Else subnormal number */
y *= two64; /* make x a normal number */
E = -64;
i = _Asm_getf(2/*_FR_D*/, y); /* and update i */
}
if (i >= ULL(7ff0000000000000)) return x+x; /* Inf or Nan */
/* Extract exponent and mantissa */
E += (i>>52)-1023;
i = i & ULL(000fffffffffffff); /* keep only mantissa */
index = (i + (ULL(1)<<(52-L-1))) >> (52-L);
/* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
if (index >= MAXINDEX){ /* corresponds to y>sqrt(2)*/
y = _Asm_setf(2/*_FR_D*/, (i | ULL(3ff0000000000000)) - ULL(0010000000000000) ); /* exponent = -1 */
E++;
}
else
y = _Asm_setf(2/*_FR_D*/, i | ULL(3ff0000000000000) ); /* exponent = 0*/
#endif /* defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64) */
/* All the previous argument reduction was exact */
/* now y holds 1+f, and E is the exponent */
index = index & INDEXMASK;
logirh = argredtable[index].logirh;
r = (double_ext) (argredtable[index].r); /* approx to 1/y.d */
z = y*r - 1. ; /* even without an FMA, all exact */
if(E==0)
roundtestmask=ACCURATE_TO_61_BITS;
else
roundtestmask=ACCURATE_TO_62_BITS;
#ifdef ESTRIN
/* Estrin polynomial evaluation */
double_ext z2,z4, p01, p23, p45, p67, p03, p47,p07;
z2 = z*z; p67 = c6 + z*c7; p45 = c4 + z*c5; p23 = c2 + z*c3; p01 = logirh + z;
z4 = z2*z2; p47 = p45 + z2*p67; p03 = p01 + z2*p23;
p07 = p03 + z4*p47;
logde = p07 + E*log2H;
#endif
#ifdef PATERSON
double_ext z4,z2,t0,t1,t2,t3,t4,t5,t6,t7,t8;
z2 = z * z; t1 = z + ps_alpha; t2 = z + ps_beta; t3 = c3 * z + c2; t4 = z + logirh;
z4 = z2 * z2; t5 = z2 + ps_c; t6 = t3 * z2 + t4;
t7 = t5 * t1 + t2; t0 = z4 * c7; t8 = t7 * t0 + t6;
logde = t8 + E*log2H;
#endif
#if 0 /* to time the first step only */
BACK_TO_DOUBLE_MODE; return (double)t;
#endif
/* To test the second step only, comment out the following line */
DE_TEST_AND_RETURN_RD(logde, roundtestmask);
log_accurate(&th, &tl, z, E, index);
RETURN_SUM_ROUNDED_DOWN(th, tl);
}
double log_ru(double x) {
double_ext logirh, r, y, z, th, tl, logde;
#if defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64)
db_number xdb;
int E, index, roundtestmask;
#else
int64_t E, i;
uint64_t index, roundtestmask;
double_ext c1,c2,c3,c4,c5,c6,c7;
#endif
E=0;
#if defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64)
xdb.d=x;
/* Filter cases */
if (xdb.i[HI] < 0x00100000){ /* x < 2^(-1022) */
if (((xdb.i[HI] & 0x7fffffff)|xdb.i[LO])==0) return -1.0/0.0; /* log(+/-0) = -Inf */
if (xdb.i[HI] < 0) return (x-x)/0.0; /* log(-x) = Nan */
/* Else subnormal number */
E = -64;
xdb.d *= two64; /* make x a normal number */
}
if (xdb.i[HI] >= 0x7ff00000) return x+x; /* Inf or Nan */
DOUBLE_EXTENDED_MODE; /* This one should be overlapped with following integer computation */
/* Extract exponent and mantissa */
E += (xdb.i[HI]>>20)-1023; /* extract the exponent */
index = (xdb.i[HI] & 0x000fffff);
xdb.i[HI] = index | 0x3ff00000; /* do exponent = 0 */
index = (index + (1<<(20-L-1))) >> (20-L);
/* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
if (index >= MAXINDEX){ /* corresponds to y>sqrt(2)*/
index = index & INDEXMASK;
xdb.i[HI] -= 0x00100000;
E++;
}
y = xdb.d;
#else /* defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64) */
/* Here come the code specific to Itanium processor */
PREFETCH_POLY_QUICK; /* defined in log-de.h */
y=x;
i = _Asm_getf(2/*_FR_D*/, y); /* Cast y to a 64-bit integer */
/* Filter special cases */
if (i<(int64_t)ULL(0010000000000000)){ /* equivalent to : x < 2^(-1022) */
if ((i & ULL(7fffffffffffffff))==0) return -1.0/0.0; /* log(+/-0) = -Inf */
if (i<0) return (x-x)/0.0; /* log(-x) = Nan */
/* Else subnormal number */
y *= two64; /* make x a normal number */
E = -64;
i = _Asm_getf(2/*_FR_D*/, y); /* and update i */
}
if (i >= ULL(7ff0000000000000)) return x+x; /* Inf or Nan */
/* Extract exponent and mantissa */
E += (i>>52)-1023;
i = i & ULL(000fffffffffffff); /* keep only mantissa */
index = (i + (ULL(1)<<(52-L-1))) >> (52-L);
/* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
if (index >= MAXINDEX){ /* corresponds to y>sqrt(2)*/
y = _Asm_setf(2/*_FR_D*/, (i | ULL(3ff0000000000000)) - ULL(0010000000000000) ); /* exponent = -1 */
index = index & INDEXMASK;
E++;
}
else
y = _Asm_setf(2/*_FR_D*/, i | ULL(3ff0000000000000) ); /* exponent = 0*/
#endif /* defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64) */
/* All the previous argument reduction was exact */
/* now y holds 1+f, and E is the exponent */
logirh = argredtable[index].logirh;
r = (double_ext) (argredtable[index].r); /* approx to 1/y.d */
z = y*r - 1. ; /* even without an FMA, all exact */
if(E==0)
roundtestmask=ACCURATE_TO_61_BITS;
else
roundtestmask=ACCURATE_TO_62_BITS;
#ifdef ESTRIN
/* Estrin polynomial evaluation */
double_ext z2,z4, p01, p23, p45, p67, p03, p47,p07;
z2 = z*z; p67 = c6 + z*c7; p45 = c4 + z*c5; p23 = c2 + z*c3; p01 = logirh + z;
z4 = z2*z2; p47 = p45 + z2*p67; p03 = p01 + z2*p23;
p07 = p03 + z4*p47;
logde = p07 + E*log2H;
#endif
#ifdef PATERSON
double_ext z4,z2,t0,t1,t2,t3,t4,t5,t6,t7,t8;
z2 = z * z; t1 = z + ps_alpha; t2 = z + ps_beta; t3 = c3 * z + c2; t4 = z + logirh;
z4 = z2 * z2; t5 = z2 + ps_c; t6 = t3 * z2 + t4;
t7 = t5 * t1 + t2; t0 = z4 * c7; t8 = t7 * t0 + t6;
logde = t8 + E*log2H;
#endif
#if 0 /* to time the first step only */
BACK_TO_DOUBLE_MODE; return (double)t;
#endif
/* To test the second step only, comment out the following line */
DE_TEST_AND_RETURN_RU(logde, roundtestmask);
log_accurate(&th, &tl, z, E, index);
RETURN_SUM_ROUNDED_UP(th, tl);
}
double log_rz(double x) {
if (x>1.0)
return log_rd(x);
else
return log_ru(x);
}