Why is this factorial algorithm not exact

Question

Sorry, it seems to me stupid to ask about this, and I'm ready to lose half of my questions when asking about this, but why does this algorithm not work? It works to a certain extent. After number 13, factorials are a bit gone. For example, the numbers do not completely match hundreds of thousands of places on.

#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 conclusion:

 Factorial of 13 is 6227020800.000000 

Here is an example of inaccurate output:

 Factorial of 14 is 87178289152.000000 

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

14 87178291200

I changed the return type to float to get a more accurate result, but I got this code for the most part from here: https://www.tutorialspoint.com/cprogramming/c_recursion.htm

EDIT: If I switch to the return type to double the output, accurate to 21. I use the% Lf line formatter to output to the printf function.

Answer 1

Someone posted a similar question a while ago. The consensus was that if you are writing it for work, use a large library of numbers (e.g. GMP), and if this programming exercise writes out a solution using an array of characters.

For instance:

 /* fact50.c calculate a table of factorials from 0! to 50! by keeping a running sum of character digits */ #include <stdio.h> #include <string.h> int main (void) { printf ("\n Table of Factorials\n\n"); // length of arrays = 65 character digits char str[] = "00000000000000000000000000000000000000000000000000000000000000000"; char sum[] = "00000000000000000000000000000000000000000000000000000000000000001"; const int len = strlen (str); int index; for ( int i = 0; i <= 50; ++i ) { memcpy (str, sum, len); for ( int j = 1; j <= i - 1; ++j ) { index = len - 1; int carry = 0; do { int digit = (sum[index] - '0') + (str[index] - '0') + carry; carry = 0; if ( digit > 9 ) { carry = 1; digit %= 10; } sum[index] = digit + '0'; --index; } while ( index >= 0 ); } printf ("%2i! = ", i); for ( index = 0; sum[index] == '0'; ++index ) printf ("%c", '.'); for ( ; index < len; ++index ) printf ("%c", sum[index]); printf ("\n"); } return 0; } 

Answer 2

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

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

Answer 3

Why is this factorial algorithm not accurate

There is nothing wrong with your algorithm . Simply used data types have a limit on the maximum amount that they can store. This will be a problem, no matter which algorithm you choose. You can change data types from float to something like a long double to hold on to something more. But in the end, it will fail anyway as soon as the factor value exceeds the capacity of this data type. In my opinion, you should put the condition in your factorial function to return without calculating anything if the argument passed in the argument is greater than the value that your selected data type can support.

Answer 4

float can represent a wider range of numbers than int , but it cannot represent all the values ​​in this range - as you approach the edge of the range (i.e. as the values ​​increase), the gap between the values ​​presented becomes 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, the IEEE-754 with one precision float gives you a bit more accuracy).

Having said that the float should be able to represent all integer values ​​up to 2 ²⁴ so I'll bet that the problem is calling printf - the float parameters "advance" to double , so there is a change in representation, and this may take into account the lost accuracy.

Try changing the return type of factorial to double and see if that helps.

<gratis rant>

Every time I see a recursive factorial function, I want to scream. Recursion in this particular case does not improve code clarity or performance compared to an iterative solution:

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

and actually can lead to worse performance due to the overhead of so many function calls.

Yes, the definition of factorial is recursive, but the implementation of a factorial function is optional. The same goes for Fibonacci sequences. There is even a closed-loop solution for Fibonacci numbers

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

which does not require any cycle in the first place.

Recursion is great for algorithms that split their data into a relatively small number of subsets of the same size (Quicksort, tree traversal, etc.). Something like this, where the partition is N-1 subsets of 1 element? Not so much.

</ grantuit rant>

Answer 5

OP meets float precision. For a typical float , integer values ​​above 16777216.0f are not all exactly representable. Some factorial results above this point are accurately presented.

Let's try it with different types.
At 11! float results exceed 16777216.0f and exactly match.
At 14! float result is inaccurate due to limited accuracy.
At 23! double result is inaccurate due to limited accuracy.

At 21! answer exceeds my uintmax_t range.
At 35! The answer exceeds my float range.
In 171! The answer exceeds my double range.

The string representation is infinitely accurate until it reaches the limits in the buffer.

 #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...