Stage simple bignum multiplication

This commit is contained in:
Justin Ethier 2022-06-15 21:50:00 -04:00
parent 4833788223
commit 5826b53607
3 changed files with 88 additions and 27 deletions

View file

@ -866,29 +866,29 @@ typedef struct {
#define C_WORD_SIZE 32
#define C_HALF_WORD_SIZE 16
// #endif
//
// /* Tunable performance-related constants */
// #ifndef C_KARATSUBA_THRESHOLD
// /* This defines when we'll switch from schoolbook to Karatsuba
// * multiplication. The smallest of the two numbers determines the
// * switch. It is pretty high right now because it generates a bit
// * more garbage and GC overhead dominates the algorithmic performance
// * gains. If the GC is improved, this can be readjusted.
// */
// # define C_KARATSUBA_THRESHOLD 70
// #endif
// #ifndef C_BURNIKEL_ZIEGLER_THRESHOLD
// /* This defines when to switch from schoolbook to Burnikel-Ziegler
// * division. It creates even more garbage than Karatsuba :(
// */
// # define C_BURNIKEL_ZIEGLER_THRESHOLD 300
// #endif
// #ifndef C_RECURSIVE_TO_STRING_THRESHOLD
// /* This threshold is in terms of the expected string length. */
// # define C_RECURSIVE_TO_STRING_THRESHOLD 750
// #endif
//
// /* These might fit better in runtime.c? */
/* Tunable performance-related constants */
#ifndef C_KARATSUBA_THRESHOLD
/* This defines when we'll switch from schoolbook to Karatsuba
* multiplication. The smallest of the two numbers determines the
* switch. It is pretty high right now because it generates a bit
* more garbage and GC overhead dominates the algorithmic performance
* gains. If the GC is improved, this can be readjusted.
*/
# define C_KARATSUBA_THRESHOLD 70
#endif
//#ifndef C_BURNIKEL_ZIEGLER_THRESHOLD
///* This defines when to switch from schoolbook to Burnikel-Ziegler
// * division. It creates even more garbage than Karatsuba :(
// */
//# define C_BURNIKEL_ZIEGLER_THRESHOLD 300
//#endif
//#ifndef C_RECURSIVE_TO_STRING_THRESHOLD
///* This threshold is in terms of the expected string length. */
//# define C_RECURSIVE_TO_STRING_THRESHOLD 750
//#endif
#define Cyc_fitsinbignumhalfdigitp(n) (C_BIGNUM_DIGIT_HI_HALF(n) == 0)
#define C_BIGNUM_DIGIT_LENGTH C_WORD_SIZE
#define C_BIGNUM_HALF_DIGIT_LENGTH C_HALF_WORD_SIZE
@ -901,6 +901,17 @@ typedef struct {
// #define C_MOST_POSITIVE_32_BIT_FIXNUM 0x3fffffff
// #define C_MOST_NEGATIVE_FIXNUM (-C_MOST_POSITIVE_FIXNUM - 1)
/* The bignum digit representation is fullword- little endian, so on
* LE machines the halfdigits are numbered in the same order. On BE
* machines, we must swap the odd and even positions.
*/
//#ifdef C_BIG_ENDIAN
#define C_uhword_ref(x, p) ((uint16_t *)(x))[(p)^1]
//#else
//#define C_uhword_ref(x, p) ((C_uhword *)(x))[(p)]
//#endif
#define C_uhword_set(x, p, d) (C_uhword_ref(x,p) = (d))
/**
* @brief Complex number
*/
@ -1610,6 +1621,7 @@ object Cyc_int2bignum2(gc_thread_data *data, int n);
string_type *bignum2string(void *data, bignum2_type *bn, int base);
object _str_to_bignum(void *data, char *str, char *str_end, int radix);
object bignum2_plus_unsigned(void *data, bignum2_type *x, bignum2_type *y, int negp);
object bignum_times_bignum_unsigned(void *data, object x, object y, int negp);
/* Remaining GC prototypes that require objects to be defined */
void *gc_alloc_from_bignum(gc_thread_data *data, bignum_type *src);

View file

@ -2581,8 +2581,7 @@ static uint32_t bignum_digits_destructive_scale_down(uint32_t *start, uint32_t *
// return carry; /* This would end up as most significant digit if it fit */
//}
// TODO: static C_regparm void
bignum_digits_multiply(object x, object y, object result)
static void bignum_digits_multiply(object x, object y, object result)
{
uint32_t product,
*xd = C_bignum_digits(x),
@ -2598,8 +2597,8 @@ bignum_digits_multiply(object x, object y, object result)
if (yj == 0) continue;
carry = 0;
for (i = 0; i < length_x; ++i) {
product = (C_uword)C_uhword_ref(xd, i) * yj +
(C_uword)C_uhword_ref(rd, i + j) + carry;
product = (uint32_t)C_uhword_ref(xd, i) * yj +
(uint32_t)C_uhword_ref(rd, i + j) + carry;
C_uhword_set(rd, i + j, product);
carry = C_BIGNUM_DIGIT_HI_HALF(product);
}
@ -2607,6 +2606,31 @@ bignum_digits_multiply(object x, object y, object result)
}
}
// TODO: static
object bignum_times_bignum_unsigned(void *data, object x, object y, int negp)
{
object res = boolean_f;
uint32_t size;
if (C_bignum_size(y) < C_bignum_size(x)) { /* Ensure size(x) <= size(y) */
object z = x;
x = y;
y = z;
}
// TODO:
//if (C_bignum_size(x) >= C_KARATSUBA_THRESHOLD)
// res = bignum_times_bignum_karatsuba(ptr, x, y, negp);
if (res == boolean_f) {
size = C_bignum_size(x) + C_bignum_size(y);
res = gc_alloc_bignum2(data, size);
C_bignum_sign(res) = negp;
bignum_digits_multiply(x, y, res);
res = C_bignum_simplify(res);
}
return res;
}
///* "small" is either a number that fits a halfdigit, or a power of two */
//static C_regparm void

View file

@ -32,6 +32,15 @@ if(is_value_type(result)) {
}
return_closcall1(data, k, result);
")
(define-c test-times
"(void *data, int argc, closure _, object k, object fx1, object fx2)"
" object x = Cyc_int2bignum2(data, obj_obj2int(fx1));
object y = Cyc_int2bignum2(data, obj_obj2int(fx2));
int negp = 0;
// TODO: set negp
object result = bignum_times_bignum_unsigned(data, x, y, negp);
return_closcall1(data, k, result);
")
(define-c test-bn
"(void *data, int argc, closure _, object k, object fx)"
" object bn = Cyc_int2bignum2(data, obj_obj2int(fx));
@ -60,6 +69,22 @@ if(is_value_type(result)) {
(write (test-str2bn "123454354534523454243999" 16))
(newline)
(write "multiplication")
(newline)
(map
(lambda (row)
(write row)
(newline))
(list
(test-times 1 1)
(test-times 1 2)
(test-times -1 2)
(test-times #x0FFFffff #x0FFFffff)
(test-times #x2FFFffff #x2FFFffff)
))
(write "general")
(newline)
(map
(lambda (row)
(write row)