82 |
4877 4889 4903 4909 4919 4931 4933 4937 4943 4951 4957 4967 4969 4973 |
4877 4889 4903 4909 4919 4931 4933 4937 4943 4951 4957 4967 4969 4973 |
83 |
4987 4993 4999 5003]) |
4987 4993 4999 5003]) |
84 |
|
|
85 |
|
;; The variable math-prime-factors-finished is set by calcFunc-prfac to |
86 |
|
;; indicate whether factoring is complete, and used by calcFunc-factors, |
87 |
|
;; calcFunc-totient and calcFunc-moebius. |
88 |
|
(defvar math-prime-factors-finished) |
89 |
|
|
90 |
;;; Combinatorics |
;;; Combinatorics |
91 |
|
|
92 |
(defun calc-gcd (arg) |
(defun calc-gcd (arg) |
200 |
(res (math-prime-test n iters))) |
(res (math-prime-test n iters))) |
201 |
(calc-report-prime-test res)))) |
(calc-report-prime-test res)))) |
202 |
|
|
203 |
|
(defvar calc-verbose-nextprime nil) |
204 |
|
|
205 |
(defun calc-next-prime (iters) |
(defun calc-next-prime (iters) |
206 |
(interactive "p") |
(interactive "p") |
207 |
(calc-slow-wrapper |
(calc-slow-wrapper |
393 |
(if (math-evenp temp) |
(if (math-evenp temp) |
394 |
even |
even |
395 |
(math-div (calcFunc-fact n) even)))) |
(math-div (calcFunc-fact n) even)))) |
396 |
(list 'calcFunc-dfact max)))) |
(list 'calcFunc-dfact n)))) |
397 |
((equal n '(var inf var-inf)) n) |
((equal n '(var inf var-inf)) n) |
398 |
(t (calc-record-why 'natnump n) |
(t (calc-record-why 'natnump n) |
399 |
(list 'calcFunc-dfact n)))) |
(list 'calcFunc-dfact n)))) |
491 |
(math-stirling-number n m 0)) |
(math-stirling-number n m 0)) |
492 |
|
|
493 |
(defvar math-stirling-cache (vector [[1]] [[1]])) |
(defvar math-stirling-cache (vector [[1]] [[1]])) |
494 |
|
|
495 |
|
;; The variable math-stirling-local-cache is local to |
496 |
|
;; math-stirling-number, but is used by math-stirling-1 |
497 |
|
;; and math-stirling-2, which are called by math-stirling-number. |
498 |
|
(defvar math-stirling-local-cache) |
499 |
|
|
500 |
(defun math-stirling-number (n m k) |
(defun math-stirling-number (n m k) |
501 |
(or (math-num-natnump n) (math-reject-arg n 'natnump)) |
(or (math-num-natnump n) (math-reject-arg n 'natnump)) |
502 |
(or (math-num-natnump m) (math-reject-arg m 'natnump)) |
(or (math-num-natnump m) (math-reject-arg m 'natnump)) |
506 |
(or (integerp m) (math-reject-arg m 'fixnump)) |
(or (integerp m) (math-reject-arg m 'fixnump)) |
507 |
(if (< n m) |
(if (< n m) |
508 |
0 |
0 |
509 |
(let ((cache (aref math-stirling-cache k))) |
(let ((math-stirling-local-cache (aref math-stirling-cache k))) |
510 |
(while (<= (length cache) n) |
(while (<= (length math-stirling-local-cache) n) |
511 |
(let ((i (1- (length cache))) |
(let ((i (1- (length math-stirling-local-cache))) |
512 |
row) |
row) |
513 |
(setq cache (vconcat cache (make-vector (length cache) nil))) |
(setq math-stirling-local-cache |
514 |
(aset math-stirling-cache k cache) |
(vconcat math-stirling-local-cache |
515 |
(while (< (setq i (1+ i)) (length cache)) |
(make-vector (length math-stirling-local-cache) nil))) |
516 |
(aset cache i (setq row (make-vector (1+ i) nil))) |
(aset math-stirling-cache k math-stirling-local-cache) |
517 |
|
(while (< (setq i (1+ i)) (length math-stirling-local-cache)) |
518 |
|
(aset math-stirling-local-cache i (setq row (make-vector (1+ i) nil))) |
519 |
(aset row 0 0) |
(aset row 0 0) |
520 |
(aset row i 1)))) |
(aset row i 1)))) |
521 |
(if (= k 1) |
(if (= k 1) |
523 |
(math-stirling-2 n m))))) |
(math-stirling-2 n m))))) |
524 |
|
|
525 |
(defun math-stirling-1 (n m) |
(defun math-stirling-1 (n m) |
526 |
(or (aref (aref cache n) m) |
(or (aref (aref math-stirling-local-cache n) m) |
527 |
(aset (aref cache n) m |
(aset (aref math-stirling-local-cache n) m |
528 |
(math-add (math-stirling-1 (1- n) (1- m)) |
(math-add (math-stirling-1 (1- n) (1- m)) |
529 |
(math-mul (- 1 n) (math-stirling-1 (1- n) m)))))) |
(math-mul (- 1 n) (math-stirling-1 (1- n) m)))))) |
530 |
|
|
531 |
(defun math-stirling-2 (n m) |
(defun math-stirling-2 (n m) |
532 |
(or (aref (aref cache n) m) |
(or (aref (aref math-stirling-local-cache n) m) |
533 |
(aset (aref cache n) m |
(aset (aref math-stirling-local-cache n) m |
534 |
(math-add (math-stirling-2 (1- n) (1- m)) |
(math-add (math-stirling-2 (1- n) (1- m)) |
535 |
(math-mul m (math-stirling-2 (1- n) m)))))) |
(math-mul m (math-stirling-2 (1- n) m)))))) |
536 |
|
|
542 |
|
|
543 |
;;; Produce a random 10-bit integer, with (random) if no seed provided, |
;;; Produce a random 10-bit integer, with (random) if no seed provided, |
544 |
;;; or else with Numerical Recipes algorithm ran3 / Knuth 3.2.2-A. |
;;; or else with Numerical Recipes algorithm ran3 / Knuth 3.2.2-A. |
545 |
|
|
546 |
|
(defvar var-RandSeed nil) |
547 |
|
(defvar math-random-cache nil) |
548 |
|
(defvar math-gaussian-cache nil) |
549 |
|
|
550 |
(defun math-init-random-base () |
(defun math-init-random-base () |
551 |
(if (and (boundp 'var-RandSeed) var-RandSeed) |
(if var-RandSeed |
552 |
(if (eq (car-safe var-RandSeed) 'vec) |
(if (eq (car-safe var-RandSeed) 'vec) |
553 |
nil |
nil |
554 |
(if (Math-integerp var-RandSeed) |
(if (Math-integerp var-RandSeed) |
575 |
(random t) |
(random t) |
576 |
(setq var-RandSeed nil |
(setq var-RandSeed nil |
577 |
math-random-cache nil |
math-random-cache nil |
|
i 0 |
|
578 |
math-random-shift -4) ; assume RAND_MAX >= 16383 |
math-random-shift -4) ; assume RAND_MAX >= 16383 |
579 |
;; This exercises the random number generator and also helps |
;; This exercises the random number generator and also helps |
580 |
;; deduce a better value for RAND_MAX. |
;; deduce a better value for RAND_MAX. |
581 |
(while (< (setq i (1+ i)) 30) |
(let ((i 0)) |
582 |
(if (> (lsh (math-abs (random)) math-random-shift) 4095) |
(while (< (setq i (1+ i)) 30) |
583 |
(setq math-random-shift (1- math-random-shift))))) |
(if (> (lsh (math-abs (random)) math-random-shift) 4095) |
584 |
|
(setq math-random-shift (1- math-random-shift)))))) |
585 |
(setq math-last-RandSeed var-RandSeed |
(setq math-last-RandSeed var-RandSeed |
586 |
math-gaussian-cache nil)) |
math-gaussian-cache nil)) |
587 |
|
|
603 |
;;; Avoid various pitfalls that may lurk in the built-in (random) function! |
;;; Avoid various pitfalls that may lurk in the built-in (random) function! |
604 |
;;; Shuffling algorithm from Numerical Recipes, section 7.1. |
;;; Shuffling algorithm from Numerical Recipes, section 7.1. |
605 |
(defun math-random-digit () |
(defun math-random-digit () |
606 |
(let (i) |
(let (i math-random-last) |
607 |
(or (and (boundp 'var-RandSeed) (eq var-RandSeed math-last-RandSeed)) |
(or (eq var-RandSeed math-last-RandSeed) |
608 |
(math-init-random-base)) |
(math-init-random-base)) |
609 |
(or math-random-cache |
(or math-random-cache |
610 |
(progn |
(progn |
619 |
(aset math-random-cache i (math-random-base)) |
(aset math-random-cache i (math-random-base)) |
620 |
(>= math-random-last 1000))) |
(>= math-random-last 1000))) |
621 |
math-random-last)) |
math-random-last)) |
|
(setq math-random-cache nil) |
|
622 |
|
|
623 |
;;; Produce an N-digit random integer. |
;;; Produce an N-digit random integer. |
624 |
(defun math-random-digits (n) |
(defun math-random-digits (n) |
658 |
(setq math-gaussian-cache (cons calc-internal-prec |
(setq math-gaussian-cache (cons calc-internal-prec |
659 |
(math-mul v1 fac))) |
(math-mul v1 fac))) |
660 |
(math-mul v2 fac)))))) |
(math-mul v2 fac)))))) |
|
(setq math-gaussian-cache nil) |
|
661 |
|
|
662 |
;;; Produce a random integer or real 0 <= N < MAX. |
;;; Produce a random integer or real 0 <= N < MAX. |
663 |
(defun calcFunc-random (max) |
(defun calcFunc-random (max) |
783 |
;;; (nil unknown) if non-prime with no known factors, |
;;; (nil unknown) if non-prime with no known factors, |
784 |
;;; (t) if prime, |
;;; (t) if prime, |
785 |
;;; (maybe N P) if probably prime (after N iters with probability P%) |
;;; (maybe N P) if probably prime (after N iters with probability P%) |
786 |
|
(defvar math-prime-test-cache '(-1)) |
787 |
|
|
788 |
|
(defvar math-prime-test-cache-k) |
789 |
|
(defvar math-prime-test-cache-q) |
790 |
|
(defvar math-prime-test-cache-nm1) |
791 |
|
|
792 |
(defun math-prime-test (n iters) |
(defun math-prime-test (n iters) |
793 |
(if (and (Math-vectorp n) (cdr n)) |
(if (and (Math-vectorp n) (cdr n)) |
794 |
(setq n (nth (1- (length n)) n))) |
(setq n (nth (1- (length n)) n))) |
873 |
(1- iters) |
(1- iters) |
874 |
0))) |
0))) |
875 |
res)) |
res)) |
|
(defvar math-prime-test-cache '(-1)) |
|
876 |
|
|
877 |
(defun calcFunc-prime (n &optional iters) |
(defun calcFunc-prime (n &optional iters) |
878 |
(or (math-num-integerp n) (math-reject-arg n 'integerp)) |
(or (math-num-integerp n) (math-reject-arg n 'integerp)) |
988 |
(if (Math-realp n) |
(if (Math-realp n) |
989 |
(calcFunc-nextprime (math-trunc n) iters) |
(calcFunc-nextprime (math-trunc n) iters) |
990 |
(math-reject-arg n 'integerp)))) |
(math-reject-arg n 'integerp)))) |
|
(setq calc-verbose-nextprime nil) |
|
991 |
|
|
992 |
(defun calcFunc-prevprime (n &optional iters) |
(defun calcFunc-prevprime (n &optional iters) |
993 |
(if (Math-integerp n) |
(if (Math-integerp n) |