-
Notifications
You must be signed in to change notification settings - Fork 0
/
belln.lisp
35 lines (32 loc) · 1.35 KB
/
belln.lisp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
(defun belln-cutoff (n)
(let* ((bn n) (log_bn (if (zerop bn) 0 (log bn)))
(log_pow (* n log_bn)) (log_fac (* bn (- log_bn 1))))
(loop while (>= (- log_pow log_fac) -2)
do (setf bn (1+ bn) log_bn (log bn)
log_pow (* n log_bn) log_fac (+ log_fac log_bn)))
bn))
(defun belln-bsplit (i j n)
(let (aa pp qq tt)
(if (= (- j i) 1)
(setq aa (expt i n) pp 1 qq i tt aa)
(let ((m (ash (+ i j) -1)))
(multiple-value-bind (tl ql pl) (belln-bsplit i m n)
(multiple-value-bind (tr qr pr) (belln-bsplit m j n)
(setq pp (* pl pr)
qq (* ql qr)
tt (+ (* qr tl) (* pl tr)))))))
(values tt qq pp)))
(defun integer-belln (n)
(let ((bn (belln-cutoff n)))
(multiple-value-bind (tt qq) (belln-bsplit 1 (1+ bn) n)
(multiple-value-bind (te qe) (split-taylor-e 0 (taylor-e-size (- (integer-length tt) (integer-length qq) -12)))
(values (round (* tt qe) (* qq te)))))))
(defun integer-belln-triangle (n)
(let ((up (make-array (1+ n) :initial-element 1))
(down (make-array (1+ n) :initial-element 1)))
(loop for i from 1 below n do
(loop for j from 1 to i do
(setf (aref down j) (+ (aref down (- j 1)) (aref up (- j 1)))))
(rotatef up down)
(setf (aref down 0) (aref up i)))
(aref down 0)))