; sums of powers
(define (range . args)
(case (length args)
((1) (range 0 (car args) (if (negative? (car args)) -1 1)))
((2) (range (car args) (cadr args) (if (< (car args) (cadr args)) 1 -1)))
((3) (let ((le? (if (negative? (caddr args)) >= <=)))
(let loop ((x(car args)) (xs '()))
(if (le? (cadr args) x)
(reverse xs)
(loop (+ x (caddr args)) (cons x xs))))))
(else (error 'range "unrecognized arguments"))))
(define (but-last xs)
(let loop ((xs xs) (zs '()))
(if (null? (cdr xs)) (reverse zs)
(loop (cdr xs) (cons (car xs) zs)))))
(define (a i j)
(if (zero? i) (/ (+ j 1))
(* (+ j 1) (- (a (- i 1) j) (a (- i 1) (+ j 1))))))
(define (b n) (a (- n 1) 0))
(define (akiyama-tanigawa limit)
(let loop ((i 0) (bs '()) (xs (map / (range 1 (+ limit 2)))))
(if (null? (cdr xs)) (reverse bs)
(loop (+ i 1) (cons (car xs) bs)
(map * (range 1 (- limit i -1))
(map - (but-last xs) (cdr xs)))))))
(display (b 19)) (newline)
(display (akiyama-tanigawa 19)) (newline)
(define (sum-powers m n)
(do ((i 1 (+ i 1)) (s 0 (+ s (expt i m)))) ((< n i) s)))
(define (choose n k)
(if (zero? k) 1
(* n (/ k) (choose (- n 1) (- k 1)))))
(define (bernoulli-formula m n)
(let ((m+1 (+ m 1)))
(let loop ((k 0) (bs (akiyama-tanigawa m+1)) (s 0))
(if (< m k) (/ s m+1)
(loop (+ k 1) (cdr bs)
(+ s (* (choose m+1 k) (car bs)
(expt n (- m+1 k)))))))))
(display (sum-powers 10 1000)) (newline)
(display (bernoulli-formula 10 1000)) (newline)