Changeset 30553 in project
 Timestamp:
 03/10/14 21:30:36 (7 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

release/4/numbers/trunk/numbers.scm
r30545 r30553 94 94 (define (badinteger loc x) (##sys#signalhook #:typeerror loc "bad argument type  not an integer" x)) 95 95 (define (badnatural loc x) (##sys#signalhook #:typeerror loc "bad argument type  must be an nonnegative integer" x)) 96 (define (badpositive loc x) (##sys#signalhook #:typeerror loc "bad argument type  must be a positive (nonzero) integer" x)) 96 97 (define (badcomplex/o loc x) (##sys#signalhook #:typeerror loc "bad argument type  complex number has no ordering" x)) 97 98 (define (badbase loc x) (##sys#signalhook #:typeerror loc "bad argument type  not a valid base" x)) … … 1571 1572 ;; from Wikipedia ;) https://en.wikipedia.org/wiki/Nth_root_algorithm 1572 1573 (define (%exactintegernthroot loc k n) 1573 (if (or (eq? 0 k) (eq? 1 k) (eq? 1 n)) ; Maybe call exactintegersqrt on n=2? 1574 (values k 0) 1575 (let ((len (integerlength k))) 1576 (if (fx< len n) ; Idea from Gambit: 2^{len1} <= k < 2^{len} 1577 (values 1 ( k 1)) ; Since we know x >= 2, we know x^{n} can't exist 1578 ;; Set initial guess to (at least) 2^ceil(ceil(log2(k))/n) 1579 (let* ((shiftamount (ceiling (/ (fx+ len 1) n))) 1580 (g0 (arithmeticshift 1 shiftamount)) 1581 (n1 (% n 1))) 1582 (let lp ((g0 g0) 1583 (g1 (%quotient loc (%+ (%* n1 g0) (%quotient loc k (%integerpower g0 n1))) n))) 1584 (if (%< g1 g0 loc) 1585 (lp g1 (%quotient loc (%+ (%* n1 g1) (%quotient loc k (%integerpower g1 n1))) n)) 1586 (values g0 (% k (%integerpower g0 n)))))))))) 1574 (cond 1575 ((or (eq? 0 k) (eq? 1 k) (eq? 1 n)) ; Maybe call exactintegersqrt on n=2? 1576 (values k 0)) 1577 ((< n 1) 1578 (badpositive loc n)) 1579 (else 1580 (let ((len (integerlength k))) 1581 (if (< len n) ; Idea from Gambit: 2^{len1} <= k < 2^{len} 1582 (values 1 ( k 1)) ; Since we know x >= 2, we know x^{n} can't exist 1583 ;; Set initial guess to (at least) 2^ceil(ceil(log2(k))/n) 1584 (let* ((shiftamount (ceiling (/ (fx+ len 1) n))) 1585 (g0 (arithmeticshift 1 shiftamount)) 1586 (n1 (% n 1))) 1587 (let lp ((g0 g0) 1588 (g1 (%quotient loc (%+ (%* n1 g0) (%quotient loc k (%integerpower g0 n1))) n))) 1589 (if (%< g1 g0 loc) 1590 (lp g1 (%quotient loc (%+ (%* n1 g1) (%quotient loc k (%integerpower g1 n1))) n)) 1591 (values g0 (% k (%integerpower g0 n))))))))))) 1587 1592 1588 1593 (define (exactintegernthroot k n)
Note: See TracChangeset
for help on using the changeset viewer.