RFR: JDK-8304028: Port fdlibm IEEEremainder to Java
Joe Darcy
darcy at openjdk.org
Tue Mar 21 06:19:48 UTC 2023
On Tue, 21 Mar 2023 06:11:57 GMT, Joe Darcy <darcy at openjdk.org> wrote:
> Last but not least, a port of fdlibm IEEEremainder from C to Java. I plan to write some more implementation-specific tests around decision points in the FDLIBM algorithm, but I wanted to get the bulk of the changes out for review first.
>
> Note that since IEEEremainder was the last native method in StrictMath.java, the StrictMath.c file needed to be deleted (or modified) since StrictMath.h was no longer generated as part of the build. (StrictMath.c was one of the file deleted as part of JDK-8302801).
>
> For testing, Mach 5 tier 1 through 3 were successful (other than an unrelated test failure that was problem listed) and the exhaustive test was locally run and passed with "16, 16" to increase the testing density.
Diffs to aid review:
$ diff -w Remainder.c Remainder.translit.java
1,25c1,6
< /* __ieee754_remainder(x,p)
< * Return :
< * returns x REM p = x - [x/p]*p as if in infinite
< * precise arithmetic, where [x/p] is the (infinite bit)
< * integer nearest x/p (in half way case choose the even one).
< * Method :
< * Based on fmod() return x-[x/p]chopped*p exactlp.
< */
<
< #include "fdlibm.h"
<
< #ifdef __STDC__
< static const double zero = 0.0;
< #else
< static double zero = 0.0;
< #endif
<
<
< #ifdef __STDC__
< double __ieee754_remainder(double x, double p)
< #else
< double __ieee754_remainder(x,p)
< double x,p;
< #endif
< {
---
> private static final class IEEEremainder {
> private static final double zero = 0.0;
> private static double one = 1.0;
> private static double[] Zero = {0.0, -0.0,};
>
> static double compute(double x, double p) {
27c8
< unsigned sx,lx,lp;
---
> /*unsigned*/ int sx,lx,lp;
48,49c29,30
< x = fabs(x);
< p = fabs(p);
---
> x = Math.abs(x);
> p = Math.abs(p);
62c43,44
< __HI(x) ^= sx;
---
> // __HI(x) ^= sx;
> x = __HI(x, __HI(x)^sx);
65,85c47,48
< /*
< * __ieee754_fmod(x,y)
< * Return x mod y in exact arithmetic
< * Method: shift and subtract
< */
<
< #include "fdlibm.h"
<
< #ifdef __STDC__
< static const double one = 1.0, Zero[] = {0.0, -0.0,};
< #else
< static double one = 1.0, Zero[] = {0.0, -0.0,};
< #endif
<
< #ifdef __STDC__
< double __ieee754_fmod(double x, double y)
< #else
< double __ieee754_fmod(x,y)
< double x,y ;
< #endif
< {
---
>
> private static double __ieee754_fmod(double x, double y) {
87c50
< unsigned lx,ly,lz;
---
> /*unsigned*/ int lx,ly,lz;
102c65,66
< if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
---
> // if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
> if((hx<hy)||(Integer.compareUnsigned(lx,ly) < 0)) return x; /* |x|<|y| return x */
104c68
< return Zero[(unsigned)sx>>31]; /* |x|=|y| return x*0*/
---
> return Zero[/*(unsigned)*/sx>>>31]; /* |x|=|y| return x*0*/ // unsigned shift
131c95
< hx = (hx<<n)|(lx>>(32-n));
---
> hx = (hx<<n)|(lx >>> (32-n)); // unsigned shift
143c107
< hy = (hy<<n)|(ly>>(32-n));
---
> hy = (hy<<n)|(ly >>> (32-n)); // unsigned shift
153,155c117,121
< while(n--) {
< hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
< if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
---
> while(n-- != 0) {
> hz=hx-hy;lz=lx-ly;
> // if(lx<ly) hz -= 1;
> if(Integer.compareUnsigned(lx, ly) < 0) hz -= 1;
> if(hz<0){hx = hx+hx+(lx >>> 31); lx = lx+lx;} // unsigned shift
158,159c124,126
< return Zero[(unsigned)sx>>31];
< hx = hz+hz+(lz>>31); lx = lz+lz;
---
> return Zero[/*(unsigned)*/sx>>>31]; // unsigned shift
> hx = hz+hz+(lz >>> 31); // unsigned shift
> lx = lz+lz;
162c129,131
< hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
---
> hz=hx-hy;lz=lx-ly;
> // if(lx<ly) hz -= 1;
> if(Integer.compareUnsigned(lx, ly) < 0) hz -= 1;
167c136
< return Zero[(unsigned)sx>>31];
---
> return Zero[/*(unsigned)*/sx >>> 31]; // unsigned shift
169c138
< hx = hx+hx+(lx>>31); lx = lx+lx;
---
> hx = hx+hx+(lx >>> 31); lx = lx+lx; // unsigned shift
174,175c143,146
< __HI(x) = hx|sx;
< __LO(x) = lx;
---
> // __HI(x) = hx|sx;
> x = __HI(x, hx|sx);
> // __LO(x) = lx;
> x = __LO(x, lx);
179c150
< lx = (lx>>n)|((unsigned)hx<<(32-n));
---
> lx = (lx >>> n)|(/*(unsigned)*/hx<<(32-n)); // unsigned shift
182c153,154
< lx = (hx<<(32-n))|(lx>>n); hx = sx;
---
> lx = (hx<<(32-n))|(lx >>> n); // unsigned shift
> hx = sx;
186,187c158,161
< __HI(x) = hx|sx;
< __LO(x) = lx;
---
> // __HI(x) = hx|sx;
> x = __HI(x, hx|sx);
> // __LO(x) = lx;
> x = __LO(x, lx);
191a166
> }
$ diff -w Remainder.translit.java Remainder.fdlibm.java
1,4c1,2
< private static final class IEEEremainder {
< private static final double zero = 0.0;
< private static double one = 1.0;
< private static double[] Zero = {0.0, -0.0,};
---
> static final class IEEEremainder {
> private IEEEremainder() {throw new UnsupportedOperationException();}
11,24c9,15
< hx = __HI(x); /* high word of x */
< lx = __LO(x); /* low word of x */
< hp = __HI(p); /* high word of p */
< lp = __LO(p); /* low word of p */
< sx = hx&0x80000000;
< hp &= 0x7fffffff;
< hx &= 0x7fffffff;
<
< /* purge off exception values */
< if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */
< if((hx>=0x7ff00000)|| /* x not finite */
< ((hp>=0x7ff00000)&& /* p is NaN */
< (((hp-0x7ff00000)|lp)!=0)))
< return (x*p)/(x*p);
---
> hx = __HI(x); // high word of x
> lx = __LO(x); // low word of x
> hp = __HI(p); // high word of p
> lp = __LO(p); // low word of p
> sx = hx & 0x8000_0000;
> hp &= 0x7fff_ffff;
> hx &= 0x7fff_ffff;
25a17,24
> // purge off exception values
> if ((hp | lp) == 0) {// p = 0
> return (x*p)/(x*p);
> }
> if ((hx >= 0x7ff0_0000) || // not finite
> ((hp >= 0x7ff0_0000) && // p is NaN
> (((hp - 0x7ff0_0000) | lp) != 0)))
> return (x*p)/(x*p);
27,28c26,31
< if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */
< if (((hx-hp)|(lx-lp))==0) return zero*x;
---
> if (hp <= 0x7fdf_ffff) { // now x < 2p
> x = __ieee754_fmod(x, p + p);
> }
> if (((hx - hp) | (lx - lp)) == 0) {
> return 0.0*x;
> }
31c34
< if (hp<0x00200000) {
---
> if (hp < 0x0020_0000) {
34c37,39
< if(x+x>=p) x -= p;
---
> if (x + x >= p) {
> x -= p;
> }
40c45,47
< if(x>=p_half) x -= p;
---
> if (x >= p_half) {
> x -= p;
> }
43d49
< // __HI(x) ^= sx;
49c55
< int n,hx,hy,hz,ix,iy,sx,i;
---
> int n, hx, hy, hz, ix, iy, sx;
52,62c58,68
< hx = __HI(x); /* high word of x */
< lx = __LO(x); /* low word of x */
< hy = __HI(y); /* high word of y */
< ly = __LO(y); /* low word of y */
< sx = hx&0x80000000; /* sign of x */
< hx ^=sx; /* |x| */
< hy &= 0x7fffffff; /* |y| */
<
< /* purge off exception values */
< if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
< ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */
---
> hx = __HI(x); // high word of x
> lx = __LO(x); // low word of x
> hy = __HI(y); // high word of y
> ly = __LO(y); // low word of y
> sx = hx & 0x8000_0000; // sign of x
> hx ^= sx; // |x|
> hy &= 0x7fff_ffff; // |y|
>
> // purge off exception values
> if((hy | ly) == 0 || (hx >= 0x7ff0_0000)|| // y = 0, or x not finite
> ((hy | ((ly | -ly) >> 31)) > 0x7ff0_0000)) // or y is NaN
65,68c71,72
< // if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
< if((hx<hy)||(Integer.compareUnsigned(lx,ly) < 0)) return x; /* |x|<|y| return x */
< if(lx==ly)
< return Zero[/*(unsigned)*/sx>>>31]; /* |x|=|y| return x*0*/ // unsigned shift
---
> if ((hx < hy) || (Integer.compareUnsigned(lx, ly) < 0)) { // |x| < |y| return x
> return x;
70,76c74,75
<
< /* determine ix = ilogb(x) */
< if(hx<0x00100000) { /* subnormal x */
< if(hx==0) {
< for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
< } else {
< for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
---
> if (lx == ly) {
> return signedZero(sx); // |x| = |y| return x*0
78,85d76
< } else ix = (hx>>20)-1023;
<
< /* determine iy = ilogb(y) */
< if(hy<0x00100000) { /* subnormal y */
< if(hy==0) {
< for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
< } else {
< for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
87d77
< } else iy = (hy>>20)-1023;
89c79,82
< /* set up {hx,lx}, {hy,ly} and align y to x */
---
> ix = ilogb(hx, lx);
> iy = ilogb(hy, ly);
>
> // set up {hx, lx}, {hy, ly} and align y to x
91,92c84,85
< hx = 0x00100000|(0x000fffff&hx);
< else { /* subnormal x, shift x to normal */
---
> hx = 0x0010_0000 | (0x000f_ffff & hx);
> else { // subnormal x, shift x to normal
103,104c96,97
< hy = 0x00100000|(0x000fffff&hy);
< else { /* subnormal y, shift y to normal */
---
> hy = 0x0010_0000 | (0x000f_ffff & hy);
> else { // subnormal y, shift y to normal
115c108
< /* fix point fmod */
---
> // fix point fmod
118,124c111,122
< hz=hx-hy;lz=lx-ly;
< // if(lx<ly) hz -= 1;
< if(Integer.compareUnsigned(lx, ly) < 0) hz -= 1;
< if(hz<0){hx = hx+hx+(lx >>> 31); lx = lx+lx;} // unsigned shift
< else {
< if((hz|lz)==0) /* return sign(x)*0 */
< return Zero[/*(unsigned)*/sx>>>31]; // unsigned shift
---
> hz = hx - hy;
> lz = lx - ly;
> if (Integer.compareUnsigned(lx, ly) < 0) {
> hz -= 1;
> }
> if (hz < 0){
> hx = hx + hx +(lx >>> 31); // unsigned shift
> lx = lx + lx;
> } else {
> if ((hz | lz) == 0) { // return sign(x)*0
> return signedZero(sx);
> }
129,138c127,143
< hz=hx-hy;lz=lx-ly;
< // if(lx<ly) hz -= 1;
< if(Integer.compareUnsigned(lx, ly) < 0) hz -= 1;
< if(hz>=0) {hx=hz;lx=lz;}
<
< /* convert back to floating value and restore the sign */
< if((hx|lx)==0) /* return sign(x)*0 */
< return Zero[/*(unsigned)*/sx >>> 31]; // unsigned shift
< while(hx<0x00100000) { /* normalize x */
< hx = hx+hx+(lx >>> 31); lx = lx+lx; // unsigned shift
---
> hz = hx - hy;
> lz = lx - ly;
> if (Integer.compareUnsigned(lx, ly) < 0) {
> hz -= 1;
> }
> if (hz >= 0) {
> hx = hz;
> lx = lz;
> }
>
> // convert back to floating value and restore the sign
> if ((hx | lx) == 0) { // return sign(x)*0
> return signedZero(sx);
> }
> while (hx < 0x0010_0000) { // normalize x
> hx = hx + hx + (lx >>> 31); // unsigned shift
> lx = lx + lx;
141,147c146,149
< if(iy>= -1022) { /* normalize output */
< hx = ((hx-0x00100000)|((iy+1023)<<20));
< // __HI(x) = hx|sx;
< x = __HI(x, hx|sx);
< // __LO(x) = lx;
< x = __LO(x, lx);
< } else { /* subnormal output */
---
> if( iy >= -1022) { // normalize output
> hx = ((hx - 0x0010_0000) | ((iy + 1023) << 20));
> x = __HI_LO(hx | sx, lx);
> } else { // subnormal output
156c158,184
< lx = hx>>(n-32); hx = sx;
---
> lx = hx >>(n - 32);
> hx = sx;
> }
> x = __HI_LO(hx | sx, lx);
> x *= 1.0; // create necessary signal
> }
> return x; // exact output
> }
>
> /*
> * Return a double zero with the same sign as the int argument.
> */
> private static double signedZero(int sign) {
> return +0.0 * ( (double)sign);
> }
>
> private static int ilogb(int hz, int lz) {
> int iz, i;
> if (hz < 0x0010_0000) { // subnormal z
> if (hz == 0) {
> for (iz = -1043, i = lz; i > 0; i <<= 1) {
> iz -= 1;
> }
> } else {
> for (iz = -1022, i = (hz << 11); i > 0; i <<= 1) {
> iz -= 1;
> }
158,162c186,187
< // __HI(x) = hx|sx;
< x = __HI(x, hx|sx);
< // __LO(x) = lx;
< x = __LO(x, lx);
< x *= one; /* create necessary signal */
---
> } else {
> iz = (hz >> 20) - 1023;
164c189
< return x; /* exact output */
---
> return iz;
-------------
PR Comment: https://git.openjdk.org/jdk/pull/13113#issuecomment-1477333025
More information about the core-libs-dev
mailing list