RFR: JDK-8301396: Port fdlibm expm1 to Java
Raffaello Giulietti
rgiulietti at openjdk.org
Fri Feb 3 12:52:49 UTC 2023
On Thu, 2 Feb 2023 22:18:32 GMT, Joe Darcy <darcy at openjdk.org> wrote:
> Next on the FDLIBM C -> Java port, expm1.
> Coming soon, hyperbolic transcendentals (sinh, cosh, tanh)!
>
> For expm1, the C vs transliteration port show the usual kind of differences, beside formatting of the constants, the use of the __HI macro on the left-hand side needs to be replaced by a method call and an assignment, as seen below:
>
>
> <
> < #include "fdlibm.h"
> <
> < #ifdef __STDC__
> < static const double
> < #else
> < static double
> < #endif
> < one = 1.0,
> < huge = 1.0e+300,
> < tiny = 1.0e-300,
> < o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */
> < ln2_hi = 6.93147180369123816490e-01,/* 0x3fe62e42, 0xfee00000 */
> < ln2_lo = 1.90821492927058770002e-10,/* 0x3dea39ef, 0x35793c76 */
> < invln2 = 1.44269504088896338700e+00,/* 0x3ff71547, 0x652b82fe */
> ---
>> static class Expm1 {
>> private static final double one = 1.0;
>> private static final double huge = 1.0e+300;
>> private static final double tiny = 1.0e-300;
>> private static final double o_threshold = 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */
>> private static final double ln2_hi = 6.93147180369123816490e-01; /* 0x3fe62e42, 0xfee00000 */
>> private static final double ln2_lo = 1.90821492927058770002e-10; /* 0x3dea39ef, 0x35793c76 */
>> private static final double invln2 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */
> 111,115c104,108
> < Q1 = -3.33333333333331316428e-02, /* BFA11111 111110F4 */
> < Q2 = 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */
> < Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
> < Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
> < Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
> ---
>> private static final double Q1 = -3.33333333333331316428e-02; /* BFA11111 111110F4 */
>> private static final double Q2 = 1.58730158725481460165e-03; /* 3F5A01A0 19FE5585 */
>> private static final double Q3 = -7.93650757867487942473e-05; /* BF14CE19 9EAADBB7 */
>> private static final double Q4 = 4.00821782732936239552e-06; /* 3ED0CFCA 86E65239 */
>> private static final double Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
> 117,123c110
> < #ifdef __STDC__
> < double expm1(double x)
> < #else
> < double expm1(x)
> < double x;
> < #endif
> < {
> ---
>> static double compute(double x) {
> 126c113
> < unsigned hx;
> ---
>> /*unsigned*/ int hx;
> 157c144
> < k = invln2*x+((xsb==0)?0.5:-0.5);
> ---
>> k = (int)(invln2*x+((xsb==0)?0.5:-0.5));
> 188c175
> < __HI(y) += (k<<20); /* add k to y's exponent */
> ---
>> y = __HI(y, __HI(y) + (k<<20)); /* add k to y's exponent */
> 193c180
> < __HI(t) = 0x3ff00000 - (0x200000>>k); /* t=1-2^-k */
> ---
>> t = __HI(t, 0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */
> 195c182
> < __HI(y) += (k<<20); /* add k to y's exponent */
> ---
>> y = __HI(y, __HI(y) + (k<<20)); /* add k to y's exponent */
> 197c184
> < __HI(t) = ((0x3ff-k)<<20); /* 2^-k */
> ---
>> t = __HI(t, ((0x3ff-k)<<20)); /* 2^-k */
> 200c187
> < __HI(y) += (k<<20); /* add k to y's exponent */
> ---
>> y = __HI(y, __HI(y) + (k<<20)); /* add k to y's exponent */
> 205c192
> <
> ---
>> }
>
>
> When comparing the transliteration port and the more idiomatic port, there were no surprising or notable differences:
>
>
> $ diff -w Expm1.translit.java Expm1.fdlibm.java
> 99,102c99,102
> < private static final double o_threshold = 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */
> < private static final double ln2_hi = 6.93147180369123816490e-01; /* 0x3fe62e42, 0xfee00000 */
> < private static final double ln2_lo = 1.90821492927058770002e-10; /* 0x3dea39ef, 0x35793c76 */
> < private static final double invln2 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */
> ---
>> private static final double o_threshold = 0x1.62e42fefa39efp9; // 7.09782712893383973096e+02
>> private static final double ln2_hi = 0x1.62e42feep-1; // 6.93147180369123816490e-01
>> private static final double ln2_lo = 0x1.a39ef35793c76p-33; // 1.90821492927058770002e-10
>> private static final double invln2 = 0x1.71547652b82fep0; // 1.44269504088896338700e+00
> 104,108c104,108
> < private static final double Q1 = -3.33333333333331316428e-02; /* BFA11111 111110F4 */
> < private static final double Q2 = 1.58730158725481460165e-03; /* 3F5A01A0 19FE5585 */
> < private static final double Q3 = -7.93650757867487942473e-05; /* BF14CE19 9EAADBB7 */
> < private static final double Q4 = 4.00821782732936239552e-06; /* 3ED0CFCA 86E65239 */
> < private static final double Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
> ---
>> private static final double Q1 = -0x1.11111111110f4p-5; // -3.33333333333331316428e-02
>> private static final double Q2 = 0x1.a01a019fe5585p-10; // 1.58730158725481460165e-03
>> private static final double Q3 = -0x1.4ce199eaadbb7p-14; // -7.93650757867487942473e-05
>> private static final double Q4 = 0x1.0cfca86e65239p-18; // 4.00821782732936239552e-06
>> private static final double Q5 = -0x1.afdb76e09c32dp-23; // -2.01099218183624371326e-07
> 116,118c116,118
> < xsb = hx&0x80000000; /* sign bit of x */
> < if(xsb==0) y=x; else y= -x; /* y = |x| */
> < hx &= 0x7fffffff; /* high word of |x| */
> ---
>> xsb = hx & 0x8000_0000; /* sign bit of x */
>> y = Math.abs(x);
>> hx &= 0x7fff_ffff; /* high word of |x| */
> 121,124c121,124
> < if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */
> < if(hx >= 0x40862E42) { /* if |x|>=709.78... */
> < if(hx>=0x7ff00000) {
> < if(((hx&0xfffff)|__LO(x))!=0)
> ---
>> if (hx >= 0x4043_687A) { /* if |x| >= 56*ln2 */
>> if (hx >= 0x4086_2E42) { /* if |x| >= 709.78... */
>> if (hx >= 0x7ff_00000) {
>> if (((hx & 0xf_ffff) | __LO(x)) != 0) {
> 126c126,131
> < else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
> ---
>> } else {
>> return (xsb == 0)? x : -1.0; /* exp(+-inf)={inf,-1} */
>> }
>> }
>> if (x > o_threshold) {
>> return huge*huge; /* overflow */
> 128d132
> < if(x > o_threshold) return huge*huge; /* overflow */
> 131c135
> < if(x+tiny<0.0) /* raise inexact */
> ---
>> if (x + tiny < 0.0) { /* raise inexact */
> 134a139
>> }
> 137,142c142,152
> < if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
> < if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
> < if(xsb==0)
> < {hi = x - ln2_hi; lo = ln2_lo; k = 1;}
> < else
> < {hi = x + ln2_hi; lo = -ln2_lo; k = -1;}
> ---
>> if (hx > 0x3fd6_2e42) { /* if |x| > 0.5 ln2 */
>> if (hx < 0x3FF0_A2B2) { /* and |x| < 1.5 ln2 */
>> if (xsb == 0) {
>> hi = x - ln2_hi;
>> lo = ln2_lo;
>> k = 1;}
>> else {
>> hi = x + ln2_hi;
>> lo = -ln2_lo;
>> k = -1;
>> }
> 151,152c161
> < }
> < else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */
> ---
>> } else if (hx < 0x3c90_0000) { /* when |x|<2**-54, return x */
> 154a164,165
>> } else {
>> k = 0;
> 156d166
> < else k = 0;
> 164,165c174,176
> < if(k==0) return x - (x*e-hxs); /* c is 0 */
> < else {
> ---
>> if (k == 0) {
>> return x - (x*e - hxs); /* c is 0 */
>> } else {
> 168c179,181
> < if(k== -1) return 0.5*(x-e)-0.5;
> ---
>> if (k == -1) {
>> return 0.5*(x - e) - 0.5;
>> }
> 170,171c183,187
> < if(x < -0.25) return -2.0*(e-(x+0.5));
> < else return one+2.0*(x-e);
> ---
>> if(x < -0.25) {
>> return -2.0*(e - (x + 0.5));
>> } else {
>> return one + 2.0*(x - e);
>> }
> 180c196
> < t = __HI(t, 0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */
> ---
>> t = __HI(t, 0x3ff0_0000 - (0x2_00000 >> k)); /* t=1-2^-k */
With the corrections by @turbanoff LGTM.
-------------
PR: https://git.openjdk.org/jdk/pull/12394
More information about the core-libs-dev
mailing list