; 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)

43867/798
(1 1/2 1/6 0 -1/30 0 1/42 0 -1/30 0 5/66 0 -691/2730 0 7/6 0 -3617/510 0 43867/798)
91409924241424243424241924242500
91409924241424243424241924242500