On Tue, 21 Mar 2023 06:11:57 GMT, Joe Darcy <[email protected]> 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
