74 lines
1.8 KiB
Text
74 lines
1.8 KiB
Text
(load "@lib/frac.l")
|
|
|
|
(de fact (N)
|
|
(cache '(NIL) N
|
|
(if (=0 N) 1 (apply * (range 1 N))) ) )
|
|
|
|
(de binomial (N K)
|
|
(frac
|
|
(/
|
|
(fact N)
|
|
(* (fact (- N K)) (fact K)) )
|
|
1 ) )
|
|
|
|
(de A (N M)
|
|
(let Sum (0 . 1)
|
|
(for X M
|
|
(setq Sum
|
|
(f+
|
|
Sum
|
|
(f*
|
|
(binomial (+ N 3) (- N (* X 6)))
|
|
(berno (- N (* X 6)) ) ) ) ) )
|
|
Sum ) )
|
|
|
|
(de berno (N)
|
|
(cache '(NIL) N
|
|
(cond
|
|
((=0 N) (1 . 1))
|
|
((= 1 N) (-1 . 2))
|
|
((bit? 1 N) (0 . 1))
|
|
(T
|
|
(case (% N 6)
|
|
(0
|
|
(f/
|
|
(f-
|
|
(frac (+ N 3) 3)
|
|
(A N (/ N 6)) )
|
|
(binomial (+ N 3) N) ) )
|
|
(2
|
|
(f/
|
|
(f-
|
|
(frac (+ N 3) 3)
|
|
(A N (/ (- N 2) 6)) )
|
|
(binomial (+ N 3) N) ) )
|
|
(4
|
|
(f/
|
|
(f-
|
|
(f* (-1 . 1) (frac (+ N 3) 6))
|
|
(A N (/ (- N 4) 6)) )
|
|
(binomial (+ N 3) N) ) ) ) ) ) ) )
|
|
|
|
(de berno-brute (N)
|
|
(cache '(NIL) N
|
|
(let Sum (0 . 1)
|
|
(cond
|
|
((=0 N) (1 . 1))
|
|
((= 1 N) (-1 . 2))
|
|
((bit? 1 N) (0 . 1))
|
|
(T
|
|
(for (X 0 (> N X) (inc X))
|
|
(setq Sum
|
|
(f+
|
|
Sum
|
|
(f* (binomial (inc N) X) (berno-brute X)) ) ) )
|
|
(f/ (f* (-1 . 1) Sum) (binomial (inc N) N)) ) ) ) ) )
|
|
|
|
(for (N 0 (> 62 N) (inc N))
|
|
(if (or (= N 1) (not (bit? 1 N)))
|
|
(tab (2 4 -60) N " => " (sym (berno N))) ) )
|
|
|
|
(for (N 0 (> 400 N) (inc N))
|
|
(test (berno N) (berno-brute N)) )
|
|
|
|
(bye)
|