Commit 015558a3df8fae2686c7b6e2f049da3f4b1a2cd8

Authored by Raymond Toy
1 parent a5a4c7ac

Fix a divide by zero error in s-exp-integral-e for v = 1. We need to

skip the first term in the series.
Showing 1 changed file with 8 additions and 5 deletions
qd-gamma.lisp
... ... @@ -257,7 +257,7 @@
257 257 ;; This works ok, but has problems for z > 3 where sometimes the
258 258 ;; result is greater than 1.
259 259 #+nil
260   -(defun erf (z)
  260 +(defun cf-erf (z)
261 261 (let* ((z2 (* z z))
262 262 (twoz2 (* 2 z2)))
263 263 (* (/ (* 2 z)
... ... @@ -504,10 +504,13 @@
504 504 (- (psi v) (log z)))
505 505 (loop for k from 0
506 506 for term = 1 then (* term (/ -z k))
507   - for sum = (/ (- 1 v)) then (+ sum (let ((denom (- k n-1)))
508   - (if (zerop denom)
509   - 0
510   - (/ term denom))))
  507 + for sum = (if (zerop n-1)
  508 + 0
  509 + (/ (- 1 v)))
  510 + then (+ sum (let ((denom (- k n-1)))
  511 + (if (zerop denom)
  512 + 0
  513 + (/ term denom))))
511 514 when (< (abs term) (* (abs sum) eps))
512 515 return sum)))
513 516 (loop for k from 0
... ...