~ chicken-core (chicken-5) 98bac9be2fd53bec014f273bca0b389915129cf5
commit 98bac9be2fd53bec014f273bca0b389915129cf5
Author: Peter Bex <peter@more-magic.net>
AuthorDate: Sat Jan 31 16:47:00 2015 +0100
Commit: Peter Bex <peter@more-magic.net>
CommitDate: Sun May 31 14:14:25 2015 +0200
Fix sqrt and signum to accept extended numbers (behaving like CL on cplxnums).
Add exact-integer-sqrt and exact-integer-nth-root for exact root finding.
diff --git a/chicken.h b/chicken.h
index 1b52ef33..bb1fef12 100644
--- a/chicken.h
+++ b/chicken.h
@@ -1276,6 +1276,7 @@ extern double trunc(double);
/* XXX TODO: This should probably be renamed C_u_fixnum_abs or something */
#define C_fixnum_abs(n) C_fix(abs(C_unfix(n)))
#define C_a_i_fixnum_abs(ptr, n, x) (((x) & C_INT_SIGN_BIT) ? C_a_i_fixnum_negate((ptr), (n), (x)) : (x))
+#define C_i_fixnum_signum(x) ((x) == C_fix(0) ? (x) : (((x) & C_INT_SIGN_BIT) ? C_fix(-1) : C_fix(1)))
#define C_i_fixnum_length(x) C_fix(C_ilen(((x) & C_INT_SIGN_BIT) ? ~C_unfix(x) : C_unfix(x)))
#define C_flonum_equalp(n1, n2) C_mk_bool(C_flonum_magnitude(n1) == C_flonum_magnitude(n2))
@@ -1289,6 +1290,7 @@ extern double trunc(double);
#define C_a_i_flonum_times(ptr, c, n1, n2) C_flonum(ptr, C_flonum_magnitude(n1) * C_flonum_magnitude(n2))
#define C_a_i_flonum_quotient(ptr, c, n1, n2) C_flonum(ptr, C_flonum_magnitude(n1) / C_flonum_magnitude(n2))
#define C_a_i_flonum_negate(ptr, c, n) C_flonum(ptr, -C_flonum_magnitude(n))
+#define C_a_u_i_flonum_signum(ptr, n, x) (C_flonum_magnitude(x) == 0.0 ? (x) : ((C_flonum_magnitude(x) < 0.0) ? C_flonum(ptr, -1.0) : C_flonum(ptr, 1.0)))
#define C_a_i_address_to_pointer(ptr, c, addr) C_mpointer(ptr, (void *)C_num_to_unsigned_int(addr))
#define C_a_i_pointer_to_address(ptr, c, pptr) C_unsigned_int_to_num(ptr, (unsigned int)C_c_pointer_nn(pptr))
@@ -1875,6 +1877,7 @@ C_fctimport void C_ccall C_invalid_procedure(int c, C_word self, ...) C_noret;
C_fctexport void C_ccall C_stop_timer(C_word c, C_word closure, C_word k) C_noret;
C_fctexport void C_ccall C_abs(C_word c, C_word self, C_word k, C_word x) C_noret;
C_fctexport void C_ccall C_u_integer_abs(C_word c, C_word self, C_word k, C_word x) C_noret;
+C_fctexport void C_ccall C_signum(C_word c, C_word self, C_word k, C_word x) C_noret;
C_fctexport void C_ccall C_apply(C_word c, C_word closure, C_word k, C_word fn, ...) C_noret;
C_fctexport void C_ccall C_do_apply(C_word n, C_word closure, C_word k) C_noret;
C_fctexport void C_ccall C_call_cc(C_word c, C_word closure, C_word k, C_word cont) C_noret;
@@ -2841,6 +2844,11 @@ C_inline C_word C_i_flonum_max(C_word x, C_word y)
return xf > yf ? x : y;
}
+C_inline C_word C_u_i_integer_signum(C_word x)
+{
+ if (x & C_FIXNUM_BIT) return C_i_fixnum_signum(x);
+ else return (C_bignum_negativep(x) ? C_fix(-1) : C_fix(1));
+}
C_inline C_word
C_a_i_flonum_quotient_checked(C_word **ptr, int c, C_word n1, C_word n2)
diff --git a/chicken.import.scm b/chicken.import.scm
index 47ab3c4e..1e324257 100644
--- a/chicken.import.scm
+++ b/chicken.import.scm
@@ -74,6 +74,8 @@
errno
error
exact-integer?
+ exact-integer-sqrt
+ exact-integer-nth-root
exit
exit-handler
expand
diff --git a/library.scm b/library.scm
index bd5a6cf6..0d84496e 100644
--- a/library.scm
+++ b/library.scm
@@ -319,6 +319,11 @@ EOF
(unless (##core#inline "C_i_exact_integerp" x)
(##sys#error-bad-exact-integer x (and (pair? loc) (car loc))) ) )
+(define (##sys#check-exact-uinteger x . loc)
+ (when (or (not (##core#inline "C_i_exact_integerp" x))
+ (##core#inline "C_i_integer_negativep" x))
+ (##sys#error-bad-exact-uinteger x (and (pair? loc) (car loc))) ) )
+
(define (##sys#check-real x . loc)
(unless (##core#inline "C_i_realp" x)
(##sys#error-bad-real x (and (pair? loc) (car loc))) ) )
@@ -454,6 +459,10 @@ EOF
(##sys#error-hook
(foreign-value "C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR" int) loc arg))
+(define (##sys#error-bad-exact-uinteger arg #!optional loc)
+ (##sys#error-hook
+ (foreign-value "C_BAD_ARGUMENT_TYPE_NO_UINTEGER_ERROR" int) loc arg))
+
(define (##sys#error-bad-inexact arg #!optional loc)
(##sys#error-hook
(foreign-value "C_CANT_REPRESENT_INEXACT_ERROR" int) loc arg))
@@ -1073,10 +1082,13 @@ EOF
(##sys#integer-times
b/g1 (##sys#integer-quotient d g2)))))))))
-(define (signum n)
- (cond ((> n 0) (if (##sys#exact? n) 1 1.0))
- ((< n 0) (if (##sys#exact? n) -1 -1.0))
- (else (if (##sys#exact? n) 0 0.0) ) ) )
+(define (##sys#extended-signum x)
+ (cond
+ ((ratnum? x) (##core#inline "C_u_i_integer_signum" (%ratnum-numerator x)))
+ ((cplxnum? x) (make-polar 1 (angle x)))
+ (else (##sys#error-bad-number x 'signum))))
+
+(define signum (##core#primitive "C_signum"))
(define (flonum->ratnum x)
;; Try to multiply by two until we reach an integer
@@ -1143,8 +1155,8 @@ EOF
;; minimum flonum exponent.
(let* ((s (min (fx- flonum-precision 1)
(fx- e minimum-denorm-flonum-expt)))
- (normalized (##sys#/-2 (arithmetic-shift n s) d))
- (r (round normalized))
+ (norm (##sys#/-2 (##sys#integer-shift n s) d))
+ (r (round norm))
(fraction (exact->inexact r))
(exp (fx- e s)))
(let ((res (fp* fraction (expt 2.0 exp))))
@@ -1152,15 +1164,15 @@ EOF
(scale (lambda (n d) ; Here, 1/2 <= n/d < 2 [N3]
(if (##sys#<-2 n d) ; n/d < 1?
;; Scale left [N3]; only needed once (see note in M3)
- (rnd (arithmetic-shift n 1) d (fx- e 1))
+ (rnd (##sys#integer-shift n 1) d (fx- e 1))
;; Already normalized
(rnd n d e)))))
;; After this step, which shifts the smaller number to
;; align with the larger, "f" in algorithm N is represented
;; in the procedures above by n/d.
(if (negative? e)
- (scale (arithmetic-shift an (##sys#--2 0 e)) d1)
- (scale an (arithmetic-shift d1 e)))))
+ (scale (##sys#integer-shift an (##sys#--2 0 e)) d1)
+ (scale an (##sys#integer-shift d1 e)))))
((cplxnum? x)
(make-complex (exact->inexact (%cplxnum-real x))
(exact->inexact (%cplxnum-imag x))))
@@ -1596,15 +1608,100 @@ EOF
(##core#inline_allocate
("C_a_i_atan" 4) (exact->inexact n)))))
-;; TODO: replace this with an actual sqrt implementation
-(define (##sys#sqrt/loc loc x)
- (##core#inline_allocate ("C_a_i_sqrt" 4) n))
+;; This is "Karatsuba Square Root" as described by Paul Zimmermann,
+;; which is 3/2K(n) + O(n log n) for an input of 2n words, where K(n)
+;; is the number of operations performed by Karatsuba multiplication.
+(define (##sys#exact-integer-sqrt a)
+ ;; Because we assume a3b+a2 >= b^2/4, we must check a few edge cases:
+ (if (and (fixnum? a) (fx<= a 4))
+ (case a
+ ((0 1) (values a 0))
+ ((2) (values 1 1))
+ ((3) (values 1 2))
+ ((4) (values 2 0))
+ (else (error "this should never happen")))
+ (let*-values
+ (((len/4) (fxshr (fx+ (integer-length a) 1) 2))
+ ((len/2) (fxshl len/4 1))
+ ((s^ r^) (##sys#exact-integer-sqrt
+ (##sys#integer-shift a (fxneg len/2))))
+ ((mask) (##sys#--2 (##sys#integer-shift 1 len/4) 1))
+ ((a0) (##sys#integer-bitwise-and a mask))
+ ((a1) (##sys#integer-bitwise-and
+ (##sys#integer-shift a (fxneg len/4)) mask))
+ ((q u) (##sys#integer-quotient&remainder
+ (##sys#+-2 (arithmetic-shift r^ len/4) a1)
+ (##sys#integer-shift s^ 1)))
+ ((s) (##sys#+-2 (##sys#integer-shift s^ len/4) q))
+ ((r) (##sys#+-2 (##sys#integer-shift u len/4)
+ (##sys#--2 a0 (##sys#*-2 q q)))))
+ (if (negative? r)
+ (values (##sys#--2 s 1)
+ (##sys#--2 (##sys#+-2 r (##sys#integer-shift s 1)) 1))
+ (values s r)))))
+
+(define (exact-integer-sqrt x)
+ (##sys#check-exact-uinteger x 'exact-integer-sqrt)
+ (##sys#exact-integer-sqrt x))
+
+;; This procedure is so large because it tries very hard to compute
+;; exact results if at all possible.
+(define (##sys#sqrt/loc loc n)
+ (cond ((cplxnum? n) ; Must be checked before we call "negative?"
+ (let ((p (##sys#/-2 (angle n) 2))
+ (m (##core#inline_allocate ("C_a_i_sqrt" 4) (magnitude n))) )
+ (make-complex (##sys#*-2 m (cos p)) (##sys#*-2 m (sin p)) ) ))
+ ((negative? n)
+ (make-complex .0 (##core#inline_allocate
+ ("C_a_i_sqrt" 4) (exact->inexact (##sys#negate n)))))
+ ((exact-integer? n)
+ (receive (s^2 r) (##sys#exact-integer-sqrt n)
+ (if (eq? 0 r)
+ s^2
+ (##core#inline_allocate ("C_a_i_sqrt" 4) (exact->inexact n)))))
+ ((ratnum? n) ; Try to compute exact sqrt (we already know n is positive)
+ (receive (ns^2 nr) (##sys#exact-integer-sqrt (%ratnum-numerator n))
+ (if (eq? nr 0)
+ (receive (ds^2 dr)
+ (##sys#exact-integer-sqrt (%ratnum-denominator n))
+ (if (eq? dr 0)
+ (##sys#/-2 ns^2 ds^2)
+ (##sys#sqrt/loc loc (exact->inexact n))))
+ (##sys#sqrt/loc loc (exact->inexact n)))))
+ (else (##core#inline_allocate ("C_a_i_sqrt" 4) (exact->inexact n)))))
(define (sqrt x) (##sys#sqrt/loc 'sqrt x))
-;; TODO: unimplemented
+(define (exact-integer-nth-root k n)
+ (##sys#check-exact-uinteger k 'exact-integer-nth-root)
+ (##sys#check-exact-uinteger n 'exact-integer-nth-root)
+ (##sys#exact-integer-nth-root/loc 'exact-integer-nth-root k n))
+
+;; Generalized Newton's algorithm for positive integers, with a little help
+;; from Wikipedia ;) https://en.wikipedia.org/wiki/Nth_root_algorithm
(define (##sys#exact-integer-nth-root/loc loc k n)
- (error "not yet implemented"))
+ (if (or (eq? 0 k) (eq? 1 k) (eq? 1 n)) ; Maybe call exact-integer-sqrt on n=2?
+ (values k 0)
+ (let ((len (integer-length k)))
+ (if (##sys#<-2 len n) ; Idea from Gambit: 2^{len-1} <= k < 2^{len}
+ (values 1 (##sys#--2 k 1)) ; Since x >= 2, we know x^{n} can't exist
+ ;; Set initial guess to (at least) 2^ceil(ceil(log2(k))/n)
+ (let* ((shift-amount (inexact->exact (ceiling (/ (fx+ len 1) n))))
+ (g0 (arithmetic-shift 1 shift-amount))
+ (n-1 (##sys#--2 n 1)))
+ (let lp ((g0 g0)
+ (g1 (quotient
+ (##sys#+-2
+ (##sys#*-2 n-1 g0)
+ (quotient k (##sys#integer-power g0 n-1)))
+ n)))
+ (if (##sys#<-2 g1 g0)
+ (lp g1 (quotient
+ (##sys#+-2
+ (##sys#*-2 n-1 g1)
+ (quotient k (##sys#integer-power g1 n-1)))
+ n))
+ (values g0 (##sys#--2 k (##sys#integer-power g0 n))))))))))
(define (##sys#integer-power base e)
(define (square x) (##sys#*-2 x x))
@@ -1614,7 +1711,7 @@ EOF
(cond
((eq? e2 0) res)
((even? e2) ; recursion is faster than iteration here
- (##sys#*-2 res (square (lp 1 (arithmetic-shift e2 -1)))))
+ (##sys#*-2 res (square (lp 1 (##sys#integer-shift e2 -1)))))
(else
(lp (##sys#*-2 res base) (##sys#--2 e2 1)))))))
@@ -1758,7 +1855,7 @@ EOF
(bex (fx- (fx- (integer-length mant) (integer-length scl))
flonum-precision)))
(if (fx< bex 0)
- (let* ((num (arithmetic-shift mant (fxneg bex)))
+ (let* ((num (##sys#integer-shift mant (fxneg bex)))
(quo (round-quotient num scl)))
(cond ((> (integer-length quo) flonum-precision)
;; Too many bits of quotient; readjust
diff --git a/manual/Unit library b/manual/Unit library
index 05e1a588..953e25d6 100644
--- a/manual/Unit library
+++ b/manual/Unit library
@@ -43,6 +43,25 @@ set, or {{#f}} otherwise. The rightmost/least-significant bit is bit 0.
Returns the number of bits needed to represent the exact integer N in
2's complement notation.
+==== exact-integer-sqrt
+
+<procedure>(exact-integer-sqrt K)</procedure>
+
+Returns two values {{s}} and {{r}}, where {{s^2 + r = K}} and {{K < (s+1)^2}}.
+In other words, {{s}} is the closest square root we can find that's equal to or
+smaller than {{K}}, and {{r}} is the rest if {{K}} isn't a neat square of two numbers.
+
+This procedure is compatible with the R7RS specification.
+
+==== exact-integer-nth-root
+
+<procedure>(exact-integer-nth-root K N)</procedure>
+
+Like {{exact-integer-sqrt}}, but with any base value. Calculates
+{{\sqrt[N]{K}}}, the {{N}}th root of {{K}} and returns two values
+{{s}} and {{r}} where {{s^N + r = K}} and {{K < (s+1)^N}}.
+
+
==== bignum?
<procedure>(bignum? X)</procedure>
@@ -277,10 +296,12 @@ finite if both the real and imaginary components are finite.
<procedure>(signum N)</procedure>
-Returns {{1}} if {{N}} is positive, {{-1}} if {{N}}
-is negative or {{0}} if {{N}} is zero. {{signum}} is exactness preserving.
-
+For real numbers, returns {{1}} if {{N}} is positive, {{-1}} if {{N}}
+is negative or {{0}} if {{N}} is zero. {{signum}} is exactness
+preserving.
+For complex numbers, returns a complex number of the same angle but
+with magnitude 1.
=== File Input/Output
diff --git a/runtime.c b/runtime.c
index ddad4e91..57d1e014 100644
--- a/runtime.c
+++ b/runtime.c
@@ -832,7 +832,7 @@ static C_PTABLE_ENTRY *create_initial_ptable()
{
/* IMPORTANT: hardcoded table size -
this must match the number of C_pte calls + 1 (NULL terminator)! */
- C_PTABLE_ENTRY *pt = (C_PTABLE_ENTRY *)C_malloc(sizeof(C_PTABLE_ENTRY) * 77);
+ C_PTABLE_ENTRY *pt = (C_PTABLE_ENTRY *)C_malloc(sizeof(C_PTABLE_ENTRY) * 78);
int i = 0;
if(pt == NULL)
@@ -900,6 +900,7 @@ static C_PTABLE_ENTRY *create_initial_ptable()
C_pte(C_integer_to_string);
C_pte(C_flonum_to_string);
/* IMPORTANT: have you read the comments at the start and the end of this function? */
+ C_pte(C_signum);
C_pte(C_abs);
C_pte(C_u_integer_abs);
C_pte(C_negate);
@@ -5493,6 +5494,25 @@ void C_ccall C_u_integer_abs(C_word c, C_word self, C_word k, C_word x)
}
}
+void C_ccall C_signum(C_word c, C_word self, C_word k, C_word x)
+{
+ if (c != 3) {
+ C_bad_argc_2(c, 3, self);
+ } else if (x & C_FIXNUM_BIT) {
+ C_kontinue(k, C_i_fixnum_signum(x));
+ } else if (C_immediatep(x)) {
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "signum", x);
+ } else if (C_block_header(x) == C_FLONUM_TAG) {
+ C_word *a = C_alloc(C_SIZEOF_FLONUM);
+ C_kontinue(k, C_a_u_i_flonum_signum(&a, 1, x));
+ } else if (C_header_bits(x) == C_BIGNUM_TYPE) {
+ C_kontinue(k, C_bignum_negativep(x) ? C_fix(-1) : C_fix(1));
+ } else {
+ try_extended_number("\003sysextended-signum", 2, k, x);
+ }
+}
+
+
/* XXX TODO OBSOLETE: This can be removed after recompiling c-platform.scm */
C_regparm C_word C_fcall C_a_i_abs(C_word **a, int c, C_word x)
{
diff --git a/types.db b/types.db
index e5819d7e..1a4d949a 100644
--- a/types.db
+++ b/types.db
@@ -849,6 +849,12 @@
(arithmetic-shift (#(procedure #:clean #:enforce #:foldable) arithmetic-shift (integer fixnum) integer)
((integer fixnum) (##sys#integer-shift #(1) #(2))))
+(exact-integer-nth-root (#(procedure #:clean #:enforce #:foldable) exact-integer-nth-root (integer integer) integer integer)
+ ((integer integer) (##sys#exact-integer-nth-root/loc 'exact-integer-nth-root #(1) #(2))))
+
+(exact-integer-sqrt (#(procedure #:clean #:enforce #:foldable) exact-integer-sqrt (integer) integer integer)
+ ((integer) (##sys#exact-integer-sqrt #(1))))
+
(bignum? (#(procedure #:pure #:predicate bignum) bignum? (*) boolean))
(bit-set? (#(procedure #:clean #:enforce #:foldable) bit-set? (integer integer) boolean)
@@ -1230,7 +1236,16 @@
(setter (#(procedure #:clean #:enforce) setter (procedure) procedure))
(signal (procedure signal (*) . *))
-(signum (#(procedure #:clean #:enforce) signum (number) number))
+
+(signum (#(procedure #:clean #:enforce) signum (number) (or fixnum float cplxnum))
+ ((fixnum) (fixnum) (##core#inline "C_i_fixnum_signum" #(1)))
+ ((integer) (fixnum) (##core#inline "C_u_i_integer_signum" #(1)))
+ ((float) (float)
+ (##core#inline_allocate ("C_a_u_i_flonum_signum" 4) #(1)))
+ ((ratnum) (fixnum)
+ (##core#inline "C_u_i_integer_signum" (##sys#slot #(1) '1)))
+ ((cplxnum) ((or float cplxnum)) (##sys#extended-signum #(1))))
+
(software-type (#(procedure #:pure) software-type () symbol))
(software-version (#(procedure #:pure) software-version () symbol))
(string->blob (#(procedure #:clean #:enforce) string->blob (string) blob))
@@ -1262,7 +1277,6 @@
(with-exception-handler
(#(procedure #:enforce) with-exception-handler ((procedure (*) . *) (procedure () . *)) . *))
-
;; chicken (internal)
(##sys#foreign-char-argument (#(procedure #:clean #:enforce) ##sys#foreign-char-argument (char) char)
Trap