in commons-math-core/src/main/java/org/apache/commons/math4/core/jdkmath/AccurateMath.java [2063:2272]
private static void reducePayneHanek(double x, double[] result) {
/* Convert input double to bits */
long inbits = Double.doubleToRawLongBits(x);
int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
/* Convert to fixed point representation */
inbits &= 0x000fffffffffffffL;
inbits |= 0x0010000000000000L;
/* Normalize input to be between 0.5 and 1.0 */
exponent++;
inbits <<= 11;
/* Based on the exponent, get a shifted copy of recip2pi */
long shpi0;
long shpiA;
long shpiB;
int idx = exponent >> 6;
int shift = exponent - (idx << 6);
if (shift != 0) {
shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx - 1] << shift);
shpi0 |= RECIP_2PI[idx] >>> (64 - shift);
shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx + 1] >>> (64 - shift));
shpiB = (RECIP_2PI[idx + 1] << shift) | (RECIP_2PI[idx + 2] >>> (64 - shift));
} else {
shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx - 1];
shpiA = RECIP_2PI[idx];
shpiB = RECIP_2PI[idx + 1];
}
/* Multiply input by shpiA */
long a = inbits >>> 32;
long b = inbits & 0xffffffffL;
long c = shpiA >>> 32;
long d = shpiA & 0xffffffffL;
long ac = a * c;
long bd = b * d;
long bc = b * c;
long ad = a * d;
long prodB = bd + (ad << 32);
long prodA = ac + (ad >>> 32);
boolean bita = (bd & 0x8000000000000000L) != 0;
boolean bitb = (ad & 0x80000000L) != 0;
boolean bitsum = (prodB & 0x8000000000000000L) != 0;
/* Carry */
if ((bita && bitb) ||
((bita || bitb) && !bitsum)) {
prodA++;
}
bita = (prodB & 0x8000000000000000L) != 0;
bitb = (bc & 0x80000000L) != 0;
prodB += bc << 32;
prodA += bc >>> 32;
bitsum = (prodB & 0x8000000000000000L) != 0;
/* Carry */
if ((bita && bitb) ||
((bita || bitb) && !bitsum)) {
prodA++;
}
/* Multiply input by shpiB */
c = shpiB >>> 32;
d = shpiB & 0xffffffffL;
ac = a * c;
bc = b * c;
ad = a * d;
/* Collect terms */
ac += (bc + ad) >>> 32;
bita = (prodB & 0x8000000000000000L) != 0;
bitb = (ac & 0x8000000000000000L) != 0;
prodB += ac;
bitsum = (prodB & 0x8000000000000000L) != 0;
/* Carry */
if ((bita && bitb) ||
((bita || bitb) && !bitsum)) {
prodA++;
}
/* Multiply by shpi0 */
c = shpi0 >>> 32;
d = shpi0 & 0xffffffffL;
bd = b * d;
bc = b * c;
ad = a * d;
prodA += bd + ((bc + ad) << 32);
/*
* prodA, prodB now contain the remainder as a fraction of PI. We want this as a fraction of
* PI/2, so use the following steps:
* 1.) multiply by 4.
* 2.) do a fixed point muliply by PI/4.
* 3.) Convert to floating point.
* 4.) Multiply by 2
*/
/* This identifies the quadrant */
int intPart = (int)(prodA >>> 62);
/* Multiply by 4 */
prodA <<= 2;
prodA |= prodB >>> 62;
prodB <<= 2;
/* Multiply by PI/4 */
a = prodA >>> 32;
b = prodA & 0xffffffffL;
c = PI_O_4_BITS[0] >>> 32;
d = PI_O_4_BITS[0] & 0xffffffffL;
ac = a * c;
bd = b * d;
bc = b * c;
ad = a * d;
long prod2B = bd + (ad << 32);
long prod2A = ac + (ad >>> 32);
bita = (bd & 0x8000000000000000L) != 0;
bitb = (ad & 0x80000000L) != 0;
bitsum = (prod2B & 0x8000000000000000L) != 0;
/* Carry */
if ((bita && bitb) ||
((bita || bitb) && !bitsum)) {
prod2A++;
}
bita = (prod2B & 0x8000000000000000L) != 0;
bitb = (bc & 0x80000000L) != 0;
prod2B += bc << 32;
prod2A += bc >>> 32;
bitsum = (prod2B & 0x8000000000000000L) != 0;
/* Carry */
if ((bita && bitb) ||
((bita || bitb) && !bitsum)) {
prod2A++;
}
/* Multiply input by pio4bits[1] */
c = PI_O_4_BITS[1] >>> 32;
d = PI_O_4_BITS[1] & 0xffffffffL;
ac = a * c;
bc = b * c;
ad = a * d;
/* Collect terms */
ac += (bc + ad) >>> 32;
bita = (prod2B & 0x8000000000000000L) != 0;
bitb = (ac & 0x8000000000000000L) != 0;
prod2B += ac;
bitsum = (prod2B & 0x8000000000000000L) != 0;
/* Carry */
if ((bita && bitb) ||
((bita || bitb) && !bitsum)) {
prod2A++;
}
/* Multiply inputB by pio4bits[0] */
a = prodB >>> 32;
b = prodB & 0xffffffffL;
c = PI_O_4_BITS[0] >>> 32;
d = PI_O_4_BITS[0] & 0xffffffffL;
ac = a * c;
bc = b * c;
ad = a * d;
/* Collect terms */
ac += (bc + ad) >>> 32;
bita = (prod2B & 0x8000000000000000L) != 0;
bitb = (ac & 0x8000000000000000L) != 0;
prod2B += ac;
bitsum = (prod2B & 0x8000000000000000L) != 0;
/* Carry */
if ((bita && bitb) ||
((bita || bitb) && !bitsum)) {
prod2A++;
}
/* Convert to double */
double tmpA = (prod2A >>> 12) / TWO_POWER_52; // High order 52 bits
double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low bits
double sumA = tmpA + tmpB;
double sumB = -(sumA - tmpA - tmpB);
/* Multiply by PI/2 and return */
result[0] = intPart;
result[1] = sumA * 2.0;
result[2] = sumB * 2.0;
}