Commit fe8cffb5ee8e7161addf7586ecaee00682c6bf1b

Authored by Raymond Toy
1 parent e10e402d

Clean up implemenation of s-exp-integral-e a bit.

Showing 1 changed file with 8 additions and 7 deletions
qd-gamma.lisp
... ... @@ -497,18 +497,19 @@
497 497 (if (and (realp v)
498 498 (= v (ftruncate v)))
499 499 ;; v is an integer
500   - (let ((n (truncate v)))
501   - (- (* (/ (expt -z (- v 1))
  500 + (let* ((n (truncate v))
  501 + (n-1 (1- n)))
  502 + (- (* (/ (expt -z n-1)
502 503 (gamma v))
503 504 (- (psi v) (log z)))
504 505 (loop for k from 0
505   - for term = 1 then (* term (/ -z k))
506   - for sum = (/ (- 1 v)) then (+ sum (let ((denom (+ k 1 (- n))))
  506 + for term = 1 then (* term (/ -z k))
  507 + for sum = (/ (- 1 v)) then (+ sum (let ((denom (- k n-1)))
507 508 (if (zerop denom)
508 509 0
509   - (/ term (+ k 1 -v)))))
510   - when (< (abs term) (* (abs sum) eps))
511   - return sum)))
  510 + (/ term denom))))
  511 + when (< (abs term) (* (abs sum) eps))
  512 + return sum)))
512 513 (loop for k from 0
513 514 for term = 1 then (* term (/ -z k))
514 515 for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
... ...