41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
1 ;;; calc-comb.el --- combinatoric functions for Calc
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
2
|
64325
|
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
|
79702
|
4 ;; 2005, 2006, 2007, 2008 Free Software Foundation, Inc.
|
41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
5
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
6 ;; Author: David Gillespie <daveg@synaptics.com>
|
77465
|
7 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
|
40785
|
8
|
|
9 ;; This file is part of GNU Emacs.
|
|
10
|
94654
|
11 ;; GNU Emacs is free software: you can redistribute it and/or modify
|
76595
|
12 ;; it under the terms of the GNU General Public License as published by
|
94654
|
13 ;; the Free Software Foundation, either version 3 of the License, or
|
|
14 ;; (at your option) any later version.
|
40785
|
15
|
76595
|
16 ;; GNU Emacs is distributed in the hope that it will be useful,
|
|
17 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
18 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
19 ;; GNU General Public License for more details.
|
|
20
|
|
21 ;; You should have received a copy of the GNU General Public License
|
94654
|
22 ;; along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
|
40785
|
23
|
41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
24 ;;; Commentary:
|
40785
|
25
|
41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
26 ;;; Code:
|
40785
|
27
|
|
28 ;; This file is autoloaded from calc-ext.el.
|
58649
|
29
|
40785
|
30 (require 'calc-ext)
|
|
31 (require 'calc-macs)
|
|
32
|
41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
33 (defconst math-primes-table
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
34 [2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
35 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
36 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
37 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
38 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
39 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
40 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
41 719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
42 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
43 953 967 971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
44 1051 1061 1063 1069 1087 1091 1093 1097 1103 1109 1117 1123 1129 1151
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
45 1153 1163 1171 1181 1187 1193 1201 1213 1217 1223 1229 1231 1237 1249
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
46 1259 1277 1279 1283 1289 1291 1297 1301 1303 1307 1319 1321 1327 1361
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
47 1367 1373 1381 1399 1409 1423 1427 1429 1433 1439 1447 1451 1453 1459
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
48 1471 1481 1483 1487 1489 1493 1499 1511 1523 1531 1543 1549 1553 1559
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
49 1567 1571 1579 1583 1597 1601 1607 1609 1613 1619 1621 1627 1637 1657
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
50 1663 1667 1669 1693 1697 1699 1709 1721 1723 1733 1741 1747 1753 1759
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
51 1777 1783 1787 1789 1801 1811 1823 1831 1847 1861 1867 1871 1873 1877
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
52 1879 1889 1901 1907 1913 1931 1933 1949 1951 1973 1979 1987 1993 1997
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
53 1999 2003 2011 2017 2027 2029 2039 2053 2063 2069 2081 2083 2087 2089
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
54 2099 2111 2113 2129 2131 2137 2141 2143 2153 2161 2179 2203 2207 2213
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
55 2221 2237 2239 2243 2251 2267 2269 2273 2281 2287 2293 2297 2309 2311
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
56 2333 2339 2341 2347 2351 2357 2371 2377 2381 2383 2389 2393 2399 2411
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
57 2417 2423 2437 2441 2447 2459 2467 2473 2477 2503 2521 2531 2539 2543
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
58 2549 2551 2557 2579 2591 2593 2609 2617 2621 2633 2647 2657 2659 2663
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
59 2671 2677 2683 2687 2689 2693 2699 2707 2711 2713 2719 2729 2731 2741
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
60 2749 2753 2767 2777 2789 2791 2797 2801 2803 2819 2833 2837 2843 2851
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
61 2857 2861 2879 2887 2897 2903 2909 2917 2927 2939 2953 2957 2963 2969
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
62 2971 2999 3001 3011 3019 3023 3037 3041 3049 3061 3067 3079 3083 3089
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
63 3109 3119 3121 3137 3163 3167 3169 3181 3187 3191 3203 3209 3217 3221
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
64 3229 3251 3253 3257 3259 3271 3299 3301 3307 3313 3319 3323 3329 3331
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
65 3343 3347 3359 3361 3371 3373 3389 3391 3407 3413 3433 3449 3457 3461
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
66 3463 3467 3469 3491 3499 3511 3517 3527 3529 3533 3539 3541 3547 3557
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
67 3559 3571 3581 3583 3593 3607 3613 3617 3623 3631 3637 3643 3659 3671
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
68 3673 3677 3691 3697 3701 3709 3719 3727 3733 3739 3761 3767 3769 3779
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
69 3793 3797 3803 3821 3823 3833 3847 3851 3853 3863 3877 3881 3889 3907
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
70 3911 3917 3919 3923 3929 3931 3943 3947 3967 3989 4001 4003 4007 4013
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
71 4019 4021 4027 4049 4051 4057 4073 4079 4091 4093 4099 4111 4127 4129
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
72 4133 4139 4153 4157 4159 4177 4201 4211 4217 4219 4229 4231 4241 4243
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
73 4253 4259 4261 4271 4273 4283 4289 4297 4327 4337 4339 4349 4357 4363
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
74 4373 4391 4397 4409 4421 4423 4441 4447 4451 4457 4463 4481 4483 4493
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
75 4507 4513 4517 4519 4523 4547 4549 4561 4567 4583 4591 4597 4603 4621
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
76 4637 4639 4643 4649 4651 4657 4663 4673 4679 4691 4703 4721 4723 4729
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
77 4733 4751 4759 4783 4787 4789 4793 4799 4801 4813 4817 4831 4861 4871
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
78 4877 4889 4903 4909 4919 4931 4933 4937 4943 4951 4957 4967 4969 4973
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
79 4987 4993 4999 5003])
|
40785
|
80
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
81 ;; The variable math-prime-factors-finished is set by calcFunc-prfac to
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
82 ;; indicate whether factoring is complete, and used by calcFunc-factors,
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
83 ;; calcFunc-totient and calcFunc-moebius.
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
84 (defvar math-prime-factors-finished)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
85
|
40785
|
86 ;;; Combinatorics
|
|
87
|
|
88 (defun calc-gcd (arg)
|
|
89 (interactive "P")
|
|
90 (calc-slow-wrapper
|
41047
|
91 (calc-binary-op "gcd" 'calcFunc-gcd arg)))
|
40785
|
92
|
|
93 (defun calc-lcm (arg)
|
|
94 (interactive "P")
|
|
95 (calc-slow-wrapper
|
41047
|
96 (calc-binary-op "lcm" 'calcFunc-lcm arg)))
|
40785
|
97
|
|
98 (defun calc-extended-gcd ()
|
|
99 (interactive)
|
|
100 (calc-slow-wrapper
|
41047
|
101 (calc-enter-result 2 "egcd" (cons 'calcFunc-egcd (calc-top-list-n 2)))))
|
40785
|
102
|
|
103 (defun calc-factorial (arg)
|
|
104 (interactive "P")
|
|
105 (calc-slow-wrapper
|
41047
|
106 (calc-unary-op "fact" 'calcFunc-fact arg)))
|
40785
|
107
|
|
108 (defun calc-gamma (arg)
|
|
109 (interactive "P")
|
|
110 (calc-slow-wrapper
|
41047
|
111 (calc-unary-op "gmma" 'calcFunc-gamma arg)))
|
40785
|
112
|
|
113 (defun calc-double-factorial (arg)
|
|
114 (interactive "P")
|
|
115 (calc-slow-wrapper
|
41047
|
116 (calc-unary-op "dfac" 'calcFunc-dfact arg)))
|
40785
|
117
|
|
118 (defun calc-choose (arg)
|
|
119 (interactive "P")
|
|
120 (calc-slow-wrapper
|
|
121 (if (calc-is-hyperbolic)
|
|
122 (calc-binary-op "perm" 'calcFunc-perm arg)
|
41047
|
123 (calc-binary-op "chos" 'calcFunc-choose arg))))
|
40785
|
124
|
|
125 (defun calc-perm (arg)
|
|
126 (interactive "P")
|
|
127 (calc-hyperbolic-func)
|
41047
|
128 (calc-choose arg))
|
40785
|
129
|
|
130 (defvar calc-last-random-limit '(float 1 0))
|
|
131 (defun calc-random (n)
|
|
132 (interactive "P")
|
|
133 (calc-slow-wrapper
|
|
134 (if n
|
|
135 (calc-enter-result 0 "rand" (list 'calcFunc-random
|
|
136 (calc-get-random-limit
|
|
137 (prefix-numeric-value n))))
|
|
138 (calc-enter-result 1 "rand" (list 'calcFunc-random
|
|
139 (calc-get-random-limit
|
41047
|
140 (calc-top-n 1)))))))
|
40785
|
141
|
|
142 (defun calc-get-random-limit (val)
|
|
143 (if (eq val 0)
|
|
144 calc-last-random-limit
|
41047
|
145 (setq calc-last-random-limit val)))
|
40785
|
146
|
|
147 (defun calc-rrandom ()
|
|
148 (interactive)
|
|
149 (calc-slow-wrapper
|
|
150 (setq calc-last-random-limit '(float 1 0))
|
41047
|
151 (calc-enter-result 0 "rand" (list 'calcFunc-random '(float 1 0)))))
|
40785
|
152
|
|
153 (defun calc-random-again (arg)
|
|
154 (interactive "p")
|
|
155 (calc-slow-wrapper
|
|
156 (while (>= (setq arg (1- arg)) 0)
|
|
157 (calc-enter-result 0 "rand" (list 'calcFunc-random
|
41047
|
158 calc-last-random-limit)))))
|
40785
|
159
|
|
160 (defun calc-shuffle (n)
|
|
161 (interactive "P")
|
|
162 (calc-slow-wrapper
|
|
163 (if n
|
|
164 (calc-enter-result 1 "shuf" (list 'calcFunc-shuffle
|
|
165 (prefix-numeric-value n)
|
|
166 (calc-get-random-limit
|
|
167 (calc-top-n 1))))
|
|
168 (calc-enter-result 2 "shuf" (list 'calcFunc-shuffle
|
|
169 (calc-top-n 1)
|
|
170 (calc-get-random-limit
|
41047
|
171 (calc-top-n 2)))))))
|
40785
|
172
|
|
173 (defun calc-report-prime-test (res)
|
|
174 (cond ((eq (car res) t)
|
|
175 (calc-record-message "prim" "Prime (guaranteed)"))
|
|
176 ((eq (car res) nil)
|
|
177 (if (cdr res)
|
|
178 (if (eq (nth 1 res) 'unknown)
|
|
179 (calc-record-message
|
|
180 "prim" "Non-prime (factors unknown)")
|
|
181 (calc-record-message
|
|
182 "prim" "Non-prime (%s is a factor)"
|
|
183 (math-format-number (nth 1 res))))
|
|
184 (calc-record-message "prim" "Non-prime")))
|
|
185 (t
|
|
186 (calc-record-message
|
|
187 "prim" "Probably prime (%d iters; %s%% chance of error)"
|
|
188 (nth 1 res)
|
|
189 (let ((calc-float-format '(fix 2)))
|
41047
|
190 (math-format-number (nth 2 res)))))))
|
40785
|
191
|
|
192 (defun calc-prime-test (iters)
|
|
193 (interactive "p")
|
|
194 (calc-slow-wrapper
|
|
195 (let* ((n (calc-top-n 1))
|
|
196 (res (math-prime-test n iters)))
|
41047
|
197 (calc-report-prime-test res))))
|
40785
|
198
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
199 (defvar calc-verbose-nextprime nil)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
200
|
40785
|
201 (defun calc-next-prime (iters)
|
|
202 (interactive "p")
|
|
203 (calc-slow-wrapper
|
|
204 (let ((calc-verbose-nextprime t))
|
|
205 (if (calc-is-inverse)
|
|
206 (calc-enter-result 1 "prvp" (list 'calcFunc-prevprime
|
|
207 (calc-top-n 1) (math-abs iters)))
|
|
208 (calc-enter-result 1 "nxtp" (list 'calcFunc-nextprime
|
41047
|
209 (calc-top-n 1) (math-abs iters)))))))
|
40785
|
210
|
|
211 (defun calc-prev-prime (iters)
|
|
212 (interactive "p")
|
|
213 (calc-invert-func)
|
41047
|
214 (calc-next-prime iters))
|
40785
|
215
|
|
216 (defun calc-prime-factors (iters)
|
|
217 (interactive "p")
|
|
218 (calc-slow-wrapper
|
|
219 (let ((res (calcFunc-prfac (calc-top-n 1))))
|
|
220 (if (not math-prime-factors-finished)
|
|
221 (calc-record-message "pfac" "Warning: May not be fully factored"))
|
41047
|
222 (calc-enter-result 1 "pfac" res))))
|
40785
|
223
|
|
224 (defun calc-totient (arg)
|
|
225 (interactive "P")
|
|
226 (calc-slow-wrapper
|
41047
|
227 (calc-unary-op "phi" 'calcFunc-totient arg)))
|
40785
|
228
|
|
229 (defun calc-moebius (arg)
|
|
230 (interactive "P")
|
|
231 (calc-slow-wrapper
|
41047
|
232 (calc-unary-op "mu" 'calcFunc-moebius arg)))
|
40785
|
233
|
|
234
|
|
235 (defun calcFunc-gcd (a b)
|
|
236 (if (Math-messy-integerp a)
|
|
237 (setq a (math-trunc a)))
|
|
238 (if (Math-messy-integerp b)
|
|
239 (setq b (math-trunc b)))
|
|
240 (cond ((and (Math-integerp a) (Math-integerp b))
|
|
241 (math-gcd a b))
|
|
242 ((Math-looks-negp a)
|
|
243 (calcFunc-gcd (math-neg a) b))
|
|
244 ((Math-looks-negp b)
|
|
245 (calcFunc-gcd a (math-neg b)))
|
|
246 ((Math-zerop a) b)
|
|
247 ((Math-zerop b) a)
|
|
248 ((and (Math-ratp a)
|
|
249 (Math-ratp b))
|
|
250 (math-make-frac (math-gcd (if (eq (car-safe a) 'frac) (nth 1 a) a)
|
|
251 (if (eq (car-safe b) 'frac) (nth 1 b) b))
|
|
252 (calcFunc-lcm
|
|
253 (if (eq (car-safe a) 'frac) (nth 2 a) 1)
|
|
254 (if (eq (car-safe b) 'frac) (nth 2 b) 1))))
|
|
255 ((not (Math-integerp a))
|
|
256 (calc-record-why 'integerp a)
|
|
257 (list 'calcFunc-gcd a b))
|
|
258 (t
|
|
259 (calc-record-why 'integerp b)
|
41047
|
260 (list 'calcFunc-gcd a b))))
|
40785
|
261
|
|
262 (defun calcFunc-lcm (a b)
|
|
263 (let ((g (calcFunc-gcd a b)))
|
|
264 (if (Math-numberp g)
|
|
265 (math-div (math-mul a b) g)
|
41047
|
266 (list 'calcFunc-lcm a b))))
|
40785
|
267
|
|
268 (defun calcFunc-egcd (a b) ; Knuth section 4.5.2
|
|
269 (cond
|
|
270 ((not (Math-integerp a))
|
|
271 (if (Math-messy-integerp a)
|
|
272 (calcFunc-egcd (math-trunc a) b)
|
|
273 (calc-record-why 'integerp a)
|
|
274 (list 'calcFunc-egcd a b)))
|
|
275 ((not (Math-integerp b))
|
|
276 (if (Math-messy-integerp b)
|
|
277 (calcFunc-egcd a (math-trunc b))
|
|
278 (calc-record-why 'integerp b)
|
|
279 (list 'calcFunc-egcd a b)))
|
|
280 (t
|
|
281 (let ((u1 1) (u2 0) (u3 a)
|
|
282 (v1 0) (v2 1) (v3 b)
|
|
283 t1 t2 q)
|
|
284 (while (not (eq v3 0))
|
|
285 (setq q (math-idivmod u3 v3)
|
|
286 t1 (math-sub u1 (math-mul v1 (car q)))
|
|
287 t2 (math-sub u2 (math-mul v2 (car q)))
|
|
288 u1 v1 u2 v2 u3 v3
|
|
289 v1 t1 v2 t2 v3 (cdr q)))
|
41047
|
290 (list 'vec u3 u1 u2)))))
|
40785
|
291
|
|
292
|
|
293 ;;; Factorial and related functions.
|
|
294
|
81537
|
295 (defconst math-small-factorial-table
|
82440
ddd464fd7cb4
(math-small-factorial-table): Remove `eval-when-compile's.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
296 (vector 1 1 2 6 24 120 720 5040 40320 362880 3628800 39916800
|
ddd464fd7cb4
(math-small-factorial-table): Remove `eval-when-compile's.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
297 (math-read-number-simple "479001600")
|
ddd464fd7cb4
(math-small-factorial-table): Remove `eval-when-compile's.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
298 (math-read-number-simple "6227020800")
|
ddd464fd7cb4
(math-small-factorial-table): Remove `eval-when-compile's.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
299 (math-read-number-simple "87178291200")
|
ddd464fd7cb4
(math-small-factorial-table): Remove `eval-when-compile's.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
300 (math-read-number-simple "1307674368000")
|
ddd464fd7cb4
(math-small-factorial-table): Remove `eval-when-compile's.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
301 (math-read-number-simple "20922789888000")
|
ddd464fd7cb4
(math-small-factorial-table): Remove `eval-when-compile's.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
302 (math-read-number-simple "355687428096000")
|
ddd464fd7cb4
(math-small-factorial-table): Remove `eval-when-compile's.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
303 (math-read-number-simple "6402373705728000")
|
ddd464fd7cb4
(math-small-factorial-table): Remove `eval-when-compile's.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
304 (math-read-number-simple "121645100408832000")
|
ddd464fd7cb4
(math-small-factorial-table): Remove `eval-when-compile's.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
305 (math-read-number-simple "2432902008176640000")))
|
81537
|
306
|
40785
|
307 (defun calcFunc-fact (n) ; [I I] [F F] [Public]
|
|
308 (let (temp)
|
|
309 (cond ((Math-integer-negp n)
|
|
310 (if calc-infinite-mode
|
|
311 '(var uinf var-uinf)
|
|
312 (math-reject-arg n 'range)))
|
|
313 ((integerp n)
|
|
314 (if (<= n 20)
|
81537
|
315 (aref math-small-factorial-table n)
|
40785
|
316 (math-factorial-iter (1- n) 2 1)))
|
|
317 ((and (math-messy-integerp n)
|
|
318 (Math-lessp n 100))
|
|
319 (math-inexact-result)
|
|
320 (setq temp (math-trunc n))
|
|
321 (if (>= temp 0)
|
|
322 (if (<= temp 20)
|
|
323 (math-float (calcFunc-fact temp))
|
|
324 (math-with-extra-prec 1
|
|
325 (math-factorial-iter (1- temp) 2 '(float 1 0))))
|
|
326 (math-reject-arg n 'range)))
|
|
327 ((math-numberp n)
|
|
328 (let* ((q (math-quarter-integer n))
|
|
329 (tn (and q (Math-lessp n 1000) (Math-lessp -1000 n)
|
|
330 (1+ (math-floor n)))))
|
|
331 (cond ((and tn (= q 2)
|
|
332 (or calc-symbolic-mode (< (math-abs tn) 20)))
|
|
333 (let ((q (if (< tn 0)
|
|
334 (math-div
|
|
335 (math-pow -2 (- tn))
|
|
336 (math-double-factorial-iter (* -2 tn) 3 1 2))
|
49598
|
337 (math-div
|
40785
|
338 (math-double-factorial-iter (* 2 tn) 3 1 2)
|
|
339 (math-pow 2 tn)))))
|
|
340 (math-mul q (if calc-symbolic-mode
|
|
341 (list 'calcFunc-sqrt '(var pi var-pi))
|
|
342 (math-sqrt-pi)))))
|
|
343 ((and tn (>= tn 0) (< tn 20)
|
|
344 (memq q '(1 3)))
|
|
345 (math-inexact-result)
|
|
346 (math-div
|
|
347 (math-mul (math-double-factorial-iter (* 4 tn) q 1 4)
|
|
348 (if (= q 1) (math-gamma-1q) (math-gamma-3q)))
|
|
349 (math-pow 4 tn)))
|
|
350 (t
|
|
351 (math-inexact-result)
|
|
352 (math-with-extra-prec 3
|
|
353 (math-gammap1-raw (math-float n)))))))
|
|
354 ((equal n '(var inf var-inf)) n)
|
|
355 (t (calc-record-why 'numberp n)
|
41047
|
356 (list 'calcFunc-fact n)))))
|
40785
|
357
|
|
358 (math-defcache math-gamma-1q nil
|
|
359 (math-with-extra-prec 3
|
|
360 (math-gammap1-raw '(float -75 -2))))
|
|
361
|
|
362 (math-defcache math-gamma-3q nil
|
|
363 (math-with-extra-prec 3
|
|
364 (math-gammap1-raw '(float -25 -2))))
|
|
365
|
|
366 (defun math-factorial-iter (count n f)
|
|
367 (if (= (% n 5) 1)
|
|
368 (math-working (format "factorial(%d)" (1- n)) f))
|
|
369 (if (> count 0)
|
|
370 (math-factorial-iter (1- count) (1+ n) (math-mul n f))
|
41047
|
371 f))
|
40785
|
372
|
|
373 (defun calcFunc-dfact (n) ; [I I] [F F] [Public]
|
|
374 (cond ((Math-integer-negp n)
|
|
375 (if (math-oddp n)
|
|
376 (if (eq n -1)
|
|
377 1
|
|
378 (math-div (if (eq (math-mod n 4) 3) 1 -1)
|
|
379 (calcFunc-dfact (math-sub -2 n))))
|
|
380 (list 'calcFunc-dfact n)))
|
|
381 ((Math-zerop n) 1)
|
|
382 ((integerp n) (math-double-factorial-iter n (+ 2 (% n 2)) 1 2))
|
|
383 ((math-messy-integerp n)
|
|
384 (let ((temp (math-trunc n)))
|
|
385 (math-inexact-result)
|
|
386 (if (natnump temp)
|
|
387 (if (Math-lessp temp 200)
|
|
388 (math-with-extra-prec 1
|
|
389 (math-double-factorial-iter temp (+ 2 (% temp 2))
|
|
390 '(float 1 0) 2))
|
|
391 (let* ((half (math-div2 temp))
|
|
392 (even (math-mul (math-pow 2 half)
|
|
393 (calcFunc-fact (math-float half)))))
|
|
394 (if (math-evenp temp)
|
|
395 even
|
|
396 (math-div (calcFunc-fact n) even))))
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
397 (list 'calcFunc-dfact n))))
|
40785
|
398 ((equal n '(var inf var-inf)) n)
|
|
399 (t (calc-record-why 'natnump n)
|
41047
|
400 (list 'calcFunc-dfact n))))
|
40785
|
401
|
|
402 (defun math-double-factorial-iter (max n f step)
|
|
403 (if (< (% n 12) step)
|
|
404 (math-working (format "dfact(%d)" (- n step)) f))
|
|
405 (if (<= n max)
|
|
406 (math-double-factorial-iter max (+ n step) (math-mul n f) step)
|
41047
|
407 f))
|
40785
|
408
|
|
409 (defun calcFunc-perm (n m) ; [I I I] [F F F] [Public]
|
|
410 (cond ((and (integerp n) (integerp m) (<= m n) (>= m 0))
|
|
411 (math-factorial-iter m (1+ (- n m)) 1))
|
|
412 ((or (not (math-num-integerp n))
|
|
413 (and (math-messy-integerp n) (Math-lessp 100 n))
|
|
414 (not (math-num-integerp m))
|
|
415 (and (math-messy-integerp m) (Math-lessp 100 m)))
|
|
416 (or (math-realp n) (equal n '(var inf var-inf))
|
|
417 (math-reject-arg n 'realp))
|
|
418 (or (math-realp m) (equal m '(var inf var-inf))
|
|
419 (math-reject-arg m 'realp))
|
|
420 (and (math-num-integerp n) (math-negp n) (math-reject-arg n 'range))
|
|
421 (and (math-num-integerp m) (math-negp m) (math-reject-arg m 'range))
|
|
422 (math-div (calcFunc-fact n) (calcFunc-fact (math-sub n m))))
|
|
423 (t
|
|
424 (let ((tn (math-trunc n))
|
|
425 (tm (math-trunc m)))
|
|
426 (math-inexact-result)
|
|
427 (or (integerp tn) (math-reject-arg tn 'fixnump))
|
|
428 (or (integerp tm) (math-reject-arg tm 'fixnump))
|
|
429 (or (and (<= tm tn) (>= tm 0)) (math-reject-arg tm 'range))
|
|
430 (math-with-extra-prec 1
|
41047
|
431 (math-factorial-iter tm (1+ (- tn tm)) '(float 1 0)))))))
|
40785
|
432
|
|
433 (defun calcFunc-choose (n m) ; [I I I] [F F F] [Public]
|
|
434 (cond ((and (integerp n) (integerp m) (<= m n) (>= m 0))
|
|
435 (if (> m (/ n 2))
|
|
436 (math-choose-iter (- n m) n 1 1)
|
|
437 (math-choose-iter m n 1 1)))
|
|
438 ((not (math-realp n))
|
|
439 (math-reject-arg n 'realp))
|
|
440 ((not (math-realp m))
|
|
441 (math-reject-arg m 'realp))
|
|
442 ((not (math-num-integerp m))
|
|
443 (if (and (math-num-integerp n) (math-negp n))
|
|
444 (list 'calcFunc-choose n m)
|
|
445 (math-div (calcFunc-fact (math-float n))
|
|
446 (math-mul (calcFunc-fact m)
|
|
447 (calcFunc-fact (math-sub n m))))))
|
|
448 ((math-negp m) 0)
|
|
449 ((math-negp n)
|
|
450 (let ((val (calcFunc-choose (math-add (math-add n m) -1) m)))
|
|
451 (if (math-evenp (math-trunc m))
|
|
452 val
|
|
453 (math-neg val))))
|
|
454 ((and (math-num-integerp n)
|
|
455 (Math-lessp n m))
|
|
456 0)
|
|
457 (t
|
|
458 (math-inexact-result)
|
|
459 (let ((tm (math-trunc m)))
|
|
460 (or (integerp tm) (math-reject-arg tm 'fixnump))
|
|
461 (if (> tm 100)
|
|
462 (math-div (calcFunc-fact (math-float n))
|
|
463 (math-mul (calcFunc-fact (math-float m))
|
|
464 (calcFunc-fact (math-float
|
|
465 (math-sub n m)))))
|
|
466 (math-with-extra-prec 1
|
41047
|
467 (math-choose-float-iter tm n 1 1)))))))
|
40785
|
468
|
|
469 (defun math-choose-iter (m n i c)
|
|
470 (if (and (= (% i 5) 1) (> i 5))
|
|
471 (math-working (format "choose(%d)" (1- i)) c))
|
|
472 (if (<= i m)
|
|
473 (math-choose-iter m (1- n) (1+ i)
|
|
474 (math-quotient (math-mul c n) i))
|
41047
|
475 c))
|
40785
|
476
|
|
477 (defun math-choose-float-iter (count n i c)
|
|
478 (if (= (% i 5) 1)
|
|
479 (math-working (format "choose(%d)" (1- i)) c))
|
|
480 (if (> count 0)
|
|
481 (math-choose-float-iter (1- count) (math-sub n 1) (1+ i)
|
|
482 (math-div (math-mul c n) i))
|
41047
|
483 c))
|
40785
|
484
|
|
485
|
|
486 ;;; Stirling numbers.
|
|
487
|
|
488 (defun calcFunc-stir1 (n m)
|
41047
|
489 (math-stirling-number n m 1))
|
40785
|
490
|
|
491 (defun calcFunc-stir2 (n m)
|
41047
|
492 (math-stirling-number n m 0))
|
40785
|
493
|
41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
494 (defvar math-stirling-cache (vector [[1]] [[1]]))
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
495
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
496 ;; The variable math-stirling-local-cache is local to
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
497 ;; math-stirling-number, but is used by math-stirling-1
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
498 ;; and math-stirling-2, which are called by math-stirling-number.
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
499 (defvar math-stirling-local-cache)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
500
|
40785
|
501 (defun math-stirling-number (n m k)
|
|
502 (or (math-num-natnump n) (math-reject-arg n 'natnump))
|
|
503 (or (math-num-natnump m) (math-reject-arg m 'natnump))
|
|
504 (if (consp n) (setq n (math-trunc n)))
|
|
505 (or (integerp n) (math-reject-arg n 'fixnump))
|
|
506 (if (consp m) (setq m (math-trunc m)))
|
|
507 (or (integerp m) (math-reject-arg m 'fixnump))
|
|
508 (if (< n m)
|
|
509 0
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
510 (let ((math-stirling-local-cache (aref math-stirling-cache k)))
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
511 (while (<= (length math-stirling-local-cache) n)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
512 (let ((i (1- (length math-stirling-local-cache)))
|
40785
|
513 row)
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
514 (setq math-stirling-local-cache
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
515 (vconcat math-stirling-local-cache
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
516 (make-vector (length math-stirling-local-cache) nil)))
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
517 (aset math-stirling-cache k math-stirling-local-cache)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
518 (while (< (setq i (1+ i)) (length math-stirling-local-cache))
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
519 (aset math-stirling-local-cache i (setq row (make-vector (1+ i) nil)))
|
40785
|
520 (aset row 0 0)
|
|
521 (aset row i 1))))
|
|
522 (if (= k 1)
|
|
523 (math-stirling-1 n m)
|
41047
|
524 (math-stirling-2 n m)))))
|
40785
|
525
|
|
526 (defun math-stirling-1 (n m)
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
527 (or (aref (aref math-stirling-local-cache n) m)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
528 (aset (aref math-stirling-local-cache n) m
|
40785
|
529 (math-add (math-stirling-1 (1- n) (1- m))
|
41047
|
530 (math-mul (- 1 n) (math-stirling-1 (1- n) m))))))
|
40785
|
531
|
|
532 (defun math-stirling-2 (n m)
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
533 (or (aref (aref math-stirling-local-cache n) m)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
534 (aset (aref math-stirling-local-cache n) m
|
40785
|
535 (math-add (math-stirling-2 (1- n) (1- m))
|
41047
|
536 (math-mul m (math-stirling-2 (1- n) m))))))
|
40785
|
537
|
41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
538 (defvar math-random-table nil)
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
539 (defvar math-last-RandSeed nil)
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
540 (defvar math-random-ptr1 nil)
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
541 (defvar math-random-ptr2 nil)
|
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
diff
changeset
|
542 (defvar math-random-shift nil)
|
40785
|
543
|
|
544 ;;; Produce a random 10-bit integer, with (random) if no seed provided,
|
|
545 ;;; or else with Numerical Recipes algorithm ran3 / Knuth 3.2.2-A.
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
546
|
59009
|
547 (defvar var-RandSeed)
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
548 (defvar math-random-cache nil)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
549 (defvar math-gaussian-cache nil)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
550
|
40785
|
551 (defun math-init-random-base ()
|
59009
|
552 (if (and (boundp 'var-RandSeed) var-RandSeed)
|
40785
|
553 (if (eq (car-safe var-RandSeed) 'vec)
|
|
554 nil
|
|
555 (if (Math-integerp var-RandSeed)
|
|
556 (let* ((seed (math-sub 161803 var-RandSeed))
|
81742
3340bf730dea
(math-small-factorial-table,math-init-random-base,math-prime-test):
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
557 (mj (1+ (math-mod seed 1000000)))
|
3340bf730dea
(math-small-factorial-table,math-init-random-base,math-prime-test):
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
558 (mk (1+ (math-mod (math-quotient seed 1000000)
|
3340bf730dea
(math-small-factorial-table,math-init-random-base,math-prime-test):
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
559 1000000)))
|
40785
|
560 (i 0))
|
|
561 (setq math-random-table (cons 'vec (make-list 55 mj)))
|
|
562 (while (<= (setq i (1+ i)) 54)
|
|
563 (let* ((ii (% (* i 21) 55))
|
|
564 (p (nthcdr ii math-random-table)))
|
|
565 (setcar p mk)
|
|
566 (setq mk (- mj mk)
|
|
567 mj (car p)))))
|
|
568 (math-reject-arg var-RandSeed "*RandSeed must be an integer"))
|
|
569 (setq var-RandSeed (list 'vec var-RandSeed)
|
|
570 math-random-ptr1 math-random-table
|
|
571 math-random-cache nil
|
|
572 math-random-ptr2 (nthcdr 31 math-random-table))
|
|
573 (let ((i 200))
|
|
574 (while (> (setq i (1- i)) 0)
|
|
575 (math-random-base))))
|
|
576 (random t)
|
|
577 (setq var-RandSeed nil
|
|
578 math-random-cache nil
|
|
579 math-random-shift -4) ; assume RAND_MAX >= 16383
|
|
580 ;; This exercises the random number generator and also helps
|
|
581 ;; deduce a better value for RAND_MAX.
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
582 (let ((i 0))
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
583 (while (< (setq i (1+ i)) 30)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
584 (if (> (lsh (math-abs (random)) math-random-shift) 4095)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
585 (setq math-random-shift (1- math-random-shift))))))
|
40785
|
586 (setq math-last-RandSeed var-RandSeed
|
41047
|
587 math-gaussian-cache nil))
|
40785
|
588
|
|
589 (defun math-random-base ()
|
|
590 (if var-RandSeed
|
|
591 (progn
|
|
592 (setq math-random-ptr1 (or (cdr math-random-ptr1)
|
|
593 (cdr math-random-table))
|
|
594 math-random-ptr2 (or (cdr math-random-ptr2)
|
|
595 (cdr math-random-table)))
|
|
596 (logand (lsh (setcar math-random-ptr1
|
|
597 (logand (- (car math-random-ptr1)
|
|
598 (car math-random-ptr2)) 524287))
|
|
599 -6) 1023))
|
41047
|
600 (logand (lsh (random) math-random-shift) 1023)))
|
40785
|
601
|
|
602
|
|
603 ;;; Produce a random digit in the range 0..999.
|
|
604 ;;; Avoid various pitfalls that may lurk in the built-in (random) function!
|
|
605 ;;; Shuffling algorithm from Numerical Recipes, section 7.1.
|
59009
|
606 (defvar math-random-last)
|
81716
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
607 (defun math-random-three-digit-number ()
|
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
608 "Return a random three digit number."
|
59009
|
609 (let (i)
|
|
610 (or (and (boundp 'var-RandSeed) (eq var-RandSeed math-last-RandSeed))
|
40785
|
611 (math-init-random-base))
|
|
612 (or math-random-cache
|
|
613 (progn
|
|
614 (setq math-random-last (math-random-base)
|
|
615 math-random-cache (make-vector 13 nil)
|
|
616 i -1)
|
|
617 (while (< (setq i (1+ i)) 13)
|
|
618 (aset math-random-cache i (math-random-base)))))
|
|
619 (while (progn
|
|
620 (setq i (/ math-random-last 79) ; 0 <= i < 13
|
|
621 math-random-last (aref math-random-cache i))
|
|
622 (aset math-random-cache i (math-random-base))
|
|
623 (>= math-random-last 1000)))
|
41047
|
624 math-random-last))
|
40785
|
625
|
|
626 ;;; Produce an N-digit random integer.
|
|
627 (defun math-random-digits (n)
|
81716
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
628 "Produce a random N digit integer."
|
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
629 (let* ((slop (% (- 3 (% n 3)) 3))
|
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
630 (i (/ (+ n slop) 3))
|
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
631 (rnum 0))
|
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
632 (while (> i 0)
|
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
633 (setq rnum
|
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
634 (math-add
|
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
635 (math-random-three-digit-number)
|
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
636 (math-mul rnum 1000)))
|
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
637 (setq i (1- i)))
|
66a18d2f078a
(math-random-digit): Rename to `math-random-three-digit-number'.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
638 (math-normalize (math-scale-right rnum slop))))
|
40785
|
639
|
|
640 ;;; Produce a uniformly-distributed random float 0 <= N < 1.
|
|
641 (defun math-random-float ()
|
|
642 (math-make-float (math-random-digits calc-internal-prec)
|
41047
|
643 (- calc-internal-prec)))
|
40785
|
644
|
|
645 ;;; Produce a Gaussian-distributed random float with mean=0, sigma=1.
|
|
646 (defun math-gaussian-float ()
|
|
647 (math-with-extra-prec 2
|
|
648 (if (and math-gaussian-cache
|
|
649 (= (car math-gaussian-cache) calc-internal-prec))
|
|
650 (prog1
|
|
651 (cdr math-gaussian-cache)
|
|
652 (setq math-gaussian-cache nil))
|
|
653 (let* ((v1 (math-add (math-mul (math-random-float) 2) -1))
|
|
654 (v2 (math-add (math-mul (math-random-float) 2) -1))
|
|
655 (r (math-add (math-sqr v1) (math-sqr v2))))
|
|
656 (while (or (not (Math-lessp r 1)) (math-zerop r))
|
|
657 (setq v1 (math-add (math-mul (math-random-float) 2) -1)
|
|
658 v2 (math-add (math-mul (math-random-float) 2) -1)
|
|
659 r (math-add (math-sqr v1) (math-sqr v2))))
|
|
660 (let ((fac (math-sqrt (math-mul (math-div (calcFunc-ln r) r) -2))))
|
|
661 (setq math-gaussian-cache (cons calc-internal-prec
|
|
662 (math-mul v1 fac)))
|
41047
|
663 (math-mul v2 fac))))))
|
40785
|
664
|
|
665 ;;; Produce a random integer or real 0 <= N < MAX.
|
|
666 (defun calcFunc-random (max)
|
|
667 (cond ((Math-zerop max)
|
|
668 (math-gaussian-float))
|
|
669 ((Math-integerp max)
|
|
670 (let* ((digs (math-numdigs max))
|
|
671 (r (math-random-digits (+ digs 3))))
|
|
672 (math-mod r max)))
|
|
673 ((Math-realp max)
|
|
674 (math-mul (math-random-float) max))
|
|
675 ((and (eq (car max) 'intv) (math-constp max)
|
|
676 (Math-lessp (nth 2 max) (nth 3 max)))
|
|
677 (if (math-floatp max)
|
|
678 (let ((val (math-add (math-mul (math-random-float)
|
|
679 (math-sub (nth 3 max) (nth 2 max)))
|
|
680 (nth 2 max))))
|
|
681 (if (or (and (memq (nth 1 max) '(0 1)) ; almost not worth
|
|
682 (Math-equal val (nth 2 max))) ; checking!
|
|
683 (and (memq (nth 1 max) '(0 2))
|
|
684 (Math-equal val (nth 3 max))))
|
|
685 (calcFunc-random max)
|
|
686 val))
|
|
687 (let ((lo (if (memq (nth 1 max) '(0 1))
|
|
688 (math-add (nth 2 max) 1) (nth 2 max)))
|
|
689 (hi (if (memq (nth 1 max) '(1 3))
|
|
690 (math-add (nth 3 max) 1) (nth 3 max))))
|
|
691 (if (Math-lessp lo hi)
|
|
692 (math-add (calcFunc-random (math-sub hi lo)) lo)
|
|
693 (math-reject-arg max "*Empty interval")))))
|
|
694 ((eq (car max) 'vec)
|
|
695 (if (cdr max)
|
|
696 (nth (1+ (calcFunc-random (1- (length max)))) max)
|
|
697 (math-reject-arg max "*Empty list")))
|
|
698 ((and (eq (car max) 'sdev) (math-constp max) (Math-realp (nth 1 max)))
|
|
699 (math-add (math-mul (math-gaussian-float) (nth 2 max)) (nth 1 max)))
|
41047
|
700 (t (math-reject-arg max 'realp))))
|
40785
|
701
|
|
702 ;;; Choose N objects at random from the set MAX without duplicates.
|
|
703 (defun calcFunc-shuffle (n &optional max)
|
|
704 (or max (setq max n n -1))
|
|
705 (or (and (Math-num-integerp n)
|
|
706 (or (natnump (setq n (math-trunc n))) (eq n -1)))
|
|
707 (math-reject-arg n 'integerp))
|
|
708 (cond ((or (math-zerop max)
|
|
709 (math-floatp max)
|
|
710 (eq (car-safe max) 'sdev))
|
|
711 (if (< n 0)
|
|
712 (math-reject-arg n 'natnump)
|
|
713 (math-simple-shuffle n max)))
|
|
714 ((and (<= n 1) (>= n 0))
|
|
715 (math-simple-shuffle n max))
|
|
716 ((and (eq (car-safe max) 'intv) (math-constp max))
|
|
717 (let ((num (math-add (math-sub (nth 3 max) (nth 2 max))
|
|
718 (cdr (assq (nth 1 max)
|
|
719 '((0 . -1) (1 . 0)
|
|
720 (2 . 0) (3 . 1))))))
|
|
721 (min (math-add (nth 2 max) (if (memq (nth 1 max) '(0 1))
|
|
722 1 0))))
|
|
723 (if (< n 0) (setq n num))
|
|
724 (or (math-posp num) (math-reject-arg max 'range))
|
|
725 (and (Math-lessp num n) (math-reject-arg n 'range))
|
|
726 (if (Math-lessp n (math-quotient num 3))
|
|
727 (math-simple-shuffle n max)
|
|
728 (if (> (* n 4) (* num 3))
|
|
729 (math-add (math-sub min 1)
|
|
730 (math-shuffle-list n num (calcFunc-index num)))
|
|
731 (let ((tot 0)
|
|
732 (m 0)
|
|
733 (vec nil))
|
|
734 (while (< m n)
|
|
735 (if (< (calcFunc-random (- num tot)) (- n m))
|
|
736 (setq vec (cons (math-add min tot) vec)
|
|
737 m (1+ m)))
|
|
738 (setq tot (1+ tot)))
|
|
739 (math-shuffle-list n n (cons 'vec vec)))))))
|
|
740 ((eq (car-safe max) 'vec)
|
|
741 (let ((size (1- (length max))))
|
|
742 (if (< n 0) (setq n size))
|
|
743 (if (and (> n (/ size 2)) (<= n size))
|
|
744 (math-shuffle-list n size (copy-sequence max))
|
|
745 (let* ((vals (calcFunc-shuffle
|
|
746 n (list 'intv 3 1 (1- (length max)))))
|
|
747 (p vals))
|
|
748 (while (setq p (cdr p))
|
|
749 (setcar p (nth (car p) max)))
|
|
750 vals))))
|
|
751 ((math-integerp max)
|
|
752 (if (math-posp max)
|
|
753 (calcFunc-shuffle n (list 'intv 2 0 max))
|
|
754 (calcFunc-shuffle n (list 'intv 1 max 0))))
|
41047
|
755 (t (math-reject-arg max 'realp))))
|
40785
|
756
|
|
757 (defun math-simple-shuffle (n max)
|
|
758 (let ((vec nil)
|
|
759 val)
|
|
760 (while (>= (setq n (1- n)) 0)
|
|
761 (while (math-member (setq val (calcFunc-random max)) vec))
|
|
762 (setq vec (cons val vec)))
|
41047
|
763 (cons 'vec vec)))
|
40785
|
764
|
|
765 (defun math-shuffle-list (n size vec)
|
|
766 (let ((j size)
|
|
767 k temp
|
|
768 (p vec))
|
|
769 (while (cdr (setq p (cdr p)))
|
|
770 (setq k (calcFunc-random j)
|
|
771 j (1- j)
|
|
772 temp (nth k p))
|
|
773 (setcar (nthcdr k p) (car p))
|
|
774 (setcar p temp))
|
41047
|
775 (cons 'vec (nthcdr (- size n -1) vec))))
|
40785
|
776
|
|
777 (defun math-member (x list)
|
|
778 (while (and list (not (equal x (car list))))
|
|
779 (setq list (cdr list)))
|
41047
|
780 list)
|
40785
|
781
|
|
782
|
|
783 ;;; Check if the integer N is prime. [X I]
|
|
784 ;;; Return (nil) if non-prime,
|
|
785 ;;; (nil N) if non-prime with known factor N,
|
|
786 ;;; (nil unknown) if non-prime with no known factors,
|
|
787 ;;; (t) if prime,
|
|
788 ;;; (maybe N P) if probably prime (after N iters with probability P%)
|
58145
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
789 (defvar math-prime-test-cache '(-1))
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
790
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
791 (defvar math-prime-test-cache-k)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
792 (defvar math-prime-test-cache-q)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
793 (defvar math-prime-test-cache-nm1)
|
be3e73786548
(math-prime-factors-finished): Declared it as a variable.
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
794
|
40785
|
795 (defun math-prime-test (n iters)
|
|
796 (if (and (Math-vectorp n) (cdr n))
|
|
797 (setq n (nth (1- (length n)) n)))
|
|
798 (if (Math-messy-integerp n)
|
|
799 (setq n (math-trunc n)))
|
|
800 (let ((res))
|
|
801 (while (> iters 0)
|
|
802 (setq res
|
|
803 (cond ((and (integerp n) (<= n 5003))
|
|
804 (list (= (math-next-small-prime n) n)))
|
|
805 ((not (Math-integerp n))
|
|
806 (error "Argument must be an integer"))
|
|
807 ((Math-integer-negp n)
|
|
808 '(nil))
|
81742
3340bf730dea
(math-small-factorial-table,math-init-random-base,math-prime-test):
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
809 ((Math-natnum-lessp n 8000000)
|
40785
|
810 (setq n (math-fixnum n))
|
|
811 (let ((i -1) v)
|
|
812 (while (and (> (% n (setq v (aref math-primes-table
|
|
813 (setq i (1+ i)))))
|
|
814 0)
|
|
815 (< (* v v) n)))
|
|
816 (if (= (% n v) 0)
|
|
817 (list nil v)
|
|
818 '(t))))
|
|
819 ((not (equal n (car math-prime-test-cache)))
|
|
820 (cond ((= (% (nth 1 n) 2) 0) '(nil 2))
|
|
821 ((= (% (nth 1 n) 5) 0) '(nil 5))
|
81594
5f2dd5474f8a
(math-init-random-base,math-prime-test): Use math-read-number-simple
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
822 (t (let ((q n) (sum 0))
|
5f2dd5474f8a
(math-init-random-base,math-prime-test): Use math-read-number-simple
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
823 (while (not (eq q 0))
|
5f2dd5474f8a
(math-init-random-base,math-prime-test): Use math-read-number-simple
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
824 (setq sum (%
|
5f2dd5474f8a
(math-init-random-base,math-prime-test): Use math-read-number-simple
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
825 (+
|
5f2dd5474f8a
(math-init-random-base,math-prime-test): Use math-read-number-simple
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
826 sum
|
5f2dd5474f8a
(math-init-random-base,math-prime-test): Use math-read-number-simple
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
827 (calcFunc-mod
|
81742
3340bf730dea
(math-small-factorial-table,math-init-random-base,math-prime-test):
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
828 q 1000000))
|
81594
5f2dd5474f8a
(math-init-random-base,math-prime-test): Use math-read-number-simple
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
829 111111))
|
5f2dd5474f8a
(math-init-random-base,math-prime-test): Use math-read-number-simple
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
830 (setq q
|
5f2dd5474f8a
(math-init-random-base,math-prime-test): Use math-read-number-simple
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
831 (math-quotient
|
81742
3340bf730dea
(math-small-factorial-table,math-init-random-base,math-prime-test):
Jay Belanger <jay.p.belanger@gmail.com>
diff
changeset
|
832 q 1000000)))
|
40785
|
833 (cond ((= (% sum 3) 0) '(nil 3))
|
|
834 ((= (% sum 7) 0) '(nil 7))
|
|
835 ((= (% sum 11) 0) '(nil 11))
|
|
836 ((= (% sum 13) 0) '(nil 13))
|
|
837 ((= (% sum 37) 0) '(nil 37))
|
|
838 (t
|
|
839 (setq math-prime-test-cache-k 1
|
|
840 math-prime-test-cache-q
|
|
841 (math-div2 n)
|
|
842 math-prime-test-cache-nm1
|
|
843 (math-add n -1))
|
|
844 (while (math-evenp
|
|
845 math-prime-test-cache-q)
|
|
846 (setq math-prime-test-cache-k
|
|
847 (1+ math-prime-test-cache-k)
|
|
848 math-prime-test-cache-q
|
|
849 (math-div2
|
|
850 math-prime-test-cache-q)))
|
|
851 (setq iters (1+ iters))
|
|
852 (list 'maybe
|
|
853 0
|
|
854 (math-sub
|
|
855 100
|
|
856 (math-div
|
|
857 '(float 232 0)
|
|
858 (math-numdigs n))))))))))
|
|
859 ((not (eq (car (nth 1 math-prime-test-cache)) 'maybe))
|
|
860 (nth 1 math-prime-test-cache))
|
|
861 (t ; Fermat step
|
|
862 (let* ((x (math-add (calcFunc-random (math-add n -2)) 2))
|
|
863 (y (math-pow-mod x math-prime-test-cache-q n))
|
|
864 (j 0))
|
|
865 (while (and (not (eq y 1))
|
|
866 (not (equal y math-prime-test-cache-nm1))
|
|
867 (< (setq j (1+ j)) math-prime-test-cache-k))
|
|
868 (setq y (math-mod (math-mul y y) n)))
|
|
869 (if (or (equal y math-prime-test-cache-nm1)
|
|
870 (and (eq y 1) (eq j 0)))
|
|
871 (list 'maybe
|
|
872 (1+ (nth 1 (nth 1 math-prime-test-cache)))
|
|
873 (math-mul (nth 2 (nth 1 math-prime-test-cache))
|
|
874 '(float 25 -2)))
|
|
875 '(nil unknown))))))
|
|
876 (setq math-prime-test-cache (list n res)
|
|
877 iters (if (eq (car res) 'maybe)
|
|
878 (1- iters)
|
|
879 0)))
|
41047
|
880 res))
|
40785
|
881
|
|
882 (defun calcFunc-prime (n &optional iters)
|
|
883 (or (math-num-integerp n) (math-reject-arg n 'integerp))
|
|
884 (or (not iters) (math-num-integerp iters) (math-reject-arg iters 'integerp))
|
|
885 (if (car (math-prime-test (math-trunc n) (math-trunc (or iters 1))))
|
|
886 1
|
41047
|
887 0))
|
40785
|
888
|
|
889 ;;; Theory: summing base-10^6 digits modulo 111111 is "casting out 999999s".
|
|
890 ;;; Initial probability that N is prime is 1/ln(N) = log10(e)/log10(N).
|
|
891 ;;; After culling [2,3,5,7,11,13,37], probability of primality is 5.36 x more.
|
|
892 ;;; Initial reported probability of non-primality is thus 100% - this.
|
|
893 ;;; Each Fermat step multiplies this probability by 25%.
|
|
894 ;;; The Fermat step is algorithm P from Knuth section 4.5.4.
|
|
895
|
|
896
|
|
897 (defun calcFunc-prfac (n)
|
|
898 (setq math-prime-factors-finished t)
|
|
899 (if (Math-messy-integerp n)
|
|
900 (setq n (math-trunc n)))
|
|
901 (if (Math-natnump n)
|
|
902 (if (Math-natnum-lessp 2 n)
|
|
903 (let (factors res p (i 0))
|
|
904 (while (and (not (eq n 1))
|
|
905 (< i (length math-primes-table)))
|
|
906 (setq p (aref math-primes-table i))
|
|
907 (while (eq (cdr (setq res (cond ((eq n p) (cons 1 0))
|
|
908 ((eq n 1) (cons 0 1))
|
|
909 ((consp n) (math-idivmod n p))
|
|
910 (t (cons (/ n p) (% n p))))))
|
|
911 0)
|
|
912 (math-working "factor" p)
|
|
913 (setq factors (nconc factors (list p))
|
|
914 n (car res)))
|
|
915 (or (eq n 1)
|
|
916 (Math-natnum-lessp p (car res))
|
|
917 (setq factors (nconc factors (list n))
|
|
918 n 1))
|
|
919 (setq i (1+ i)))
|
|
920 (or (setq math-prime-factors-finished (eq n 1))
|
|
921 (setq factors (nconc factors (list n))))
|
|
922 (cons 'vec factors))
|
|
923 (list 'vec n))
|
|
924 (if (Math-integerp n)
|
|
925 (if (eq n -1)
|
|
926 (list 'vec n)
|
|
927 (cons 'vec (cons -1 (cdr (calcFunc-prfac (math-neg n))))))
|
|
928 (calc-record-why 'integerp n)
|
41047
|
929 (list 'calcFunc-prfac n))))
|
40785
|
930
|
|
931 (defun calcFunc-totient (n)
|
|
932 (if (Math-messy-integerp n)
|
|
933 (setq n (math-trunc n)))
|
|
934 (if (Math-natnump n)
|
|
935 (if (Math-natnum-lessp n 2)
|
|
936 (if (Math-negp n)
|
|
937 (calcFunc-totient (math-abs n))
|
|
938 n)
|
|
939 (let ((factors (cdr (calcFunc-prfac n)))
|
|
940 p)
|
|
941 (if math-prime-factors-finished
|
|
942 (progn
|
|
943 (while factors
|
|
944 (setq p (car factors)
|
|
945 n (math-mul (math-div n p) (math-add p -1)))
|
|
946 (while (equal p (car factors))
|
|
947 (setq factors (cdr factors))))
|
|
948 n)
|
|
949 (calc-record-why "*Number too big to factor" n)
|
|
950 (list 'calcFunc-totient n))))
|
|
951 (calc-record-why 'natnump n)
|
41047
|
952 (list 'calcFunc-totient n)))
|
40785
|
953
|
|
954 (defun calcFunc-moebius (n)
|
|
955 (if (Math-messy-integerp n)
|
|
956 (setq n (math-trunc n)))
|
|
957 (if (and (Math-natnump n) (not (eq n 0)))
|
|
958 (if (Math-natnum-lessp n 2)
|
|
959 (if (Math-negp n)
|
|
960 (calcFunc-moebius (math-abs n))
|
|
961 1)
|
|
962 (let ((factors (cdr (calcFunc-prfac n)))
|
|
963 (mu 1))
|
|
964 (if math-prime-factors-finished
|
|
965 (progn
|
|
966 (while factors
|
|
967 (setq mu (if (equal (car factors) (nth 1 factors))
|
|
968 0 (math-neg mu))
|
|
969 factors (cdr factors)))
|
|
970 mu)
|
|
971 (calc-record-why "Number too big to factor" n)
|
|
972 (list 'calcFunc-moebius n))))
|
|
973 (calc-record-why 'posintp n)
|
41047
|
974 (list 'calcFunc-moebius n)))
|
40785
|
975
|
|
976
|
|
977 (defun calcFunc-nextprime (n &optional iters)
|
|
978 (if (Math-integerp n)
|
|
979 (if (Math-integer-negp n)
|
|
980 2
|
|
981 (if (and (integerp n) (< n 5003))
|
|
982 (math-next-small-prime (1+ n))
|
|
983 (if (math-evenp n)
|
|
984 (setq n (math-add n -1)))
|
|
985 (let (res)
|
|
986 (while (not (car (setq res (math-prime-test
|
|
987 (setq n (math-add n 2))
|
|
988 (or iters 1))))))
|
|
989 (if (and calc-verbose-nextprime
|
|
990 (eq (car res) 'maybe))
|
|
991 (calc-report-prime-test res)))
|
|
992 n))
|
|
993 (if (Math-realp n)
|
|
994 (calcFunc-nextprime (math-trunc n) iters)
|
41047
|
995 (math-reject-arg n 'integerp))))
|
40785
|
996
|
|
997 (defun calcFunc-prevprime (n &optional iters)
|
|
998 (if (Math-integerp n)
|
|
999 (if (Math-lessp n 4)
|
|
1000 2
|
|
1001 (if (math-evenp n)
|
|
1002 (setq n (math-add n 1)))
|
|
1003 (let (res)
|
|
1004 (while (not (car (setq res (math-prime-test
|
|
1005 (setq n (math-add n -2))
|
|
1006 (or iters 1))))))
|
|
1007 (if (and calc-verbose-nextprime
|
|
1008 (eq (car res) 'maybe))
|
|
1009 (calc-report-prime-test res)))
|
|
1010 n)
|
|
1011 (if (Math-realp n)
|
|
1012 (calcFunc-prevprime (math-ceiling n) iters)
|
41047
|
1013 (math-reject-arg n 'integerp))))
|
40785
|
1014
|
|
1015 (defun math-next-small-prime (n)
|
|
1016 (if (and (integerp n) (> n 2))
|
|
1017 (let ((lo -1)
|
|
1018 (hi (length math-primes-table))
|
|
1019 mid)
|
|
1020 (while (> (- hi lo) 1)
|
|
1021 (if (> n (aref math-primes-table
|
|
1022 (setq mid (ash (+ lo hi) -1))))
|
|
1023 (setq lo mid)
|
|
1024 (setq hi mid)))
|
|
1025 (aref math-primes-table hi))
|
41047
|
1026 2))
|
40785
|
1027
|
58649
|
1028 (provide 'calc-comb)
|
40785
|
1029
|
93975
|
1030 ;; arch-tag: 1d75ee9b-0815-42bd-a321-bb3dc001cc02
|
41047
|
1031 ;;; calc-comb.el ends here
|