// Lightweight BigFloat library based on SCSlib/crlibm
// NOT intended as a replacement for mpfr/gmp for which a Clay wrapper
// will eventually be available as part of the core math lib.
// This is primarily intended for use in math libs to provide high accuracy
// functions and as a high performance 'larger than Float80/128 float type'.

import math.*;
import printer.formatter.(hex);
import numbers.floats.*;

alias BF_NB_BITS = 30;
alias BF_NB_WORDS = 8;
alias BF_MAX_RANGE = 32;
alias BF_RADIX = wrapBitshl(1u,BF_NB_BITS);
alias BF_MASK_RADIX = BF_RADIX-1;

record BigFloat (
    h_word:Array[UInt,BF_NB_WORDS],// the digits, as 32 bits words
    excp:Double, // Used to store Nan,+/-0, Inf, etc
    index:Int,// Exponent in an FP format in base 2^32
    sign:Int, // The sign equals 1 or -1
);

private alias radix_two_double  = floatFromBits(wrapBitshl(1023ul + 2*BF_NB_BITS,52));
private alias radix_mtwo_double = floatFromBits(wrapBitshl(1023ul - 2*BF_NB_BITS,52));
private alias radix_rng_double  = floatFromBits(wrapBitshl(1023ul + BF_NB_BITS*BF_MAX_RANGE,52));

overload BigFloat() = BigFloat(array(0x00000000u, 0x00000000u, 0x00000000u, 0x00000000u,
             0x00000000u, 0x00000000u, 0x00000000u, 0x00000000u),0., 0, 1);

overload BitwiseMovedType?(static BigFloat) = true;

inline overload resetUnsafe(a:BigFloat) {
    a <-- BigFloat();
}

[T] BigFloat?(static T) = false;
inline overload BigFloat?(static BigFloat) = true;


alias bf_zero = BigFloat(array(0x00000000u, 0x00000000u, 0x00000000u, 0x00000000u,
    0x00000000u, 0x00000000u, 0x00000000u, 0x00000000u), 0.,  0,   1 );

alias bf_half = BigFloat(array(0x20000000u, 0x00000000u, 0x00000000u, 0x00000000u,
    0x00000000u, 0x00000000u, 0x00000000u, 0x00000000u), 0., -1,   1 );

alias bf_one = BigFloat(array(0x00000001u, 0x00000000u, 0x00000000u, 0x00000000u,
    0x00000000u, 0x00000000u, 0x00000000u, 0x00000000u), 1.,  0,   1 );

alias bf_two = BigFloat(array(0x00000002u, 0x00000000u, 0x00000000u, 0x00000000u,
    0x00000000u, 0x00000000u, 0x00000000u, 0x00000000u), 1.,  0,   1 );

alias bf_sixinv = BigFloat(array(0x0aaaaaaau, 0x2aaaaaaau, 0x2aaaaaaau, 0x2aaaaaaau,
    0x2aaaaaaau, 0x2aaaaaaau, 0x2aaaaaaau, 0x2aaaaaaau), 1.,  -1,   1 );


overload BigFloat(x:Double) : BigFloat {
    var nb, mantissa = db_number(),db_number();
    var result = BigFloat();
    if(x>=0){result.sign = 1;    nb.2 = x;}
    else    {result.sign = -1;   nb.2 = -x;}

    var exponent = bitand(nb.0[1] , 0x7ff00000);

    if (exponent == 0x7ff00000)  {
        // x = +/- Inf, s/qNAN
        result.excp = x;
        for(i in range(BF_NB_WORDS))
            result.h_word[i] = 0;
        result.index = 0;
        result.sign = 1;
    }
    else {    // Normals,  denormals, +/- 0.
        // This number is not an exception
        result.excp = 1.;
        if (exponent == 0){
            // x is a denormal number : bring it back to the normal range
            nb.2 = nb.2 * radix_two_double;      // 2^(2.BF_NB_BITS)
            exponent = bitand(nb.0[1] , 0x7ff00000);
            result.index = -2;
        }else {
            result.index = 0;
        }

        exponent = bitshr(exponent,20);  // get the actual value

        var ind = ((exponent +(100*BF_NB_BITS - 1023))/BF_NB_BITS) - 100;
        // means : = (exponent -1023 + 100*BF_NB_BITS)/BF_NB_BITS -100
        // The business with 100*BF_NB_BITS is to stay within the positive
        // range for exponent_remainder between 1 and BF_NB_BITS

        var exponent_remainder = exponent - 1022 - (BF_NB_BITS*ind);
        result.index += ind;

        // now get the mantissa and add the implicit 1 in fp. format
        mantissa.1 = bitor(bitand(nb.1 , 0x000ffffffffffffful) , 0x0010000000000000ul);

        // and spread it over the structure
        // Everything here is 64-bit arithmetic
        result.h_word[0] =  UInt32(bitshr(mantissa.1 , 53 - exponent_remainder ));

        // 11 = 64-53
        mantissa.1 = wrapBitshl(mantissa.1 , wrapAdd(exponent_remainder,11));
        result.h_word[1] = bitand(bitshr(mantissa.0[1] , 32 - BF_NB_BITS) , BF_MASK_RADIX );
        mantissa.1 =  wrapBitshl(mantissa.1 , BF_NB_BITS);
        result.h_word[2] = bitand(bitshr(mantissa.0[1] , (32 - BF_NB_BITS)), BF_MASK_RADIX);
        result.h_word[3] = 0 ;
        result.h_word[4] = 0; result.h_word[5] = 0; result.h_word[6] = 0; result.h_word[7] = 0;
    } // end if test NaN etc
    return result;
}

[T | BigFloat?(T)]
overload Double(x:T) --> result:Double {
    var nb, rndcorr = db_number(), db_number();
    var res = 0.;
    // convert the MSB digit into a double, and store it in nb.d
    nb.2 = Double(x.h_word[0]);
    // place the two next digits in lowpart
    var t1 = UInt64(x.h_word[1]);
    var lowpart = wrapBitshl(t1,BF_NB_BITS) + x.h_word[2];
    // there is at least one significant bit in nb,
    //   and at least 2*BF_NB_BITS in lowpart,
    //   so provided BF_NB_BITS >= 27
    //   together they include the 53+ guard bits to decide rounding

    // test for  s/qNan, +/- Inf, +/- 0, placed here for obscure performance reasons
    if (x.excp != 1.)
        result <-- x.excp;
    // take the exponent of nb.d (will be in [0:BF_NB_BITS])
    var expo = bitshr(bitand(nb.0[1] , 0x7ff00000),20) - 1023;
    // compute the exponent of the result
    var expofinal = expo + BF_NB_BITS*x.index;
    // Is the number not too large for the IEEE exponent range ?
    if (expofinal >  1023)
        // return an infinity
        res = radix_rng_double*radix_rng_double;
    // Is number a denormal  ?
    else if (expofinal >= -1022){
        // x is in the normal range
        var shift = expo+2*BF_NB_BITS-53;
        var roundbits = wrapBitshl(lowpart,64-shift);
        // align the rest of the mantissa to nb : shift by (2*BF_NB_BITS)-53+expo
        lowpart = bitshr(lowpart,shift);
        // Look at the last bit to decide rounding
        if (bitand(lowpart , 0x0000000000000001ul)!=0){
            // Test for the round to even case
            if(roundbits==0)
                for (i in range(3,BF_NB_WORDS))
                    roundbits = bitor(roundbits , x.h_word[i]);
            if(roundbits == 0) {
                // round to even mantissa
                if (bitand(lowpart , 0x0000000000000002ul)!=0){
                    // mantissa odd, need to add an half-ulp
                    rndcorr.0[0] = 0;
                    rndcorr.0[1] = wrapBitshl(expo-52+1023,20);    // 2^(exp-52)
                }
                else
                    rndcorr.2 = 0.0;
            }
            else { // there are round bits need to add an half-ulp
                rndcorr.0[0] = 0;
                rndcorr.0[1] = wrapBitshl(expo-52+1023,20);    // 2^(exp-52)
            }
        }
        else{
            // need to add nothing
            rndcorr.2 = 0.0;
        }
        lowpart = bitshr(lowpart,1);
        nb.1 = bitor(nb.1,lowpart); // Finish to fill the mantissa
        res  = nb.2 + rndcorr.2;  // rounded to nearest
        // now compute the exponent from the index :
        // we need to multiply res by 2^(X_IND*BF_NB_BITS)
        // First check this number won't be a denormal itself
        if(x.index*BF_NB_BITS +1023 > 0) {
            // build the double 2^(X_IND*BF_NB_BITS)
            nb.0[1] = wrapBitshl(x.index*BF_NB_BITS +1023,20);
            nb.0[0] = 0;
            res *= nb.2;     // exact multiplication
        }
        else {
            // offset the previous computation by 2^(2*BF_NB_BITS)
            // build the double 2^(X_IND*BF_NB_BITS)
            nb.0[1] = wrapBitshl(x.index*BF_NB_BITS +1023 + 2*BF_NB_BITS,20);
            nb.0[0] = 0;
            res *= radix_mtwo_double;  // exact multiplication
            res *= nb.2;                  // exact multiplication
        }

    }
    else {
        // the final number is a denormal with 52-(expfinal+1022)
        // significant bits.
        if (expofinal < -1022 - 53 )
            res = 0.;
        else {
            // align the rest of the mantissa to nb
            lowpart = bitshr(lowpart , expo+(2*BF_NB_BITS)-52);
            // Finish to fill the mantissa
            nb.1 = bitor(nb.1,lowpart);
            // this is still a normal number.
            //Now remove its exponent and add back the implicit one
            nb.1 = bitor(bitand(nb.1,0x000FFFFFFFFFFFFFul), 0x0010000000000000ul);
            // keep only the significant bits
            nb.1 = bitshr(nb.1,-1023 - expofinal);
            // Look at the last bit to decide rounding
            if (bitand(nb.0[0] , 0x00000001)==1){
                // need to add an half-ulp */
                rndcorr.1 = 1;  // this is a full ulp but we multiply by 0.5 in the end
            }else{
                // need to add nothing
                rndcorr.2 = 0.0;
            }
            res  = 0.5*(nb.2 + rndcorr.2);  // rounded to nearest
        }
    }
    // sign management
    if (x.sign < 0)
        result <-- -res;
    else
        result <-- res;
}


private inline BF_CARRY_PROPAGATE(r1,r0,tmp) {
    tmp = bitshr(UInt64(r1),BF_NB_BITS);
    r0 += tmp;
    r1 -= wrapBitshl(tmp,BF_NB_BITS);
}


// This function copies a result into another.
[T | BigFloat?(T)]
inline overload assign(ref result:T, x:T) : {
    result.hword = x.hword;
    result.excp = x.excp;
    result.index = x.index;
    result.sign = x.sign;
}

[T | Numeric?(T)]
inline overload assign(ref result:BigFloat, x:T) {
    if(x == 0){
        resetUnsafe(result.h_word);
        result.excp = 0.;
        result.index = 0;
        result.sign = 1;
    }
    else{
        result = BigFloat(x);
    }
}

[T | BigFloat?(T)]
overload multiply(x:T, y:T) --> result:T {
    result.excp = x.excp * y.excp;
    result.sign = x.sign * y.sign;
    result.index = x.index + y.index;

    /* Partial products computation */
    var x7,y7,x6,y6 = UInt64(x.h_word[7]),y.h_word[7],UInt64(x.h_word[6]),y.h_word[6];
    var x5,y5,x4,y4 = UInt64(x.h_word[5]),y.h_word[5],UInt64(x.h_word[4]),y.h_word[4];
    var x3,y3,x2,y2 = UInt64(x.h_word[3]),y.h_word[3],UInt64(x.h_word[2]),y.h_word[2];
    var x1,y1,x0,y0 = UInt64(x.h_word[1]),y.h_word[1],UInt64(x.h_word[0]),y.h_word[0];

    var r8 = x7*y1 + x6*y2 + x5*y3 + x4*y4 + x3*y5 + x2*y6 + x1*y7;
    var r7 = x7*y0 + x6*y1 + x5*y2 + x4*y3 + x3*y4 + x2*y5 + x1*y6 + x0*y7;
    var r6 = x6*y0 + x5*y1 + x4*y2 + x3*y3 + x2*y4 + x1*y5 + x0*y6;
    var r5 = x5*y0 + x4*y1 + x3*y2 + x2*y3 + x1*y4 + x0*y5;
    var r4 = x4*y0 + x3*y1 + x2*y2 + x1*y3 + x0*y4;
    var r3 = x3*y0 + x2*y1 + x1*y2 + x0*y3;
    var r2 = x2*y0 + x1*y1 + x0*y2;
    var r1 = x1*y0 + x0*y1;
    var r0 = x0*y0;

    var val,tmp = 0ul,0ul;
    /* Carry Propagate */
    BF_CARRY_PROPAGATE(r8,r7,tmp);
    BF_CARRY_PROPAGATE(r7,r6,tmp);
    BF_CARRY_PROPAGATE(r6,r5,tmp);
    BF_CARRY_PROPAGATE(r5,r4,tmp);
    BF_CARRY_PROPAGATE(r4,r3,tmp);
    BF_CARRY_PROPAGATE(r3,r2,tmp);
    BF_CARRY_PROPAGATE(r2,r1,tmp);
    BF_CARRY_PROPAGATE(r1,r0,tmp);
    BF_CARRY_PROPAGATE(r0,val,tmp);

    if(val != 0){
        /* shift all the digits ! */
        result.h_word[0] = UInt(val); result.h_word[1] = UInt(r0); result.h_word[2] = UInt(r1);  result.h_word[3] = UInt(r2);
        result.h_word[4] = UInt(r3);  result.h_word[5] = UInt(r4); result.h_word[6] = UInt(r5);  result.h_word[7] = UInt(r6);
        inc(result.index);
    }
    else {
        result.h_word[0] = UInt(r0); result.h_word[1] = UInt(r1); result.h_word[2] = UInt(r2); result.h_word[3] = UInt(r3);
        result.h_word[4] = UInt(r4); result.h_word[5] = UInt(r5); result.h_word[6] = UInt(r6); result.h_word[7] = UInt(r7);
    }
}


//Square BigFloat
[T | BigFloat?(T)]
square(x:T) : T {
    var result = BigFloat();
    result.excp = x.excp * x.excp;
    result.index = x.index + x.index;
    result.sign = 1;

    // Partial products
    var x7,x6,x5,x4 = UInt64(x.h_word[7]),UInt64(x.h_word[6]),UInt64(x.h_word[5]),UInt64(x.h_word[4]);
    var x3,x2,x1,x0 = UInt64(x.h_word[3]),UInt64(x.h_word[2]),UInt64(x.h_word[1]),UInt64(x.h_word[0]);

    var r0 =  x0*x0;
    var r1 = (x0*x1)* 2. ;
    var r2 =  x1*x1 + (x0*x2*2.);
    var r3 = (x1*x2 +  x0*x3)* 2.;
    var r4 =  x2*x2 + (x1*x3 + x0*x4)* 2.;
    var r5 = (x2*x3 +  x1*x4 + x0*x5)* 2.;
    var r6 =  x3*x3 + (x2*x4 + x1*x5 + x0*x6)* 2.;
    var r7 = (x3*x4 +  x2*x5 + x1*x6 + x0*x7)* 2.;
    var r8 =  x4*x4 + (x3*x5 + x2*x6 + x1*x7)* 2.;

    var val,tmp = 0ul,0ul;
    // Carry propagation
    BF_CARRY_PROPAGATE(r8,r7,tmp);
    BF_CARRY_PROPAGATE(r7,r6,tmp);
    BF_CARRY_PROPAGATE(r6,r5,tmp);
    BF_CARRY_PROPAGATE(r5,r4,tmp);
    BF_CARRY_PROPAGATE(r4,r3,tmp);
    BF_CARRY_PROPAGATE(r3,r2,tmp);
    BF_CARRY_PROPAGATE(r2,r1,tmp);
    BF_CARRY_PROPAGATE(r1,r0,tmp);
    BF_CARRY_PROPAGATE(r0,val,tmp);

    if(val != 0){
        // shift all the digits !
        result.h_word[0] = UInt(val); result.h_word[1] = UInt(r0); result.h_word[2] = UInt(r1);  result.h_word[3] = UInt(r2);
        result.h_word[4] = UInt(r3);  result.h_word[5] = UInt(r4); result.h_word[6] = UInt(r5);  result.h_word[7] = UInt(r6);
        inc(result.index);
    }
    else {
        result.h_word[0] = UInt(r0); result.h_word[1] = UInt(r1); result.h_word[2] = UInt(r2); result.h_word[3] = UInt(r3);
        result.h_word[4] = UInt(r4); result.h_word[5] = UInt(r5); result.h_word[6] = UInt(r6); result.h_word[7] = UInt(r7);
    }
    return result;
}



//  The function that does the work in case of an addition
//
//  do_add is the function that does the addition of two BigFloat numbers,
//    assuming that x.sign == y.sign, X_IND > Y_IND, x and y both
//    non-zero.
//
private do_add(result, x, y) : {
    var r0,r1,r2,r3,r4,r5,r6,r7 = 0,0,0,0,0,0,0,0;
    var dif = x.index - y.index;
    result.excp = x.excp + y.excp - 1.;
    result.index = x.index;
    result.sign = x.sign;

    switch (dif)
        case (0){
            r0 = x.h_word[0] + y.h_word[0]; r1 = x.h_word[1] + y.h_word[1];
            r2 = x.h_word[2] + y.h_word[2]; r3 = x.h_word[3] + y.h_word[3];
            r4 = x.h_word[4] + y.h_word[4]; r5 = x.h_word[5] + y.h_word[5];
            r6 = x.h_word[6] + y.h_word[6]; r7 = x.h_word[7] + y.h_word[7];
        }
        case (1){
            r0 = x.h_word[0];           r1 = x.h_word[1] + y.h_word[0];
            r2 = x.h_word[2] + y.h_word[1]; r3 = x.h_word[3] + y.h_word[2];
            r4 = x.h_word[4] + y.h_word[3]; r5 = x.h_word[5] + y.h_word[4];
            r6 = x.h_word[6] + y.h_word[5]; r7 = x.h_word[7] + y.h_word[6];
        }
        case (2){
            r0 = x.h_word[0];           r1 = x.h_word[1];
            r2 = x.h_word[2] + y.h_word[0]; r3 = x.h_word[3] + y.h_word[1];
            r4 = x.h_word[4] + y.h_word[2]; r5 = x.h_word[5] + y.h_word[3];
            r6 = x.h_word[6] + y.h_word[4]; r7 = x.h_word[7] + y.h_word[5];
        }
        case (3){
            r0 = x.h_word[0];           r1 = x.h_word[1];
            r2 = x.h_word[2];           r3 = x.h_word[3] + y.h_word[0];
            r4 = x.h_word[4] + y.h_word[1]; r5 = x.h_word[5] + y.h_word[2];
            r6 = x.h_word[6] + y.h_word[3]; r7 = x.h_word[7] + y.h_word[4];
        }
        case (4){
            r0 = x.h_word[0];           r1 = x.h_word[1];
            r2 = x.h_word[2];           r3 = x.h_word[3];
            r4 = x.h_word[4] + y.h_word[0]; r5 = x.h_word[5] + y.h_word[1];
            r6 = x.h_word[6] + y.h_word[2]; r7 = x.h_word[7] + y.h_word[3];
        }
        case (5){
            r0 = x.h_word[0];           r1 = x.h_word[1];
            r2 = x.h_word[2];           r3 = x.h_word[3];
            r4 = x.h_word[4];           r5 = x.h_word[5] + y.h_word[0];
            r6 = x.h_word[6] + y.h_word[1]; r7 = x.h_word[7] + y.h_word[2];
        }
        case (6){
            r0 = x.h_word[0];           r1 = x.h_word[1];
            r2 = x.h_word[2];           r3 = x.h_word[3];
            r4 = x.h_word[4];           r5 = x.h_word[5];
            r6 = x.h_word[6] + y.h_word[0]; r7 = x.h_word[7] + y.h_word[1];
        }
        case (7){
            r0 = x.h_word[0];           r1 = x.h_word[1];
            r2 = x.h_word[2];           r3 = x.h_word[3];
            r4 = x.h_word[4];           r5 = x.h_word[5];
            r6 = x.h_word[6];           r7 = x.h_word[7] + y.h_word[0];
        }
        else{// dif >= 8
            result.h_word = x.h_word;
        }

    // Carry propagation
    var carry = bitshr(r7 , BF_NB_BITS); r6 += carry;  r7 = bitand(r7 , BF_MASK_RADIX);
    carry = bitshr(r6 , BF_NB_BITS); r5 += carry;  r6 = bitand(r6 , BF_MASK_RADIX);
    carry = bitshr(r5 , BF_NB_BITS); r4 += carry;  r5 = bitand(r5 , BF_MASK_RADIX);
    carry = bitshr(r4 , BF_NB_BITS); r3 += carry;  r4 = bitand(r4 , BF_MASK_RADIX);
    carry = bitshr(r3 , BF_NB_BITS); r2 += carry;  r3 = bitand(r3 , BF_MASK_RADIX);
    carry = bitshr(r2 , BF_NB_BITS); r1 += carry;  r2 = bitand(r2 , BF_MASK_RADIX);
    carry = bitshr(r1 , BF_NB_BITS); r0 += carry;  r1 = bitand(r1 , BF_MASK_RADIX);
    carry = bitshr(r0 , BF_NB_BITS);

    if (carry!=0){
        result.h_word[7] = r6; result.h_word[6] = r5;  result.h_word[5] = r4; result.h_word[4] = r3;
        result.h_word[3] = r2; result.h_word[2] = r1;  result.h_word[1] = bitand(r0 , BF_MASK_RADIX);
        result.h_word[0] = 1 ;
        inc(result.index);
    }
    else {
        result.h_word[0] = r0; result.h_word[1] = r1; result.h_word[2] = r2; result.h_word[3] = r3;
        result.h_word[4] = r4; result.h_word[5] = r5; result.h_word[6] = r6; result.h_word[7] = r7;
    }
}


// SUBTRACTION
//
// This procedure assumes :
// -    X_IND >= Y_IND
// -   X_SIGN != Y_SIGN
//    neither x or y is zero
//  and result = x - y

private do_sub(result, x, y){
    var res = Array[Int,BF_NB_WORDS]();
    result.excp = x.excp + y.excp - 1.;
    var dif  = x.index - y.index;
    result.index = x.index;

    // The easy case
    if(dif >= BF_NB_WORDS)
        result = x;
    else { // 0 <= dif <= (BF_NB_WORDS-1)
        var carry = 0;
        if(dif == 0) {
            var i,cp = 0,0;
            while(x.h_word[i] == y.h_word[i] and i<BF_NB_WORDS)
                inc(i);
            if (x.h_word[i] > y.h_word[i])
                cp = 1;
            else if (x.h_word[i] < y.h_word[i])
                cp = -1;
            else
                cp = 0;
            if (cp == 0) {
                result = bf_zero;
                return;
            }
            else { // cp <> 0
                if (cp > 0){
                    // x > y
                    result.sign = x.sign;
                    for(i in reverseRange(BF_NB_WORDS-1,0)){
                        var s = wrapSubtract(wrapSubtract(Int(x.h_word[i]) , Int(y.h_word[i])) , carry);
                        carry = Int(bitshr(bitand(s,Int(BF_RADIX)),BF_NB_BITS));
                        res[i] = wrapAdd(bitand(s,Int(BF_RADIX)), s);
                    }
                }
                else {
                    result.sign = - x.sign;
                    for(i in reverseRange(BF_NB_WORDS-1,0)){
                        var s = Int(wrapAdd(-Int(x.h_word[i]), Int(y.h_word[i])) - carry);
                        carry = bitshr(bitand(s,Int(BF_RADIX)),BF_NB_BITS);
                        res[i] = wrapAdd(bitand(s,Int(BF_RADIX)), s);
                    }
                }
            }
        }
        else {
            // 1<=dif<(BF_NB_WORDS-1) Digits of x and y overlap but the
            // sign will be that of x
            result.sign = x.sign;
            var s,i,j = 0,(BF_NB_WORDS-1),(BF_NB_WORDS-1)-dif;
            while(i >= 0){
                if(j >= 0)
                    s = wrapSubtract(wrapSubtract(Int(x.h_word[i]) , Int(y.h_word[j])) , carry);
                else
                    s = (x.h_word[i] - carry);
                carry = bitshr(bitand(s,BF_RADIX),BF_NB_BITS);
                res[i] = wrapAdd(Int(bitand(s,BF_RADIX)), s);
                dec(i);
                dec(j);
            }
        }
        // check for cancellations
        var i = 0;
        while (res[i] == 0 and i < BF_NB_WORDS)
            inc(i);
        if(i > 0) { // cancellation, shift result
            result.index -= i;
            var j = 0;
            while(i < BF_NB_WORDS){
                result.h_word[j] = res[i];
                inc(i);
                inc(j);
            }
            while(j < BF_NB_WORDS){
                result.h_word[j] = 0;
                inc(j);
            }
        }
        else {
            for(i in range(BF_NB_WORDS))
                result.h_word[i] = UInt(res[i]);
        }
    }

}


// Addition (result is a normalised number).
[T | BigFloat?(T)]
overload add(x:T, y:T) : T {
    var result = BigFloat();
    if (x.excp == 0)
        result = y;
    else if (y.excp == 0)
        result = x;
    else if (x.sign == y.sign)
        if(x.index >= y.index)
            do_add(result,x,y);
        else
            do_add(result,y,x);
    else
        if(x.index >= y.index)
            do_sub(result,x,y);
        else
            do_sub(result,y,x);
    return result;
}


// Subtraction (result is a normalised number).
// The arguments x, y and result may point to the same memory location.
[T | BigFloat?(T)]
overload subtract(x:T, y:T) :T {
    var result = BigFloat();
    if (x.excp == 0){
        result = y;
        result.sign = -result.sign;
    }
    else if (y.excp == 0)
        result = x;
    else if (x.sign == y.sign){
        // Same sign, so it's a sub
        if(x.index >= y.index)
            do_sub(result,x,y);
        else{
            do_sub(result,y,x);
            result.sign = -result.sign;
        }
    }
    else{
        if(x.index >= y.index)
            do_add(result,x,y);
        else{
            do_add(result,y,x);
            result.sign = -result.sign;
        }
    }
    return result;
}

// Compute 1/x with a Newton scheme
[T | BigFloat?(T)]
overload recip(x:T) --> result:T {
    var tmp = x;
    tmp.index = 0;
    var app_x = Double(tmp);
    // inv is a 53-bit approximation of 1/x
    var inv = 1./app_x;
    var res = BigFloat(inv);
    res.index -= x.index;

    // First Newton Iteration
    var res1 = x * res;
    res1 = bf_two - res1;
    res *= res1;

    // Second Newton Iteration
    res1 = x * res;
    res1 = bf_two - res1;
    result <-- res * res1;
}

// Compute result = x/y;
[T | BigFloat?(T)]
overload divide(x:T, y:T) :T {
    var result = BigFloat();
    if (x.excp != 1.)
        result.excp = x.excp / y.excp;
    else{
        var res = recip(y);
        result <-- res * x;
    }
    return result;
}

// Convert an integer number in multiprecision representation
[I | Integer?(I)]
overload BigFloat(x:I) :BigFloat {
    var result = BigFloat();
    var ux = 0u;
    if(x >= 0){
        result.sign = 1;
        ux = UInt(x);
    }
    else{
        result.sign = -1;
        ux = UInt(-x);
    }
    if (ux > BF_RADIX){
        result.index = 1;
        result.h_word[0] = bitshr(ux - BF_RADIX, BF_NB_BITS);
        result.h_word[1] =  ux -  wrapBitshl(result.h_word[0] , BF_NB_BITS);
    }else {
        result.index = 0;
        result.h_word[0] = ux;
        result.h_word[1] = 0;
    }
    for(i in range(2,BF_NB_WORDS))
        result.h_word[i] = 0;

    if (x != 0)
        result.excp = 1;
    else
        result.excp = 0;
    return result;
}

overload printTo(stream, x:BigFloat){
    var d = db_number();
    d.2 = Double(x);
    println("Exception : ", x.excp);
    println("Index= ",x.index ,"   \n Sign=  ", x.sign ," \n Double value= ",d.2,"   \n Hex mantissa= ", d.0[1],"  ", d.0[0]);
    for(i in range(BF_NB_WORDS))
        println("  D ", i," : ", hex(x.h_word[i]),"  ", x.h_word[i]);
}
