### Fix #18: Add more tests for incomplete gamma functions

```Add more tests using the data provided in
#9 (comment 7019).

We first need convert the original test data into a form that we can
read by converting the numbers to #q form.```
parent 0f3c4ea4
 ... ... @@ -400,14 +400,44 @@ (setf ans (add-qd-d ans p))))) (defun make-float (sign int-part frac-part scale exp) "Create a quad-float given the SIGN, INT-PART, FRAC-PART, SCALE, and EXP. For example, if we're reading #q-123.45678q100, MAKE-FLOAT is called like (MAKE-FLOAT -1 123 45678 5 100)" (declare (type (member -1 1) sign) (type unsigned-byte int-part frac-part) (fixnum scale exp)) (let ((qd (rational-to-qd (* (+ int-part (/ frac-part (expt 10 scale))) (expt 10 exp))))) (if (minusp sign) (neg-qd qd) qd))) ;; Quick exit when the int part and frac part are zero. (when (and (zerop int-part) (zerop frac-part)) (return-from make-float (make-qd-d (* sign 0d0)))) ;; Let i, f, s, e be int-part, frac-part, scale and exp, respectively. Compute the ;; rational number ;; ;; r = (i+(f/10^s))*10^e = (10^s*i+f)/10^s*10^e ;; ;; Take base 10 log: ;; ;; log10(r) = log10(10^s*i+f) + e - s ;; ;; This will underflow if log10(r) <= ;; log10(least-positive-double-float) and overflow if log10(r) >= ;; log10(most-positive-double-float). Let's catch these cases early ;; so that we don't end up computing huge rationals for 10^(e-s). (let* ((log-res (+ (log (+ frac-part (* int-part (expt 10 scale))) 10) (- exp scale))) (qd (cond ((< log-res (log least-positive-double-float 10)) ;; Underflow. Silently set to 0. (rational-to-qd 0)) ((> log-res (log most-positive-double-float 10)) ;; Overflow. (error 'floating-point-overflow)) (t (rational-to-qd (* (+ int-part (/ frac-part (expt 10 scale))) (expt 10 exp))))))) (if (minusp sign) (neg-qd qd) qd))) ;; This seems to work, but really needs to be rewritten! (defun read-qd (stream) ... ...