[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Bug or numerical rounding/truncation errors in factorial.m function?
From: |
José Luis García Pallero |
Subject: |
Bug or numerical rounding/truncation errors in factorial.m function? |
Date: |
Sun, 3 Oct 2010 13:35:37 +0200 |
Hello,
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
factorial(18)
ans = 6402373705727994 *****************different
prod(1:18)
ans = 6402373705728000 *****************different
factorial(19)
ans = 121645100408832080 *****************different
prod(1:19)
ans = 121645100408832000 *****************different
factorial(20)
ans = 2432902008176640000
prod(1:20)
ans = 2432902008176640000
17, 18 and 19 are sufficient small for obtain the correct result via
direct multiplication. In the code of factorial.m function the
computation is performed as
x = round (gamma (n+1));
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.
--
*****************************************
José Luis García Pallero
address@hidden
(o<
/ / \
V_/_
Use Debian GNU/Linux and enjoy!
*****************************************
- Bug or numerical rounding/truncation errors in factorial.m function?,
José Luis García Pallero <=