help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: address@hidden
Subject: Re: Bug or numerical rounding/truncation errors in factorial.m function?
Date: Mon, 4 Oct 2010 02:44:41 -0700 (PDT)

For some reason, our Professor Gallagher suggests using, as a faster approximation of factorial(#), exp(gammaln(#+1))
Is this related to the accuracy issue?
octave-3.2.4.exe:4>factorial(18)
ans = 6402373705728000
octave-3.2.4.exe.5>exp(gammaln(19))
ans = 6402373705727994
Brian


From: Francesco Potortì <address@hidden>
To: José Luis García Pallero <address@hidden>
Cc: Octave bugs list <address@hidden>; Octave users list <address@hidden>
Sent: Mon, October 4, 2010 5:12:13 AM
Subject: Re: Bug or numerical rounding/truncation errors in factorial.m function?

>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/
_______________________________________________
Help-octave mailing list
address@hidden
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave


reply via email to

[Prev in Thread] Current Thread [Next in Thread]