#4668 numeric evaluation of genfact

[Barton Willis]

I think that for numeric values, Maxima defines genfactexactly as

    mygenfact(x,y,z) := block([n : floor(y)], product(x -k*z,k,0,n-1))$

The user documentation doesn't say that floor is applied to the second argument, but the source code is clear about that. But Maxima refuses numerical evaluation for many cases:

(%i33)  mygenfact(3.9, 4, 2);
(%o33)  1.5561000000000011

(%i34)  genfact(3.9,4,2);
(%o34)  genfact(3.9,4,2)

Here is a proposal for numeric evaluation of genfact. I'll need to double check I got the extension to negative n correct:

(in-package :bigfloat)

(defun genfact-numeric (x n z)
  (cond
    ((>= n 0)
     (let ((acc 1))
       (dotimes (k n acc)
         (setf acc (* acc (- x (* k z)))))))
    (t
     ;; For n < 0, use genfact(x,n,z) * genfact(x-n*z,-n, z) = 1
     (let* ((m (- n))  ; m = -n > 0
            (den (genfact-numeric (+ x (* m z)) m z)))
       (when (zerop den) ; works on all bigfloat numbers
         (maxima::merror "genfact(~M, ~M, ~M) is not defined" x n z))
       (/ 1 den)))))

(in-package :maxima)

[Robert Dodier]

For the record, the source code for genfact in src/asum.lisp shows that it simplifies only for first and last arguments which are fixnums, and second argument such that floor yields a fixnum. Probably the documentation should mention that.

I didn't find any description of what genfact is supposed to implement, exactly. I don't see any statements in the code and a web search for "generalized factorial" doesn't seem to find anything. I think maybe it would be good if we can say what it is that we believe we are implementing.

I noticed an existing function double_factorial which seems to be for float and bigfloat arguments. It is a simplifying function. Looks like "!!" does not dispatch to double_factorial for float or bigfloat arguments, but rather returns a genfact expression. Maybe for float and bigfloat arguments, "!!" should punt to double_factorial?

[Barton Willis]

User documentation:

Function: genfact (x, y, z)
Returns the generalized factorial, defined as x (x-z) (x - 2 z) ... (x - (y - 1) z). Thus, when x is an integer, genfact (x, x, 1) = x! and genfact (x, x/2, 2) = x!!.

From the source code, it's clear the intent is for the second argument to be rounded down to the next integer (floor)--otherwise we'd have a product from 0 thru a half integer for the !! case.

I'll look at double_factorial.