private static void reducePayneHanek(double x,double result[]){
  long inbits=Double.doubleToRawLongBits(x);
  int exponent=(int)((inbits >> 52) & 0x7ff) - 1023;
  inbits&=0x000fffffffffffffL;
  inbits|=0x0010000000000000L;
  exponent++;
  inbits<<=11;
  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];
  }
  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;
  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;
  if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
    prodA++;
  }
  c=shpiB >>> 32;
  d=shpiB & 0xffffffffL;
  ac=a * c;
  bc=b * c;
  ad=a * d;
  ac+=(bc + ad) >>> 32;
  bita=(prodB & 0x8000000000000000L) != 0;
  bitb=(ac & 0x8000000000000000L) != 0;
  prodB+=ac;
  bitsum=(prodB & 0x8000000000000000L) != 0;
  if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
    prodA++;
  }
  c=shpi0 >>> 32;
  d=shpi0 & 0xffffffffL;
  bd=b * d;
  bc=b * c;
  ad=a * d;
  prodA+=bd + ((bc + ad) << 32);
  int intPart=(int)(prodA >>> 62);
  prodA<<=2;
  prodA|=prodB >>> 62;
  prodB<<=2;
  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;
  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;
  if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
    prod2A++;
  }
  c=PI_O_4_BITS[1] >>> 32;
  d=PI_O_4_BITS[1] & 0xffffffffL;
  ac=a * c;
  bc=b * c;
  ad=a * d;
  ac+=(bc + ad) >>> 32;
  bita=(prod2B & 0x8000000000000000L) != 0;
  bitb=(ac & 0x8000000000000000L) != 0;
  prod2B+=ac;
  bitsum=(prod2B & 0x8000000000000000L) != 0;
  if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
    prod2A++;
  }
  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;
  ac+=(bc + ad) >>> 32;
  bita=(prod2B & 0x8000000000000000L) != 0;
  bitb=(ac & 0x8000000000000000L) != 0;
  prod2B+=ac;
  bitsum=(prod2B & 0x8000000000000000L) != 0;
  if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
    prod2A++;
  }
  double tmpA=(prod2A >>> 12) / TWO_POWER_52;
  double tmpB=(((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52;
  double sumA=tmpA + tmpB;
  double sumB=-(sumA - tmpA - tmpB);
  result[0]=intPart;
  result[1]=sumA * 2.0;
  result[2]=sumB * 2.0;
}

private static void reducePayneHanek ( double x , double result [ ] ) { long inbits = Double . doubleToRawLongBits ( x ) ; int exponent = ( int ) ( ( inbits > > 52 ) & 0x7ff ) - 1023 ; inbits &= 0x000fffffffffffffL ; inbits |= 0x0010000000000000L ; exponent ++ ; inbits <<= 11 ; 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 ] ; } 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 ; 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 ; if ( ( bita && bitb ) || ( ( bita || bitb ) && ! bitsum ) ) { prodA ++ ; } c = shpiB > > > 32 ; d = shpiB & 0xffffffffL ; ac = a * c ; bc = b * c ; ad = a * d ; ac += ( bc + ad ) > > > 32 ; bita = ( prodB & 0x8000000000000000L ) != 0 ; bitb = ( ac & 0x8000000000000000L ) != 0 ; prodB += ac ; bitsum = ( prodB & 0x8000000000000000L ) != 0 ; if ( ( bita && bitb ) || ( ( bita || bitb ) && ! bitsum ) ) { prodA ++ ; } c = shpi0 > > > 32 ; d = shpi0 & 0xffffffffL ; bd = b * d ; bc = b * c ; ad = a * d ; prodA += bd + ( ( bc + ad ) << 32 ) ; int intPart = ( int ) ( prodA > > > 62 ) ; prodA <<= 2 ; prodA |= prodB > > > 62 ; prodB <<= 2 ; 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 ; 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 ; if ( ( bita && bitb ) || ( ( bita || bitb ) && ! bitsum ) ) { prod2A ++ ; } c = PI_O_4_BITS [ 1 ] > > > 32 ; d = PI_O_4_BITS [ 1 ] & 0xffffffffL ; ac = a * c ; bc = b * c ; ad = a * d ; ac += ( bc + ad ) > > > 32 ; bita = ( prod2B & 0x8000000000000000L ) != 0 ; bitb = ( ac & 0x8000000000000000L ) != 0 ; prod2B += ac ; bitsum = ( prod2B & 0x8000000000000000L ) != 0 ; if ( ( bita && bitb ) || ( ( bita || bitb ) && ! bitsum ) ) { prod2A ++ ; } 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 ; ac += ( bc + ad ) > > > 32 ; bita = ( prod2B & 0x8000000000000000L ) != 0 ; bitb = ( ac & 0x8000000000000000L ) != 0 ; prod2B += ac ; bitsum = ( prod2B & 0x8000000000000000L ) != 0 ; if ( ( bita && bitb ) || ( ( bita || bitb ) && ! bitsum ) ) { prod2A ++ ; } double tmpA = ( prod2A > > > 12 ) / TWO_POWER_52 ; double tmpB = ( ( ( prod2A & 0xfffL ) << 40 ) + ( prod2B > > > 24 ) ) / TWO_POWER_52 / TWO_POWER_52 ; double sumA = tmpA + tmpB ; double sumB = - ( sumA - tmpA - tmpB ) ; result [ 0 ] = intPart ; result [ 1 ] = sumA * 2.0 ; result [ 2 ] = sumB * 2.0 ; }
reduce the input argument using the payne and hanek method . this is good for all inputs 0 . 0 < x < inf output is remainder after dividing by pi / 2 the result array should contain 3 numbers . result [ 0 ] is the integer portion , so mod 4 this gives the quadrant . result [ 1 ] is the upper bits of the remainder result [ 2 ] is the lower bits of the remainder
