;; 1f1.lisp
;;
;;  Confluent hypergeometric function M(a;b;z) or 1F1(a;b;z)
;;
;;  The confluent hypergeometric function is the solution to
;;  the differential equation:
;;
;;    z M"(a;b;z) + (b-z) M'(a;b;z) - a M(a;b;z) = 0    (A&S 13.1.1)
;;
;;  Evaluate it by direct summation of A&S 13.1.2
;;
;;                  a z    (a)_2 z^2          (a)_n z^n
;;   M(a;b;z) = 1 + --- +  ---------  + ... + --------- + ...
;;                   b     (b)_2 2!           (b)_n n!
;;
;;   where (a)_n is Pochhammer symbol (A&S 6.1.22)
;;         (a)_0 = 1
;;         (a)_n = a(a+1)(a+2)...(a+n-1),   
;;
;;   FIXME:  1F1 has singularities when b is a negative integer
;;           These are not checked and an error will occur. 

(in-package :maxima)

;; Confluent hypergeometric function M(a;b;z) or 1F1(a;b;z)
;; Arguments are (complex) maxima numbers and are not checked.  
;; Result is a (complex) maxima float and should be accurate
;; to machine precision.
(defun onefone (a b z)
  (let ((fpprec machine-mantissa-precision))
    ($float (onefonebf a b z))))

;; Confluent hypergeometric function M(a;b;z) or 1F1(a;b;z)
;; Arguments are complex maxima numbers and are not checked.
;; Result is a complex maxima bigfloat and should be accurate
;; to fpprec bits of precision.
(defun onefonebf (a b z)
  (let ( (old-fpprec fpprec) result one-f-one-fp (bits-accuracy 0))
    ;; Keep increasing precision until answer is sufficiently accurate
    (do ((fpprec (+ old-fpprec 9) (+ fpprec 64)))
	((> bits-accuracy old-fpprec))
	(multiple-value-setq 
	  (one-f-one bits-accuracy)
	  (one-f-one-fp (floattofpz a) (floattofpz b)(floattofpz z)))
	;; Reformat as bigfloat to capture current fpprec
	(setq one-f-one (fpztobigfloat one-f-one)))
    (setq fpprec old-fpprec)
    ($bfloat one-f-one)))

;; Confluent hypergeometric function M(a;b;z) or 1F1(a;b;z)
;; by direct summation of A&S 13.1.2
;;
;; Arguments are complex FP numbers represented as list
;; ( (mantissa exponent) (mantissa exponent))
;;
;;                  a z    (a)_2 z^2          (a)_n z^n
;;   M(a;b;z) = 1 + --- +  ---------  + ... + --------- + ...
;;                   b     (b)_2 2!           (b)_n n!
;;
;;   where (a)_n = a(a+1)(a+2)...(a+n-1),   
;;         (a)_0 = 1
;;
;;   nth partial sum s_n = s_(n-1) + t_n
;;   with
;;         t_1 = a z/b
;;         t_n = t_(n-1) * (a+n-1)*z/((b+n-1)*n)
;;
;;   To reduce number of (complex) divisions, accumulate 
;;   numerator and denominator separately.
;;
;;   s_n = N_n / D_n,
;;
;;   D_n = (b)_n n!, the denominator of t_n
;;   D_0 = 1
;;   D_n = D_(n-1) * (b+n-1)*n
;;
;;   N_0 = 1
;;   N_n = N_(n-1) * (b+n-1)*n + T_n
;;
;;   T_0 = 1
;;   T_n = T_(n-1) * (a+n-1)*z
;;
;;  Can have significant loss of precision due to cancellation of
;;  large terms.  Estimate the number of bits of precision lost
;;  by comparing exponents of largest partial sum s_n and final s_n.
;;
;;  For example, if a and b are equal and z = 0 + 140 i then 
;;  intermediate terms are of order 10^60 ~ 2^200 but result is 
;;  exp(140 i) = cos(140) + i sin(140) ~ O(1)
;;
;;  Exponent of t_n ~ (exponent of T_n) - (exponent of D_n)
;;  Complex numbers are a complication.  Just take larger of two values
;;
(defun one-f-one-fp (a b z)
  (prog
    ((sum-numerator   (floattofpz 1))
     (sum-denominator (floattofpz 1))
     (term-numerator  (floattofpz 1))
     (sum-exponent 0) ;; The first term is 1
     (max-sum-exponent 0)
     bits-accuracy ;; estimate of number of accuracy bits in result
     aterm bterm answer n-1-fp )
    (do ((n 1 (f+ n 1)))
	;; Iterate until the last term is insignificant when compared 
	;; to the largest partial sum
	((> (- max-sum-exponent fpprec) ;; Least significant bit in partial sum
	   (- (fpz-exponent term-numerator) (fpz-exponent sum-denominator) 1)))
	(setq n-1-fp (floattofpz (f- n 1)))
	(setq bterm (fpz* (floattofpz n) (fpz+ b n-1-fp)))
	(setq sum-denominator (fpz* sum-denominator bterm))
	(setq sum-numerator (fpz* sum-numerator bterm))
	(setq aterm (fpz* z (fpz+ a n-1-fp)))
	(setq term-numerator (fpz* aterm term-numerator))
	(setq sum-numerator (fpz+ sum-numerator term-numerator))
	(setq sum-exponent (- (fpz-exponent sum-numerator)
			      (fpz-exponent sum-denominator)))
	(setq max-sum-exponent (max max-sum-exponent sum-exponent)))
    (setq answer (fpz/ sum-numerator sum-denominator))
    ;; The accuracy estimate can be about three bits too big
    (setq bits-accuracy
	  (max 0 (- fpprec (- max-sum-exponent (fpz-exponent answer) -3))))
    (return (values answer bits-accuracy))))

;; Exponent of an FP number.  Return large -ve number if zero
(defun fp-exponent (x) 
  (if (equal (car x) 0) most-negative-fixnum (cadr x)))

;; Given a complex FP number, return larger exponent
(defun fpz-exponent (z)
  (max (fp-exponent (fpz-re z)) (fp-exponent (fpz-im z))))