Re: Bug or numerical rounding/truncation errors in factorial.m function?

[Francesco Potortì]

>I'm using Octave 3.2.4 in a Debian Sid box. I'm trying to compute some
>factorials and here is the results for some values:
>
>factorial(16)
>ans =  20922789888000
>prod(1:16)
>ans =  20922789888000
>
>factorial(17)
>ans =  355687428096001 *****************different
>prod(1:17)
>ans =  355687428096000 *****************different

[...]

>Is possible that gamma function has a bug or is only a problem with
>ieee arithmetic? I suggest for the computation of factorials for small
>numbers a behavior similar to the used in GSL, in the file
>specfun/gamma.c (sorry, I can't paste a link because savanna don't
>permit surfing the src). In this file, the factorials for numbers from
>0 to 170 are stored explicitly (hard coded) for speed.

Yes, in fact the function gamma could be more precise than that for
arguments 17 and 18, which can be represented exactly in a 53-bit
mantissa:

octave> factorial(18)
ans =  6402373705727994
octave> prod(1:18)
ans =  6402373705728000
octave> factorial(18)/bitmax
ans =  0.710806270035391

-- 
Francesco Potortì (ricercatore)        Voice: +39 050 315 3058 (op.2111)
ISTI - Area della ricerca CNR          Fax:   +39 050 315 2040
via G. Moruzzi 1, I-56124 Pisa         Email: address@hidden
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/