289 lines
6.7 KiB
C
289 lines
6.7 KiB
C
#include <stdint.h>
|
|
#include <string.h>
|
|
#include <stdio.h>
|
|
|
|
uint32_t one[] = { 1 };
|
|
|
|
int32_t max(int32_t a,int32_t b)
|
|
{
|
|
return a>b ? a : b;
|
|
}
|
|
|
|
int32_t big_is_zero( uint32_t* op, int32_t length)
|
|
{
|
|
for (int n=0;n<length;n++)
|
|
if (op[n])
|
|
return 0;
|
|
return -1;
|
|
}
|
|
|
|
#define DUMP(op,l) for (int ___p = l-1; ___p >= 0; ___p--) printf("%08x", op[___p]); printf("\n");
|
|
|
|
int32_t big_get_bitlength(uint32_t *op, int32_t length)
|
|
{
|
|
int32_t sign = (op[ length - 1 ] >> 31) & 0x00000001;
|
|
for (int ni = length - 1; ni >= 0; ni--)
|
|
for (int nb = 31; nb >= 0; nb--)
|
|
{
|
|
if (((op[ni] >> nb) & 0x00000001) != sign)
|
|
return (ni << 5) + nb + 2;
|
|
}
|
|
return 1;
|
|
}
|
|
|
|
int32_t big_get_length(uint32_t *op, int32_t length)
|
|
{
|
|
return (big_get_bitlength(op, length) + 31) >> 5;
|
|
}
|
|
|
|
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) & 0x1FFFFFFF))) == 0)
|
|
// {
|
|
// return 1 - n;
|
|
// }
|
|
// }
|
|
// return 0;
|
|
// } else {
|
|
// for (int n = length << 5; n > 0; n--)
|
|
// {
|
|
// if ((op[(n - 1) >> 5] & (1 << ((n - 1) & 0x1FFFFFFF))) != 0)
|
|
// {
|
|
// return n - 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 < 0)){
|
|
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)
|
|
{
|
|
int64_t c = 0;
|
|
|
|
for (int n=0;n<length_a;n++)
|
|
{
|
|
c += (int64_t)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)
|
|
{
|
|
int64_t c = 0;
|
|
for (int n=0;n<length_a;n++)
|
|
{
|
|
c += (int64_t)op_a[n] - ((n < length_b) ? op_b[n] : 0);
|
|
op_a[n] = (c & 0xffffffff);
|
|
c >>= 32;
|
|
|
|
if ((c == 0) && (n > length_b))
|
|
break;
|
|
}
|
|
}
|
|
|
|
int32_t big_cmp(uint32_t* op_a, int32_t length_a, uint32_t* op_b, int32_t length_b)
|
|
{
|
|
int32_t l = max(length_a, length_b);
|
|
uint32_t temp[l];
|
|
memcpy(temp, op_a, length_a << 2);
|
|
big_sub( temp, l, op_b, length_b );
|
|
|
|
if (temp[l-1] & 0x80000000)
|
|
return -1;
|
|
if (big_is_zero( temp, l ))
|
|
return 0;
|
|
return 1;
|
|
}
|
|
|
|
|
|
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 op_a_bak[length_a];
|
|
memcpy( op_a_bak, op_a, length_a << 2);
|
|
memset( op_a, 0x00, length_a << 2);
|
|
|
|
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_bak[l1] * op_b[l2]);
|
|
for (int ct = target; (ct < length_a) && (ui64 != 0); ct++)
|
|
{
|
|
ui64 += op_a[ct];
|
|
op_a[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);
|
|
|
|
int l = max(length_a, length_b);
|
|
|
|
uint32_t dividend[l];
|
|
uint32_t divisor[l];
|
|
|
|
memset( dividend, 0x00, l<<2 );
|
|
memset( divisor, 0x00, l<<2 );
|
|
|
|
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, l);
|
|
if (sign_b < 0)
|
|
big_twos( divisor, l);
|
|
|
|
int lg2a, lg2b, shift;
|
|
|
|
lg2a = big_get_bitlength(dividend, l);
|
|
lg2b = big_get_bitlength(divisor, l);
|
|
shift = lg2a - lg2b;
|
|
|
|
// printf("big_divmod():\n");
|
|
// printf(" dividend = ");
|
|
// DUMP(dividend, l);
|
|
// printf(" divisor = ");
|
|
// DUMP(divisor, l);
|
|
|
|
if (shift > 0)
|
|
big_shl(divisor, l, shift);
|
|
|
|
// printf(" divisor (SHIFTED) = ");
|
|
// DUMP(divisor, l);
|
|
|
|
for (int n = 0; n <= (shift); n++)
|
|
{
|
|
// printf("+++\n");
|
|
// printf(" dividend = ");
|
|
// DUMP(dividend, l);
|
|
// printf(" divisor = ");
|
|
// DUMP(divisor, l);
|
|
|
|
big_shl(op_a, length_a, 1);
|
|
|
|
if (big_cmp( divisor, l, dividend, l ) <= 0)
|
|
{
|
|
op_a[0] |= 1;
|
|
big_sub( dividend, l, divisor, l );
|
|
}
|
|
|
|
// printf(" result = ");
|
|
// DUMP(op_a, length_a);
|
|
|
|
big_shr( divisor, l, 1);
|
|
}
|
|
|
|
if ((sign_a < 0) != (sign_b < 0))
|
|
big_twos(op_a, length_a);
|
|
if (sign_a < 0)
|
|
big_twos(dividend, l);
|
|
|
|
memcpy( op_b, dividend, length_b << 2);
|
|
// printf("----------------------------------------\n");
|
|
}
|
|
|
|
int32_t big_parse(uint32_t* op, int32_t length, char *source, int32_t radix)
|
|
{
|
|
int sl = strlen(source);
|
|
memset( op, 0x00, length << 2);
|
|
for (int n=0;n<sl;n++)
|
|
{
|
|
uint32_t v = source[n] - '0';
|
|
if (v > 9)
|
|
v = source[n] - 55;
|
|
if (v > 36)
|
|
v = source[n] - 87;
|
|
|
|
big_smul(op, length, (uint32_t*)&radix, 1);
|
|
big_add(op, length, (uint32_t*)&v, 1);
|
|
}
|
|
return 0;
|
|
}
|