diff options
| -rw-r--r-- | ChangeLog | 18 | ||||
| -rw-r--r-- | arith.c | 44 | ||||
| -rw-r--r-- | eval.c | 3 | ||||
| -rw-r--r-- | lib.c | 14 | ||||
| -rw-r--r-- | lib.h | 3 | ||||
| -rw-r--r-- | mpi-patches/faster-square-root | 183 | ||||
| -rw-r--r-- | mpi-patches/series | 1 | ||||
| -rw-r--r-- | txr.1 | 6 |
8 files changed, 269 insertions, 3 deletions
@@ -1,5 +1,23 @@ 2011-12-13 Kaz Kylheku <kaz@kylheku.com> + * arith.c (highest_bit): Linkage changed to static. + (abso, isqrt): New functions. + (isqrt_fixnum): New static function. + + * eval.c (eval_init): Registered abs, sqrt and numberp instrinsics. + + * lib.c (numberp): New function. + + * lib.h (numberp, abso, isqrt): Declared. + + * mpi-patches/series: New patch added. + + * mpi-patches/faster-square-root: New patch added. + + * txr.1: Documentation stubs for new functions. + +2011-12-13 Kaz Kylheku <kaz@kylheku.com> + * arith.c (expt): Fix broken bignum x fixnum combination. 2011-12-13 Kaz Kylheku <kaz@kylheku.com> @@ -84,7 +84,7 @@ static val normalize(val bignum) } } -int highest_bit(int_ptr_t n) +static int highest_bit(int_ptr_t n) { #if SIZEOF_PTR == 8 if (n & 0x7FFFFFFF00000000) { @@ -414,6 +414,18 @@ val neg(val anum) } } +val abso(val anum) +{ + if (bignump(anum)) { + val n = make_bignum(); + mp_abs(mp(anum), mp(n)); + return n; + } else { + cnum n = c_num(anum); + return num(n < 0 ? n : n); + } +} + val mul(val anum, val bnum) { int tag_a = tag(anum); @@ -890,6 +902,36 @@ val expt(val anum, val bnum) abort(); } +static int_ptr_t isqrt_fixnum(int_ptr_t a) +{ + int_ptr_t mask = (int_ptr_t) 1 << (highest_bit(a) / 2); + int_ptr_t root = 0; + + for (; mask != 0; mask >>= 1) { + int_ptr_t next_guess = root | mask; + if (next_guess * next_guess <= a) + root = next_guess; + } + + return root; +} + +val isqrt(val anum) +{ + if (fixnump(anum)) { + cnum a = c_num(anum); + if (a < 0) + uw_throw(error_s, lit("sqrt: negative operand")); + return num_fast(isqrt_fixnum(c_num(anum))); + } else if (bignump(anum)) { + val n = make_bignum(); + if (mp_sqrt(mp(anum), mp(n)) != MP_OKAY) + uw_throw(error_s, lit("sqrt: negative operand")); + return normalize(n); + } + uw_throwf(error_s, lit("sqrt: invalid operand ~s"), anum, nao); +} + void arith_init(void) { mp_init(&NUM_MAX_MP); @@ -1156,11 +1156,14 @@ void eval_init(void) reg_fun(intern(lit("+"), user_package), func_n0v(plusv)); reg_fun(intern(lit("-"), user_package), func_n1v(minusv)); reg_fun(intern(lit("*"), user_package), func_n0v(mulv)); + reg_fun(intern(lit("abs"), user_package), func_n1(abso)); reg_fun(intern(lit("trunc"), user_package), func_n2(trunc)); reg_fun(intern(lit("mod"), user_package), func_n2(mod)); reg_fun(intern(lit("expt"), user_package), func_n0v(exptv)); + reg_fun(intern(lit("sqrt"), user_package), func_n1(isqrt)); reg_fun(intern(lit("fixnump"), user_package), func_n1(fixnump)); reg_fun(intern(lit("bignump"), user_package), func_n1(bignump)); + reg_fun(intern(lit("numberp"), user_package), func_n1(numberp)); reg_fun(intern(lit("zerop"), user_package), func_n1(zerop)); reg_fun(intern(lit(">"), user_package), func_n1v(gtv)); @@ -832,6 +832,20 @@ val bignump(val num) return (is_ptr(num) && type(num) == BGNUM) ? t : nil; } +val numberp(val num) +{ + switch (tag(num)) { + case TAG_NUM: + return t; + case TAG_PTR: + if (num->t.type == BGNUM) + return t; + /* fallthrough */ + default: + return nil; + } +} + val plusv(val nlist) { if (!nlist) @@ -363,11 +363,13 @@ val num(cnum val); cnum c_num(val num); val fixnump(val num); val bignump(val num); +val numberp(val num); val plus(val anum, val bnum); val plusv(val nlist); val minus(val anum, val bnum); val minusv(val minuend, val nlist); val neg(val num); +val abso(val num); val mul(val anum, val bnum); val mulv(val nlist); val trunc(val anum, val bnum); @@ -387,6 +389,7 @@ val maxv(val first, val rest); val minv(val first, val rest); val expt(val base, val exp); val exptv(val nlist); +val isqrt(val anum); val string_own(wchar_t *str); val string(const wchar_t *str); val string_utf8(const char *str); diff --git a/mpi-patches/faster-square-root b/mpi-patches/faster-square-root new file mode 100644 index 00000000..a19bab03 --- /dev/null +++ b/mpi-patches/faster-square-root @@ -0,0 +1,183 @@ +Index: mpi-1.8.6/mpi.c +=================================================================== +--- mpi-1.8.6.orig/mpi.c 2011-12-13 15:18:37.000000000 -0800 ++++ mpi-1.8.6/mpi.c 2011-12-13 17:21:59.000000000 -0800 +@@ -158,6 +158,9 @@ + mp_err s_mp_grow(mp_int *mp, mp_size min); /* increase allocated size */ + mp_err s_mp_pad(mp_int *mp, mp_size min); /* left pad with zeroes */ + ++int s_highest_bit_mp(mp_int *a); ++mp_err s_mp_set_bit(mp_int *a, int bit); ++ + void s_mp_clamp(mp_int *mp); /* clip leading zeroes */ + + void s_mp_exch(mp_int *a, mp_int *b); /* swap a and b in place */ +@@ -1535,77 +1538,55 @@ + + /* {{{ mp_sqrt(a, b) */ + +-/* +- mp_sqrt(a, b) +- +- Compute the integer square root of a, and store the result in b. +- Uses an integer-arithmetic version of Newton's iterative linear +- approximation technique to determine this value; the result has the +- following two properties: +- +- b^2 <= a +- (b+1)^2 >= a +- +- It is a range error to pass a negative value. +- */ + mp_err mp_sqrt(mp_int *a, mp_int *b) + { +- mp_int x, t; +- mp_err res; ++ int mask_shift; ++ mp_int root, guess, *proot = &root, *pguess = &guess; ++ mp_int guess_sqr; ++ mp_err err = MP_MEM; + + ARGCHK(a != NULL && b != NULL, MP_BADARG); + +- /* Cannot take square root of a negative value */ +- if(SIGN(a) == MP_NEG) ++ if (mp_cmp_z(b) == MP_LT) + return MP_RANGE; + +- /* Special cases for zero and one, trivial */ +- if(mp_cmp_d(a, 0) == MP_EQ || mp_cmp_d(a, 1) == MP_EQ) +- return mp_copy(a, b); +- +- /* Initialize the temporaries we'll use below */ +- if((res = mp_init_size(&t, USED(a))) != MP_OKAY) +- return res; +- +- /* Compute an initial guess for the iteration as a itself */ +- if((res = mp_init_copy(&x, a)) != MP_OKAY) +- goto X; +- +- for(;;) { +- /* t = (x * x) - a */ +- mp_copy(&x, &t); /* can't fail, t is big enough for original x */ +- if((res = mp_sqr(&t, &t)) != MP_OKAY || +- (res = mp_sub(&t, a, &t)) != MP_OKAY) +- goto CLEANUP; +- +- /* t = t / 2x */ +- s_mp_mul_2(&x); +- if((res = mp_div(&t, &x, &t, NULL)) != MP_OKAY) +- goto CLEANUP; +- s_mp_div_2(&x); +- +- /* Terminate the loop, if the quotient is zero */ +- if(mp_cmp_z(&t) == MP_EQ) +- break; +- +- /* x = x - t */ +- if((res = mp_sub(&x, &t, &x)) != MP_OKAY) +- goto CLEANUP; +- ++ if ((err = mp_init(&root))) ++ goto out; ++ if ((err = mp_init(&guess))) ++ goto cleanup_root; ++ if ((err = mp_init(&guess_sqr))) ++ goto cleanup_guess; ++ ++ for (mask_shift = s_highest_bit_mp(a) / 2; mask_shift > 0; mask_shift--) { ++ mp_int *temp; ++ ++ if ((err = mp_copy(proot, pguess))) ++ goto cleanup; ++ if ((err = s_mp_set_bit(pguess, mask_shift))) ++ goto cleanup; ++ if ((err = mp_copy(pguess, &guess_sqr))) ++ goto cleanup; ++ if ((err = s_mp_sqr(&guess_sqr))) ++ goto cleanup; ++ ++ if (s_mp_cmp(&guess_sqr, a) <= 0) { ++ temp = proot; ++ proot = pguess; ++ pguess = temp; ++ } + } + +- /* Copy result to output parameter */ +- mp_sub_d(&x, 1, &x); +- s_mp_exch(&x, b); +- +- CLEANUP: +- mp_clear(&x); +- X: +- mp_clear(&t); +- +- return res; +- +-} /* end mp_sqrt() */ ++ err = mp_copy(proot, b); ++ ++cleanup: ++ mp_clear(&guess_sqr); ++cleanup_guess: ++ mp_clear(&guess); ++cleanup_root: ++ mp_clear(&root); ++out: ++ return err; ++} + + /* }}} */ + +@@ -2541,21 +2522,9 @@ + + int mp_count_bits(mp_int *mp) + { +- int len; +- mp_digit d; +- + ARGCHK(mp != NULL, MP_BADARG); + +- len = DIGIT_BIT * (USED(mp) - 1); +- d = DIGIT(mp, USED(mp) - 1); +- +- while(d != 0) { +- ++len; +- d >>= 1; +- } +- +- return len; +- ++ return s_highest_bit_mp(mp); + } /* end mp_count_bits() */ + + /* }}} */ +@@ -3119,6 +3088,27 @@ + abort(); + } + ++int s_highest_bit_mp(mp_int *a) ++{ ++ int nd1 = USED(a) - 1; ++ return s_highest_bit(DIGIT(a, nd1)) + nd1 * MP_DIGIT_BIT; ++} ++ ++mp_err s_mp_set_bit(mp_int *a, int bit) ++{ ++ int nd = (bit + MP_DIGIT_BIT) / MP_DIGIT_BIT; ++ int nbit = bit - (nd - 1) * MP_DIGIT_BIT; ++ mp_err res; ++ ++ if (nd == 0) ++ return MP_OKAY; ++ ++ if ((res = s_mp_pad(a, nd)) != MP_OKAY) ++ return res; ++ ++ DIGIT(a, nd - 1) |= ((mp_digit) 1 << nbit); ++ return MP_OKAY; ++} + + /* {{{ s_mp_exch(a, b) */ + diff --git a/mpi-patches/series b/mpi-patches/series index 311d7fe5..0181c920 100644 --- a/mpi-patches/series +++ b/mpi-patches/series @@ -11,3 +11,4 @@ mpi-set-double-intptr fix-bad-shifts bit-search-optimizations shrink-mpi-int +faster-square-root @@ -4809,9 +4809,11 @@ The following are Lisp functions and variables built-in to TXR. .SS Functions eq, eql and equal -.SS Arithmetic functions +, -, *, trunc, mod, expt +.SS Arithmetic functions +, -, *, trunc, mod, expt, sqrt -.SS Functions fixnump, bignump +.SS Arithmetic function abs + +.SS Functions fixnump, bignump, numberp .SS Function zerop |