I have implemented a Maxima function factorial_double(z).
This function is a generaliziation of genfact() for real and complex values.
Further I have tested this function with the parser (See the code below).
It works fine. The function genfact() is still present, but is no longer used by
the parser. The testsuite has no problems with this change. There are no tests
in the testsuite which depends on the return of %genfact as a noun form.
With this changes we would have a function genfact() which is specialized to
positive integer arguments and a new function factorial_double(z) generalized to
real and complex arguments. The operator !! would be translated to the noun form
of the function factorial_double.
At last it would be nice to extend the function factorial to complex arguments
too.
Furthermore some more code can be implemented to support the derivative, give
mirror symmetry to the function or do some transformations and argument
simplifications.
So what do you think?
Here is the code for the Double factorial:
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(def-mheader |$!!| (%factorial_double))
(def-led (|$!!| 160.) (op left)
(list '$expr
(mheader '$!!)
(convert left '$expr)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun $factorial_double (z)
(simplify (list '(%factorial_double) (resimplify z))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Set the properties noun, verb, alias and reversealias
(defprop $factorial_double %factorial_double alias)
(defprop $factorial_double %factorial_double verb)
(defprop %factorial_double $factorial_double reversealias)
(defprop %factorial_double $factorial_double noun)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defprop %factorial_double simp-factorial-double operators)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun simp-factorial-double (expr z simpflag)
(oneargcheck expr)
(setq z (simpcheck (cadr expr) simpflag))
(cond
((and (fixnump z) (> z -1) (or (minusp $factlim) (< z $factlim)))
;; Positive Integer less then $factlim or $factlim is -1. Call gfact.
(gfact z (floor (/ z 2)) 2))
((and (mnump z)
(eq ($sign z) '$neg)
(zerop1 (sub ($truncate (div z 2)) (div z 2))))
;; Even negative integer or real representation. Not defined.
(merror "factorial_double(~:M) is undefined." z))
((or (integerp z) ; at this point odd negative integer
(floatp z)
(complex-number-p z))
;; Odd negative integer, float or complex float.
(complexify
(factorial-double
(complex (float ($realpart z)) (float ($imagpart z))))))
((and (not (ratnump z))
(or ($bfloatp z)
(complex-number-p z 'bigfloat-or-number-p)))
;; bigfloat or complex bigfloat.
(bfloat-factorial-double z))
(t
(eqtest (list '(%factorial_double) z) expr))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Double factorial for a complex float argument. The result is a CL complex.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun factorial-double (z)
(let ((pival (float pi)))
(*
(expt
(/ 2 pival)
(/ (- 1 (cos (* pival z))) 4))
(expt 2 (/ z 2))
(gamma-lanczos (+ 1 (/ z 2))))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Double factorial for a bigfloat or complex bigfloat argument
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun bfloat-factorial-double (z)
(let ((pival ($bfloat '$%pi)))
($rectform
(mul
(power
(div 2.0 pival)
(div (sub 1 (simplify (list '(%cos) (mul pival z)))) 4.0))
(power 2.0 (div z 2.0))
(simplify (list '(%gamma) (add 1.0 (div z 2.0))))))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
Dieter Kaiser