#28 Implementation of Math.IEEEremainder

[Wu Gansha]

The implementation of
Java_java_lang_Math_IEEEremainder is wrong in ORP, it
just calls fmod, which is corresponding to java
operator %, but doesn't do in IEEE 754 arithmetic rule.
So I reimplement this routine as the following:

inline jlong castDoubleToLong(jdouble d)
{
union {
jlong l;
jdouble d;
} U;
U.d = d;
return U.l;
} // castDoubleToLong

inline jdouble castLongToDouble(jlong l)
{
union {
jlong l;
jdouble d;
} U;
U.l = l;
return U.d;
} // castLongToDouble

inline bool is_nan(jlong long_bits)
{
jlong e = (long_bits & __INT64_C
(0x7ff0000000000000));
jlong g = (long_bits & __INT64_C
(0x000fffffffffffff));

if ((g != 0) && (e == __INT64_C
(0x7ff0000000000000)))
return true;
return false;
} // is_nan

inline bool is_positive_inf(jlong long_bits)
{
jlong e = (long_bits & __INT64_C
(0x7ff0000000000000));

if (e == __INT64_C(0x7ff0000000000000) )
return true;
return false;
} // is_positive_inf

inline bool is_negative_inf(jlong long_bits)
{
jlong e = (long_bits & __INT64_C
(0xfff0000000000000));

if (e == __INT64_C(0xfff0000000000000) )
return true;

return false;
} // is_negative_inf

inline bool is_inf(jlong long_bits)
{
return is_positive_inf(long_bits) ||
is_negative_inf(long_bits);
}

jdouble get_java_lang_Double_NAN(JNIEnv *jenv)
{
jclass double_class = jenv->FindClass
("java/lang/Double");
assert(double_class);
jfieldID NaN_id = jenv->GetStaticFieldID
(double_class, "NaN", "D");
assert(NaN_id);
jdouble NaN = jenv->GetStaticDoubleField
(double_class, NaN_id);
return NaN;
}

/******************************************************
******************
* The result of a Java floating-point remainder
operation is determined
* by the rules of IEEE arithmetic:
*
* 1.If either operand is NaN, the result is NaN.
* 2.If the result is not NaN, the sign of the
result equals the sign of
* the dividend.
* 3.If the dividend is an infinity, or the divisor
is a zero, or both,
* the result is NaN.
* 4.If the dividend is finite and the divisor is an
infinity, the result
* equals the dividend.
* 5.If the dividend is a zero and the divisor is
finite, the result equals
* the dividend.
* 6.In the remaining cases, where neither an
infinity, nor a zero, nor NaN
* is involved, the floating-point remainder r
from the division of a
* dividend n by a divisor d is defined by the
mathematical relation:
* r = n - (d * q)
* where q is the mathematical integer closest
to the exact mathematical
* value of the quotient n/d; if two mathematical
integers are equally
* close to then n is the integer that is even.
If the remainder is zero,
* its sign is the same as the sign of the first
argument.
*

*******************************************************
*******************/

/*
* Class: java_lang_Math
* Method: IEEEremainder
* Signature: (DD)D
*/
JNIEXPORT jdouble JNICALL
Java_java_lang_Math_IEEEremainder
(JNIEnv *jenv, jclass jclazz, jdouble x, jdouble y)
{
// It's wrong to use fmod to implement
IEEEremainder,
// fmod should be used to implement %
// double q = fmod(x, y);

jlong nl = castDoubleToLong(x);
jlong dl = castDoubleToLong(y);

// Step 1:
if(is_nan(nl) || is_nan(dl))
return get_java_lang_Double_NAN(jenv);

// Step 2:
jlong sign = 0x7fffffffffffffff | (nl &
0x8000000000000000);

// Step 3:
if( y == 0 || is_inf(nl) )
return get_java_lang_Double_NAN(jenv);

// Step 4:
if( !is_inf(nl) && is_inf(dl) )
return x;

// Step 5:
if( x == 0 && is_inf(dl) )
return 0;

// Step 6:
double q = x / y;
double q1 = floor(q);
double q2 = ceil(q);
double dq1 = q - q1;
double dq2 = q2 - q;
if(dq1 > dq2)
q = q2;
else
if(dq1 < dq2)
q = q1;
else
if((jlong)q1 & 0x1 == 0)
q = q1;
else
q = q2;
double r = x - q * y;
if(r == 0){
jlong rl = castDoubleToLong(r);
rl &= sign;
r = castLongToDouble(rl);
}

return r;
} // Java_java_lang_Math_IEEEremainder