
import math.native.bigfloat.*;

//
// This file include functions to compute y=x-N*pi/2 and return the last two bits of N
// in order to know which quadrant we are considering.
//
// Computed by Payne and Hanek methods. For more information
// you can read K. C. Ng research report from Sun Microsystems:
// "Argument reduction for huge argument: Good to the last bit" (July 13, 1992)
//

private alias two_over_pi = array(
   0x28be60dbu, 0x24e44152u, 0x27f09d5fu, 0x11f534ddu,
   0x3036d8a5u, 0x1993c439u, 0x0107f945u, 0x23abdebbu,
   0x31586dc9u, 0x06e3a424u, 0x374b8019u, 0x092eea09u,
   0x3464873fu, 0x21deb1cbu, 0x04a69cfbu, 0x288235f5u,
   0x0baed121u, 0x0e99c702u, 0x1ad17df9u, 0x013991d6u,
   0x0e60d4ceu, 0x1f49c845u, 0x3e2ef7e4u, 0x283b1ff8u,
   0x25fff781u, 0x1980fef2u, 0x3c462d68u, 0x0a6d1f6du,
   0x0d9fb3c9u, 0x3cb09b74u, 0x3d18fd9au, 0x1e5fea2du,
   0x1d49eeb1u, 0x3ebe5f17u, 0x2cf41ce7u, 0x378a5292u,
   0x3a9afed7u, 0x3b11f8d5u, 0x3421580cu, 0x3046fc7bu,
   0x1aeafc33u, 0x3bc209afu, 0x10d876a7u, 0x2391615eu,
   0x3986c219u, 0x199855f1u, 0x1281a102u, 0x0dffd880u);



// This number stores 211 bits of pi/2

private alias Pio2 = BigFloat(array(0x00000001u, 0x2487ed51u, 0x042d1846u, 0x26263314u,
    0x1701b839u, 0x28948127u, 0x01114cf9u, 0x23a0105du),
    1.,  0,   1 );

//
// Case X_IND = -1:
//               0
//           2 ^
//  X    :   <> |--| |--| |--|   0   0 0 0 0
//  2/Pi :   <> |--| |--| |--| |--| .....
//
//  Case X_IND = 0:
//                   0
//                2 ^
//  X    :   |--| <> |--| |--|   0   0 0 0 0
//  2/Pi :        <> |--| |--| |--| |--| .....
//
//  Case X_IND = 1:
//                       0
//                    2 ^
//  X    :   |--| |--| <> |--| |--|   0   0 0 0 0
//  2/Pi :             <> |--| |--| |--| |--| .....
//
//  Case ...
//
//  Step 1:
//
//   Compute r = X . 2/Pi where:
//    - r[0] hold the integer part. (if x>0 or the once complement integer part if x<0 )
//    - r[1] to r[SCS_NB_WORDS+2] hold the reduced part
//      the 3 extra 30 bits are here to prevent possible
//      cancellation due to a number x too close to a
//      multiple of Pi/2.
//
//  Step 2:
//   Compute result = (r[1] ... r[SCS_NB_WORDS]) . Pi/2.
//
// description of local variables :
// - ind : where to start multiplying into 2opi table
//
//

// TODO OPTIM
//   better 64-bit multiplication, see in scs_mult

rem_pio2_bf(result, x){
    var r,tmp,N,sign = Array[UInt64,BF_NB_WORDS+3](), 0ul, 0u,0;
    // result r[0],...,r[10] could store till 300 bits of precision
    // that is really enough for computing the reduced argument
    if ((x.excp != 1)or(x.index < -1)){
        result = x;
        return 0;
    }

    // Compute the product |x| * 2/Pi
    if ((x.index == -1)){
        // In this case we consider number between ]-1,+1[
        // we may use simpler algorithm such as Cody And Waite
        r[0] =  0;
        r[1] =  0;
        r[2] = (two_over_pi[0]) * x.h_word[0];
        r[3] = ((two_over_pi[0]) * x.h_word[1]+(two_over_pi[1]) * x.h_word[0]);
        if(x.h_word[2] == 0){
            for(i in range(4,(BF_NB_WORDS+3))){
                r[i] = ((two_over_pi[i-3]) * x.h_word[1]
                +(two_over_pi[i-2]) * x.h_word[0]);
            }
        }else {
            for(i in range(4,BF_NB_WORDS+3)){
                r[i] = ((two_over_pi[i-4]) * x.h_word[2]
                +(two_over_pi[i-3]) * x.h_word[1]
                +(two_over_pi[i-2]) * x.h_word[0]);
            }
        }
    }else {
        if (x.index == 0){
            r[0] =  0;
            r[1] =  (two_over_pi[0]) * x.h_word[0];
            r[2] = ((two_over_pi[0]) * x.h_word[1]
            +(two_over_pi[1]) * x.h_word[0]);
            if(x.h_word[2] == 0){
                for(i in range(3,BF_NB_WORDS+3)){
                    r[i] = ((two_over_pi[i-2]) * x.h_word[1]
                    +(two_over_pi[i-1]) * x.h_word[0]);
                }
            }else {
                for(i in range(3,BF_NB_WORDS+3)){
                    r[i] = ((two_over_pi[i-3]) * x.h_word[2]
                       +(two_over_pi[i-2]) * x.h_word[1]
                       +(two_over_pi[i-1]) * x.h_word[0]);
                }
            }
        }else {
            if (x.index == 1){
                r[0] =  (two_over_pi[0]) * x.h_word[0];
                r[1] = ((two_over_pi[0]) * x.h_word[1]
                       +(two_over_pi[1]) * x.h_word[0]);
                if(x.h_word[2] == 0){
                    for(i in range(2,BF_NB_WORDS+3)){
                        r[i] = ((two_over_pi[i-1]) * x.h_word[1]
                           +(two_over_pi[ i ]) * x.h_word[0]);
                    }
                }else {
                    for(i in range(2,BF_NB_WORDS+3)){
                      r[i] = ((two_over_pi[i-2]) * x.h_word[2]
                         +(two_over_pi[i-1]) * x.h_word[1]
                         +(two_over_pi[ i ]) * x.h_word[0]);
                    }
                }
            }else {
                if (x.index == 2){
                    r[0] = ((two_over_pi[0]) * x.h_word[1]
                     +(two_over_pi[1]) * x.h_word[0]);
                    if(x.h_word[2] == 0){
                        for(i in range(1,BF_NB_WORDS+3)){
                          r[i] = ((two_over_pi[ i ]) * x.h_word[1]
                             +(two_over_pi[i+1]) * x.h_word[0]);
                        }
                    }else {
                        for(i in range(1,BF_NB_WORDS+3)){
                            r[i] = ((two_over_pi[i-1]) * x.h_word[2]
                               +(two_over_pi[ i ]) * x.h_word[1]
                               +(two_over_pi[i+1]) * x.h_word[0]);
                        }
                    }
                }else {
                    var ind = (x.index - 3);
                    var two_over_pi_pt = Pointer[UInt](&two_over_pi[ind]);
                    if(x.h_word[2] == 0){
                        for(i in range(BF_NB_WORDS+3)){
                          r[i] = ((two_over_pi_pt[i+1]) * x.h_word[1]
                             +(two_over_pi_pt[i+2]) * x.h_word[0]);
                        }
                    }else {
                        for(i in range(BF_NB_WORDS+3)){
                            r[i] = wrapAdd(wrapMultiply(two_over_pi_pt[ i ] , x.h_word[2])
                            , wrapMultiply(two_over_pi_pt[i+1] , x.h_word[1])
                            , wrapMultiply(two_over_pi_pt[i+2] , x.h_word[0]));
                        }
                    }
                }
            }
        }
    }

    // Carry propagate
    r[BF_NB_WORDS+1] += bitshr(r[BF_NB_WORDS+2],30);
    for(i in reverseRange((BF_NB_WORDS+1),1)) {
        tmp = bitshr(r[i],30);
        r[i-1] += tmp;
        r[i] -= bitshl(tmp,30);
    }

    // The integer part is in r[0]
    N = UInt(r[0]);

    // test if the reduced part is bigger than Pi/4
    if (r[1] > (BF_RADIX)/2){
        N += 1;
        sign = -1;
        for(i in range(1,(BF_NB_WORDS+3)))
            r[i] = bitand(bitnot(UInt(r[i])) , 0x3fffffff);
    }
    else
        sign = 1;

    // Now we get the reduce argument and check for possible
    // cancellation By Kahan algorithm we will have at most 2 digits
    // of cancellations r[1] and r[2] in the worst case.
    var i = 0;
    if (r[1] == 0)
        if (r[2] == 0)
            i = 3;
        else
            i = 2;
    else
        i = 1;

    for(j in range(BF_NB_WORDS))
        result.h_word[j] = UInt(r[i+j]);

    result.excp   = 1.;
    result.index   = -i;
    result.sign = sign*x.sign;

    // Last step :
    //   Multiplication by pi/2
    result = Pio2 * result;
    return Int(x.sign*N);
}


