diff options
Diffstat (limited to 'lisp/calc/calc-math.el')
| -rw-r--r-- | lisp/calc/calc-math.el | 97 |
1 files changed, 19 insertions, 78 deletions
diff --git a/lisp/calc/calc-math.el b/lisp/calc/calc-math.el index 50c8758ace2..4ca8515989b 100644 --- a/lisp/calc/calc-math.el +++ b/lisp/calc/calc-math.el @@ -25,6 +25,8 @@ ;; This file is autoloaded from calc-ext.el. + +(require 'cl-lib) (require 'calc-ext) (require 'calc-macs) @@ -95,8 +97,7 @@ If this can't be done, return NIL." (and (<= calc-internal-prec math-emacs-precision) (math-realp x) - (let* ((fx (math-float x)) - (xpon (+ (nth 2 x) (1- (math-numdigs (nth 1 x)))))) + (let* ((xpon (+ (nth 2 x) (1- (math-numdigs (nth 1 x)))))) (and (<= math-smallest-emacs-expt xpon) (<= xpon math-largest-emacs-expt) (condition-case nil @@ -371,51 +372,15 @@ If this can't be done, return NIL." ;;; with an overestimate always works, even using truncating integer division! (defun math-isqrt (a) (cond ((Math-zerop a) a) - ((not (math-natnump a)) + ((not (natnump a)) (math-reject-arg a 'natnump)) - ((integerp a) - (math-isqrt-small a)) - (t - (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a)))))))) + (t (cl-isqrt a)))) (defun calcFunc-isqrt (a) (if (math-realp a) (math-isqrt (math-floor a)) (math-floor (math-sqrt a)))) - -;;; This returns (flag . result) where the flag is t if A is a perfect square. -(defun math-isqrt-bignum (a) ; [P.l L] - (let ((len (length a))) - (if (= (% len 2) 0) - (let* ((top (nthcdr (- len 2) a))) - (math-isqrt-bignum-iter - a - (math-scale-bignum-digit-size - (math-bignum-big - (1+ (math-isqrt-small - (+ (* (nth 1 top) math-bignum-digit-size) (car top))))) - (1- (/ len 2))))) - (let* ((top (nth (1- len) a))) - (math-isqrt-bignum-iter - a - (math-scale-bignum-digit-size - (list (1+ (math-isqrt-small top))) - (/ len 2))))))) - -(defun math-isqrt-bignum-iter (a guess) ; [l L l] - (math-working "isqrt" (cons 'bigpos guess)) - (let* ((q (math-div-bignum a guess)) - (s (math-add-bignum (car q) guess)) - (g2 (math-div2-bignum s)) - (comp (math-compare-bignum g2 guess))) - (if (< comp 0) - (math-isqrt-bignum-iter a g2) - (cons (and (= comp 0) - (math-zerop-bignum (cdr q)) - (= (% (car s) 2) 0)) - guess)))) - (defun math-zerop-bignum (a) (and (eq (car a) 0) (progn @@ -428,19 +393,6 @@ If this can't be done, return NIL." n (1- n))) a) -(defun math-isqrt-small (a) ; A > 0. [S S] - (let ((g (cond ((>= a 1000000) 10000) - ((>= a 10000) 1000) - ((>= a 100) 100) - (t 10))) - g2) - (while (< (setq g2 (/ (+ g (/ a g)) 2)) g) - (setq g g2)) - g)) - - - - ;;; Compute the square root of a number. ;;; [T N] if possible, else [F N] if possible, else [C N]. [Public] (defun math-sqrt (a) @@ -449,32 +401,24 @@ If this can't be done, return NIL." (and (math-known-nonposp a) (math-imaginary (math-sqrt (math-neg a)))) (and (integerp a) - (let ((sqrt (math-isqrt-small a))) + (let ((sqrt (cl-isqrt a))) (if (= (* sqrt sqrt) a) sqrt (if calc-symbolic-mode (list 'calcFunc-sqrt a) (math-sqrt-float (math-float a) (math-float sqrt)))))) - (and (eq (car-safe a) 'bigpos) - (let* ((res (math-isqrt-bignum (cdr a))) - (sqrt (math-normalize (cons 'bigpos (cdr res))))) - (if (car res) - sqrt - (if calc-symbolic-mode - (list 'calcFunc-sqrt a) - (math-sqrt-float (math-float a) (math-float sqrt)))))) (and (eq (car-safe a) 'frac) - (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a))))) - (num-sqrt (math-normalize (cons 'bigpos (cdr num-res)))) - (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a))))) - (den-sqrt (math-normalize (cons 'bigpos (cdr den-res))))) - (if (and (car num-res) (car den-res)) + (let* ((num-sqrt (cl-isqrt (nth 1 a))) + (num-exact (= (* num-sqrt num-sqrt) (nth 1 a))) + (den-sqrt (cl-isqrt (nth 2 a))) + (den-exact (= (* den-sqrt den-sqrt) (nth 2 a)))) + (if (and num-exact den-exact) (list 'frac num-sqrt den-sqrt) (if calc-symbolic-mode - (if (or (car num-res) (car den-res)) - (math-div (if (car num-res) + (if (or num-exact den-exact) + (math-div (if num-exact num-sqrt (list 'calcFunc-sqrt (nth 1 a))) - (if (car den-res) + (if den-exact den-sqrt (list 'calcFunc-sqrt (nth 2 a)))) (list 'calcFunc-sqrt a)) (math-sqrt-float (math-float a) @@ -482,12 +426,9 @@ If this can't be done, return NIL." (and (eq (car-safe a) 'float) (if calc-symbolic-mode (if (= (% (nth 2 a) 2) 0) - (let ((res (math-isqrt-bignum - (cdr (Math-bignum-test (nth 1 a)))))) - (if (car res) - (math-make-float (math-normalize - (cons 'bigpos (cdr res))) - (/ (nth 2 a) 2)) + (let ((res (cl-isqrt (nth 1 a)))) + (if (= (* res res) (nth 1 a)) + (math-make-float res (/ (nth 2 a) 2)) (signal 'inexact-result nil))) (signal 'inexact-result nil)) (math-sqrt-float a))) @@ -551,7 +492,7 @@ If this can't be done, return NIL." (if (null guess) (let ((ldiff (- (math-numdigs (nth 1 a)) 6))) (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff))) - (setq guess (math-make-float (math-isqrt-small + (setq guess (math-make-float (cl-isqrt (math-scale-int (nth 1 a) (- ldiff))) (/ (+ (nth 2 a) ldiff) 2))))) (math-sqrt-float-iter a guess))))) @@ -1697,7 +1638,7 @@ If this can't be done, return NIL." (while (not (Math-lessp x pow)) (setq pows (cons pow pows) pow (math-sqr pow))) - (setq n (lsh 1 (1- (length pows))) + (setq n (ash 1 (1- (length pows))) sum n pow (car pows)) (while (and (setq pows (cdr pows)) |