Why Is This Factorial Algorithm Not Accurate

Question

Sorry I feel stupid asking this and am prepared to lose half of my points asking this but why does this algorithm not work? It works up to a point. After the number 13 the factorials are a little off. For instance the numbers do not entirely match in the hundreds thousands place and onward.

#include <stdio.h>

float factorial(unsigned int i) {

    if (i <= 1) {
        return 1;
    }
    return i * factorial(i - 1);
}

int  main() {
    int i = 13;
    printf("Factorial of %d is %f\n", i, factorial(i));
    return 0;
}

Here's the output:

Factorial of 13 is 6227020800.000000

Here is an example of inaccurate output:

Factorial of 14 is 87178289152.000000

The output for the number 14 should actually be this (from mathisfun.com)

14        87,178,291,200

I changed the return type to float to obtain more accurate output but I obtained this code for the most part from here: https://www.tutorialspoint.com/cprogramming/c_recursion.htm

EDIT: If I change to the return type to double the output is accurate up to 21.I am using the %Lf string formatter for the output in the printf function.

Answer 1

Simple. float cannot accurately store integers above 16777216 without loss of precision.

int is better than float. But try long long so you can properly store 19 digits.

Answer 2

float can represent a wider range of numbers than int, but it cannot represent all the values within that range - as you approach the edge of the range (i.e., as the magnitudes of the values increase), the gap between representable values gets wider.

For example, if you cannot represent values between 0.123 and 0.124, then you also cannot represent values between 123.0 and 124.0, or 1230.0 and 1240.0, or 12300.0 and 12400.0, etc. (of course, IEEE-754 single-precision float gives you a bit more precision than that).

Having said that, float should be able to represent all integer values up to 2²⁴ exactly, so I'm going to bet the issue is in the printf call - float parameters are "promoted" to double, so there's a representation change involved, and that may account for the lost precision.

Try changing the return type of factorial to double and see if that doesn't help.

<gratuitous rant>

Every time I see a recursive factorial function I want to scream. Recursion in this particular case offers no improvement in either code clarity or performance over an iterative solution:

double fac( int x )
{
  double result = 1.0;
  while ( x )
  {
    result *= x--;
  }
  return result;
}

and can in fact result in worse performance due to the overhead of so many function calls.

Yes, the definition of a factorial is recursive, but the implementation of a factorial function doesn't have to be. Same for Fibonacci sequences. There's even a closed form solution for Fibonacci numbers

Fn = ((1 + √5)n - (1 - √5)n) / (2n * √5)

that doesn't require any looping in the first place.

Recursion's great for algorithms that partition their data into relatively few, equal-sized subsets (Quicksort, tree traversals, etc.). For something like this, where the partitioning is N-1 subsets of 1 element? Not so much.

</gratuitous rant>

Answer 3

OP is encountering the precision limits of float. For typical float, whole number values above 16777216.0f are not all exactly representable. Some factorial results above this point are exactly representable.

Let us try this with different types.
At 11!, the float results exceeds 16777216.0f and is exactly correct.
At 14!, the float result is imprecise because of limited precision.
At 23!, the double result is imprecise because of limited precision.

At 22!, the answer exceeds my uintmax_t range. (64-bit)
At 35!, the answer exceeds my float range.
At 171!, the answer exceeds my double range.

A string representation is accurate endlessly until it reaches buffer limitations.

#include <stdint.h>
#include <string.h>
#include <stdio.h>

uintmax_t factorial_uintmax(unsigned int i) {
  if (i <= 1) {
    return 1;
  }
  return i * factorial_uintmax(i - 1);
}

float factorial_float(unsigned int i) {
  if (i <= 1) {
    return 1;
  }
  return i * factorial_float(i - 1);
}

double factorial_double(unsigned int i) {
  if (i <= 1) {
    return 1;
  }
  return i * factorial_double(i - 1);
}

char * string_mult(char *y, unsigned base, unsigned x) {
  size_t len = strlen(y);
  unsigned acc = 0;
  size_t i = len;
  while (i > 0) {
    i--;
    acc += (y[i] - '0') * x;
    y[i] = acc % base + '0';
    acc /= base;
  }
  while (acc) {
    memmove(&y[1], &y[0], ++len);
    y[0] = acc % base + '0';
    acc /= base;
  }
  return y;
}

char *factorial_string(char *dest, unsigned int i) {
  strcpy(dest, "1");
  for (unsigned m = 2; m <= i; m++) {
    string_mult(dest, 10, m);
  }
  return dest;
}

void factorial_test(unsigned int i) {
  uintmax_t u = factorial_uintmax(i);
  float f = factorial_float(i);
  double d = factorial_double(i);
  char s[2000];
  factorial_string(s, i);
  printf("factorial of %3d is uintmax_t: %ju\n", i, u);
  printf("                    float:     %.0f %s\n", f, "*" + (1.0 * f == u));
  printf("                    double:    %.0f %s\n", d, "*" + (d == u));
  printf("                    string:    %s\n", s);
}

int main(void) {
  for (unsigned i = 11; i < 172; i++)
    factorial_test(i);
  return 0;
}

Output

factorial of  11 is uintmax_t: 39916800
                    float:     39916800 
                    double:    39916800 
                    string:    39916800
factorial of  12 is uintmax_t: 479001600
                    float:     479001600 
                    double:    479001600 
                    string:    479001600
factorial of  13 is uintmax_t: 6227020800
                    float:     6227020800 
                    double:    6227020800 
                    string:    6227020800
factorial of  14 is uintmax_t: 87178291200
                    float:     87178289152 *
                    double:    87178291200 
                    string:    87178291200
factorial of  20 is uintmax_t: 2432902008176640000
                    float:     2432902023163674624 *
                    double:    2432902008176640000 
                    string:    2432902008176640000
factorial of  21 is uintmax_t: 14197454024290336768
                    float:     51090940837169725440 *
                    double:    51090942171709440000 *
                    string:    51090942171709440000
factorial of  22 is uintmax_t: 17196083355034583040
                    float:     1124000724806013026304 *
                    double:    1124000727777607680000 *
                    string:    1124000727777607680000
factorial of  23 is uintmax_t: 8128291617894825984
                    float:     25852017444594485559296 *
                    double:    25852016738884978212864 *
                    string:    25852016738884976640000
factorial of  34 is uintmax_t: 4926277576697053184
                    float:     295232822996533287161359432338880069632 *
                    double:    295232799039604119555149671006000381952 *
                    string:    295232799039604140847618609643520000000
factorial of  35 is uintmax_t: 6399018521010896896
                    float:     inf *
                    double:    10333147966386144222209170348167175077888 *
                    string:    10333147966386144929666651337523200000000
factorial of 170 is uintmax_t: 0
                    float:     inf *
                    double:    72574156153079940453996357155895914678961840000000... *
                    string:    72574156153079989673967282111292631147169916812964...
factorial of 171 is uintmax_t: 0
                    float:     inf *
                    double:    inf *
                    string:    12410180702176678234248405241031039926166055775016...

Answer 4

Why Is This Factorial Algorithm Not Accurate

There's nothing wrong in your algorithm as such. It is just that the data types you use have a limit for the highest number they can store. This will be a problem no matter which algorithm you choose. You can change the data types from float to something like long double to hold something bigger. But eventually it will still start failing once the factorial value exceeds the capacity of that data type. In my opinion, you should put an a condition in your factorial function to return without calculating anything if the passed in argument is greater than a value that your chosen datatype can support.

Answer 5

As explained in the other answers it has to do with the range and precision of the different types of numbers available.

It is possible to use a character array, a string, to represent the numbers and obtain exact values. For example:

/* Fact50.c -- display a table of factorials from 0! to 50!

   50! is 65 digits in length
*/

#include <stdio.h>
#include <string.h>

#define MAXNUM      50
#define LEN         65

int multiply(int n, char fact[], int limit)
{
    int carry = 0;
    for (int i = LEN - 1; i >= limit; --i) {
        int product = (fact[i] - '0') * n + carry;
        fact[i] = product % 10 + '0';
        carry = 0;
        if (product > 9) {
            carry = product / 10;
            if (i == limit) {
                --limit;
                fact[limit] = '0';
            }
        }
    }

    return limit;
}

int main(void)
{
    printf("\n"
           "                          Table of Factorials\n"
           "\n");

    char fact[LEN + 1] = { [LEN - 1] = '1', [LEN] = '\0' };
    memset(fact, '.', LEN - 1);

    printf("%2d! = %s\n", 0, fact);
    printf("%2d! = %s\n", 1, fact);

    int limit = LEN - 1;
    for (int i = 2; i <= MAXNUM; ++i) {
        limit = multiply(i, fact, limit);
        printf("%2d! = %s\n", i, fact);
    }

    return 0;
}

program output:

                           Table of Factorials

 0! = ................................................................1
 1! = ................................................................1
 2! = ................................................................2
 3! = ................................................................6
 4! = ...............................................................24
 5! = ..............................................................120
 6! = ..............................................................720
 7! = .............................................................5040
 8! = ............................................................40320
 9! = ...........................................................362880
10! = ..........................................................3628800
11! = .........................................................39916800
12! = ........................................................479001600
13! = .......................................................6227020800
14! = ......................................................87178291200
15! = ....................................................1307674368000
16! = ...................................................20922789888000
17! = ..................................................355687428096000
18! = .................................................6402373705728000
19! = ...............................................121645100408832000
20! = ..............................................2432902008176640000
21! = .............................................51090942171709440000
22! = ...........................................1124000727777607680000
23! = ..........................................25852016738884976640000
24! = .........................................620448401733239439360000
25! = .......................................15511210043330985984000000
26! = ......................................403291461126605635584000000
27! = ....................................10888869450418352160768000000
28! = ...................................304888344611713860501504000000
29! = ..................................8841761993739701954543616000000
30! = ................................265252859812191058636308480000000
31! = ...............................8222838654177922817725562880000000
32! = .............................263130836933693530167218012160000000
33! = ............................8683317618811886495518194401280000000
34! = ..........................295232799039604140847618609643520000000
35! = ........................10333147966386144929666651337523200000000
36! = .......................371993326789901217467999448150835200000000
37! = .....................13763753091226345046315979581580902400000000
38! = ....................523022617466601111760007224100074291200000000
39! = ..................20397882081197443358640281739902897356800000000
40! = .................815915283247897734345611269596115894272000000000
41! = ...............33452526613163807108170062053440751665152000000000
42! = .............1405006117752879898543142606244511569936384000000000
43! = ............60415263063373835637355132068513997507264512000000000
44! = ..........2658271574788448768043625811014615890319638528000000000
45! = ........119622220865480194561963161495657715064383733760000000000
46! = .......5502622159812088949850305428800254892961651752960000000000
47! = .....258623241511168180642964355153611979969197632389120000000000
48! = ...12413915592536072670862289047373375038521486354677760000000000
49! = ..608281864034267560872252163321295376887552831379210240000000000
50! = 30414093201713378043612608166064768844377641568960512000000000000