Index: qd-gamma.lisp
===================================================================
--- qd-gamma.lisp	(revision 46e7193fb5b2fad35e67366ff375d9ebbb524c5c)
+++ qd-gamma.lisp	(revision d17228b86105b8a6ac652459b1100266ad0311b3)
@@ -293,5 +293,7 @@
 	    (cf-incomplete-gamma a z)
 	    (- (gamma a) (cf-incomplete-gamma-tail a z)))
-	(cf-incomplete-gamma a z))))
+	(if (< (abs z) (abs a))
+	    (cf-incomplete-gamma a z)
+	    (- (gamma a) (cf-incomplete-gamma-tail a z))))))
 
 (defun erf (z)
@@ -342,2 +344,47 @@
   (* (expt z (- v 1))
      (incomplete-gamma-tail (- 1 v) z)))
+
+(defun fresnel-s (z)
+  "Fresnel S:
+
+   S(z) = integrate(sin(%pi*t^2/2), t, 0, z) "
+  (let ((sqrt-pi (sqrt (float-pi z))))
+    (flet ((fs (z)
+	     ;; Wolfram gives
+	     ;;
+	     ;;  S(z) = (1+%i)/4*(erf(c*z) - %i*erf(conjugate(c)*z))
+	     ;;
+	     ;; where c = sqrt(%pi)/2*(1+%i).
+	     (* #c(1/4 1/4)
+		(- (erf (* #c(1/2 1/2) sqrt-pi z))
+		   (* #c(0 1)
+		      (erf (* #c(1/2 -1/2) sqrt-pi z)))))))
+      (if (realp z)
+	  ;; FresnelS is real for a real argument. And it is odd.
+	  (if (minusp z)
+	      (- (realpart (fs (- z))))
+	      (realpart (fs z)))
+	  (fs z)))))
+
+(defun fresnel-c (z)
+  "Fresnel C:
+
+   C(z) = integrate(cos(%pi*t^2/2), t, 0, z) "
+  (let ((sqrt-pi (sqrt (float-pi z))))
+    (flet ((fs (z)
+	     ;; Wolfram gives
+	     ;;
+	     ;;  C(z) = (1-%i)/4*(erf(c*z) + %i*erf(conjugate(c)*z))
+	     ;;
+	     ;; where c = sqrt(%pi)/2*(1+%i).
+	     (* #c(1/4 -1/4)
+		(+ (erf (* #c(1/2 1/2) sqrt-pi z))
+		   (* #c(0 1)
+		      (erf (* #c(1/2 -1/2) sqrt-pi z)))))))
+      (if (realp z)
+	  ;; FresnelS is real for a real argument. And it is odd.
+	  (if (minusp z)
+	      (- (realpart (fs (- z))))
+	      (realpart (fs z)))
+	  (fs z)))))
+