ln.biginteger/big.impl/bigint.c

#include <stdint.h>
#include <string.h>

uint32_t one[] = { 1 };

int32_t max(int32_t a,int32_t b)
{
    return a>b ? a : b;
}

int32_t big_get_length(uint32_t *op, int32_t length)
{
    for (; length > 1; length--){
        if (
            ((op[length-1] != 0) || ((op[length-2] & (1<<31))!=0)) && ((op[length-1] != 0xFFFFFFFF) || ((op[length-2] & (1<<31))==0))
        ){
            break;
        }
    }
    return length;
}

int32_t big_sign(uint32_t *op, int32_t length)
{
    if ((length == 1) && (op[0] == 0))
        return 0;
    else if ((op[length-1] & 0x80000000) != 0)
        return -1;
    else
        return 1;
}

int32_t big_log2(uint32_t *op, int32_t length)
{
    if (big_sign(op, length) < 0)
    {
        for (int n = (length << 5)-1; n > 0; n--)
        {
            if ((op[(n - 1) >> 5] & (1 << ((n - 1) & 0x1F))) == 0)
            {
                return 1 - n;
            }
        }
        return 0;
    } else {
        for (int n = length << 5; n > 0; n--)
        {
            if ((op[(n - 1) >> 5] & (1 << ((n - 1) & 0x1F))) != 0)
            {
                return n - 1;
            }
        }
        return 0;
    }
}

int32_t big_cmp(uint32_t* op_a, int32_t length_a, uint32_t* op_b, int32_t length_b)
{
    int l = max(length_a, length_b);
    for (int n=l-1; n >= 0; n--)
    {
        int a = n < length_a ? op_a[n] : 0;
        int b = n < length_b ? op_b[n] : 0;

        if (a < b)
            return -1;
        else if (a > b)
            return 1;
    }
    return 0;
}

void big_shl(uint32_t *op, int32_t length, int32_t n)
{
    if (n > 0)
    {
        int step = n >> 5;
        n &= 0x1F;

        for (int i = length; i > 0; i--)
        {
            int b1 = i - 1 - step;
            int b2 = b1 - 1;

            if (b1 >= 0){
                op[i - 1] = (op[b1] << n);
                if ((n != 0)&&(b2 >= 0))
                {
                    op[i - 1] |= ((op[b2] >> (32 - n)));
                }
            } else {
                op[i - 1] = 0;
            }
        }
    }
}

void big_shr(uint32_t *op, int32_t length, int32_t n)
{
    if (n > 0)
    {
        int32_t s = big_sign(op, length);
        int step = n >> 5;
        n &= 0x1F;

        for (int i = 0; i < length; i++)
        {
            int b1 = i + step;
            int b2 = b1 + 1;

            if (b1 < length)
            {
                op[i] = (op[b1] >> n);
                if ((b2 == length) && (s)){
                    op[i] |= ((0xffffffff << (32 - n)));
                } else if ((n != 0) && (b2 < length))
                {
                    op[i] |= ((op[b2] << (32 - n)));
                }
            } else {
                op[i] = s ? 0xFFFFFFFF : 0;
            }
        }
    }
}


void big_add(uint32_t* op_a, int32_t length_a, uint32_t* op_b, int32_t length_b)
{
    uint64_t c = 0;

    for (int n=0;n<length_a;n++)
    {
        c += op_a[n] + ((n < length_b) ? op_b[n] : 0);
        op_a[n] = (c & 0xffffffff);
        c >>= 32;

        if ((c == 0) && (n > length_b))
            break;
    }
}
void big_sub(uint32_t* op_a, int32_t length_a, uint32_t* op_b, int32_t length_b)
{
    uint64_t c = 0;

    for (int n=0;n<length_a;n++)
    {
        c += op_a[n] - ((n < length_b) ? op_b[n] : 0);
        op_a[n] = (c & 0xffffffff);
        c >>= 32;

        if ((c == 0) && (n > length_b))
            break;
    }
}

void big_twos(uint32_t *op, int32_t length)
{
    for (int n=0;n<length;n++)
        op[n] = ~op[n];
    big_add(op, length, one, 1);
}

void big_smul(uint32_t* op_a, int32_t length_a, uint32_t* op_b, int32_t length_b, uint32_t* result)
{
    int32_t l = length_a + length_b;

    for (int l1 = 0; l1 < length_a; l1++)
    {
        for (int l2 = 0; l2 < length_b; l2++)
        {
            int target = l1 + l2;
            uint64_t ui64 = ((uint64_t)op_a[l1] * op_b[l2]);
            for (int ct = target; (ct < l) && (ui64 != 0); ct++)
            {
                ui64 += result[ct];
                result[ct] = (uint32_t)ui64;
                ui64 >>= 32;
            }
        }
    }
}

void big_divmod(uint32_t* op_a, int32_t length_a, uint32_t* op_b, int32_t length_b)
{
    int sign_a = big_sign(op_a, length_a);
    int sign_b = big_sign(op_b, length_b);

    uint32_t dividend[length_a];
    uint32_t divisor[length_a];

    memset( divisor, 0x00, sizeof(divisor) );

    memcpy( dividend, op_a, length_a << 2);
    memcpy( divisor, op_b, length_b << 2);

    memset( op_a, 0x00, length_a << 2);

    if (sign_a < 0)
        big_twos( dividend, length_a);
    if (sign_b < 0)
        big_twos( divisor, length_a);

    int lg2a, lg2b, shift;

    lg2a = big_log2(dividend, length_a);
    lg2b = big_log2(divisor, length_a);
    shift = lg2a - lg2b;

    if (shift > 0)
        big_shl(divisor, length_a, shift);

    for (int n = 0; n <= (shift); n++)
    {
        big_shl(op_a, length_a, 1);

        if (big_cmp( divisor, length_a, dividend, length_a ) <= 0)
        {
            op_a[0] |= 1;
            big_sub( dividend, length_a, divisor, length_a );
        }

        big_shr( divisor, length_b, 1);
    }

    if ((sign_a < 0) != (sign_b < 0))
        big_twos(op_a, length_a);
    if (sign_a < 0)
        big_twos(dividend, length_a);

    memcpy( op_b, dividend, length_b << 2);
}