After an evaluation, GNOME has moved from Bugzilla to GitLab. Learn more about GitLab.
No new issues can be reported in GNOME Bugzilla anymore.
To report an issue in a GNOME project, go to GNOME GitLab.
Do not go to GNOME Gitlab for: Bluefish, Doxygen, GnuCash, GStreamer, java-gnome, LDTP, NetworkManager, Tomboy.
Bug 410005 - FACT(17) is odd
FACT(17) is odd
Status: RESOLVED FIXED
Product: Gnumeric
Classification: Applications
Component: General
git master
Other All
: High normal
: ---
Assigned To: Jody Goldberg
Jody Goldberg
Depends on:
Blocks:
 
 
Reported: 2007-02-20 13:29 UTC by Mick Twohig
Modified: 2007-02-20 18:17 UTC
See Also:
GNOME target: ---
GNOME version: ---



Description Mick Twohig 2007-02-20 13:29:07 UTC
Evaluation of 17! gave value of 355687428096001 which is out by one. Similar errors were observed for factorials of 18, 19, 21 to 37.

For all n>= 5, n! is divisible by 10. Modulo 10 division of each factorial in the worksheet was expected to yield 0, but was found to be non-zero for the ranges given.

17! was the last evaluation with a fully expanded result, from 18! onwards, exponential notation is used.

The machine this was tested on was a Pentium 3. Similar calculations at the time in OOCalc gave correct results.
Comment 1 Martin Burton 2007-02-20 17:27:11 UTC
Aha, a little poking around in the gnumeric source code reveals the problem:

static GnmValue *gnumeric_fact (FunctionEvalInfo *ei, GnmValue const * const *argv)
{
        gnm_float x = value_get_as_float (argv[0]);
        gboolean x_is_integer = (x == gnm_floor (x));

        if (x < 0 && x_is_integer)
                return value_new_error_NUM (ei->pos);
        if (x > 12 || !x_is_integer) {
                int sign;
                gnm_float tmp = gnm_lgamma_r (x + 1, &sign);
                gnm_float res = sign * gnm_exp (tmp);
                if (x_is_integer)
                        res = gnm_floor (res + 0.5);  /* Round, just in case. */
                return value_new_float (res);
        } else
                return value_new_int (fact (x));
}

which shows that gnumeric  uses the gamma/exponent method of calculating the factorial (for any factorial above 12).  Which only returns an approximation.  It then tries to fudge the approximation by rounding, which actually just makes things worse!

The gnm_* math functions are mapped through #DEFINEs to the standard math library functions so I wrote a little test program  that shows what's going on:

#include <stdio.h>
#include <math.h>

int main(void) {
        int x = 17;
        int sign;
        double tmp;
        double res;

        tmp = lgamma_r(x+1, &sign);
        res = sign * exp(tmp);

        printf("Factorial (gamma function calculated) of %d = %f\n",x,res);

        res = floor(res + 0.5);

        printf("and floored to: %f\n",res);

}

Which results in:

martin@apollo:~/src/math$ ./fact 
Factorial (gamma function calculated) of 17 = 355687428096000.687500
and floored to: 355687428096001.000000
martin@apollo:~/src/math$


Guess a more accurate factorial calculation is required! :-)
Comment 2 Morten Welinder 2007-02-20 18:03:21 UTC
Good catch.

I think what we should do is to use a table up to the point where the
numbers get so large that they cannot be represented exactly.  At that
point, rounding becomes irrelevant.
Comment 3 Morten Welinder 2007-02-20 18:10:09 UTC
17! is off by -1
18! is off by +6
19! is off by -80.

Which just tells us what we knew: using lgamma is suspect since exp can't
help blow up rounding errors.  I, naively, thought the rounding saved us.
(And it probably does <17.)

Comment 4 Morten Welinder 2007-02-20 18:17:14 UTC
This problem has been fixed in the development version. The fix will be available in the next major software release. Thank you for your bug report.