On Tue, 21 Mar 2023 06:11:57 GMT, Joe Darcy <da...@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