7.145 Bugfix of the factorial function in Maple 7 (5.4.02)

7.145.1 Helmut Kahovec

The well known bug in the factorial function of Maple7 can be fixed as shown below. The local procedures '_factorial' and 'fact' in the original module definition are flawed and may be corrected as indicated. Note that in ’fact’ the difference N-n may also be negative.

> test:=module() 
  local _factorial,cacheN,cacheF,fact; 
    _factorial:=proc(n::algebraic) 
    local i; 
      if type(n,'integer') then 
        if n<0 then 
          NumericEvent( 
            'real_to_complex', 
            NumericEvent( 
              'division_by_zero', 
              infinity+infinity*I 
            ) 
          ) 
        elif n=0 or n=1 then 
          1 
        elif n<50 then 
          mul(i,i=2..n) 
        elif n=cacheN then 
          cacheF 
        elif n>cacheN then 
          cacheF:=cacheF*fact(1+cacheN,n); 
          cacheN:=n; 
          cacheF 
#       elif 49*cacheN<50*n then               <=== deleted === 
#         cacheF:=iquo(cacheF,fact(1+n,n));    <=== deleted === 
#         cacheN:=n;                           <=== deleted === 
#         cacheF                               <=== deleted === 
        else 
          cacheF:=fact(2,n); 
          cacheN:=n; 
          cacheF 
        end if 
      elif type(n,'complex(float)') then 
        if n=0. then 
          n+1.0 
        else 
          evalf(evalf[Digits+2](n*GAMMA(n))) 
        end if 
      else 
        ':-factorial(args)' 
      end if 
    end proc; 
    cacheN,cacheF:=1000,1000!; 
    fact:=proc(n,N) 
    local q; 
      if n=N then                   # <=== corrected === 
        n                           # <=== corrected === 
      elif N>n and N-n<=32 then     # <=== corrected === 
        mul(q,q=n..N)               # <=== corrected === 
      elif n>N and n-N<=32 then     # <=== corrected === 
        mul(q,q=N..n)               # <=== corrected === 
      else 
        q:=iquo(N+n,2); 
        fact(n,q)*fact(1+q,N) 
      end if 
    end proc 
  end module:
 

> unprotect('Factorial'); 
> Factorial:=op(3,eval(test))[1]; 
 
                       Factorial := _factorial 
 
> protect('Factorial'); 
> testlib:="./test/lib": 
> mkdir("./test"); 
> mkdir(testlib); 
> march('create',testlib,100); 
> savelibname:=testlib: 
> restart; 
> savelib('Factorial');
 

We now restart the Maple7 session and demonstrate the bug in the original factorial function.

> restart; 
> libname; 
 
                     "C:\\Programme\\Maple7/lib" 
 
> 1001!/1000!,1002!/1000!,1002!/1001!,1003!/1000!,1003!/1001!,1003!/1002!; 
 
                           1, 1, 1, 1, 1, 1
 

We now restart the Maple7 session again and appropriately extend ’libname’. The factorial function seems to give correct results now:

> restart; 
> libname:="./test/lib",libname: libname; 
 
              "./test/lib", "C:\\Programme\\Maple7/lib" 
 
> 1001!/1000!,1002!/1000!,1002!/1001!,1003!/1000!,1003!/1001!,1003!/1002!; 
 
            1001, 1003002, 1002, 1006011006, 1005006, 1003
 

It is corrected with Maple 8 (U. Klein)