Error using FACT()

[CV]

Hi

I tried to use fact() - factorial of a number and see that this function is limited to 21 max (from CT contrib: ctmath2.c).

But if you try fact(21), the result is:
51.090.942.171.709.430.000 xharbour + pelles C -> WRONG
51.090.942.171.709.400.000 EXCEL 2013! -> WRONG
51.090.942.171.709.440.000 harbour + mingw -> RIGHT
51.090.942.171.709.440.000 calculator from W10 -> RIGHT

Anyway, as it is limited to 21, does anyone know a way to calculate factorial() without any limit or at least use a greater limit like 30 (without rounding issues)?

I know: the problem relies on the double precision math that harbour/C uses.

I tried myself with my own function (with a for-loop or recursively) and obviously I get the same wrong result.
If I try with greater numbers than 21 I get rounding problems (the proper upper limit should be 20 for the CT source file, if someone can fix it).

A must will be to have the gamma function, but this is a totally different story.

This program was compiled with xharbour and the C compiler used was "PELLES C" from xharbour.com, don't know if it is present whith other C compiler.

I posted wrongly this message in the harbour forum, you can inspect here:
https://groups.google.com/d/msg/harbour-users/FAtI__V6ssI/07aPi0s5BAAJ

Regards
Claudio Voskian

[Enrico Maria Giordano]

Il 28/05/2020 14:58, CV ha scritto:

> Hi
>
> I tried to use fact() - factorial of a number and see that this function is limited to 21 max (from CT contrib: ctmath2.c).
>
> But if you try fact(21), the result is:
> 51.090.942.171.709.430.000 xharbour + pelles C -> WRONG
> 51.090.942.171.709.400.000 EXCEL 2013! -> WRONG
> 51.090.942.171.709.440.000 harbour + mingw -> RIGHT
> 51.090.942.171.709.440.000 calculator from W10 -> RIGHT

51.090.942.171.709.440.000 xharbour + BCC -> RIGHT

EMG

http://www.emagsoftware.it
http://www.emagsoftware.it/emgmusic
http://www.emagsoftware.it/spectrum
http://www.emagsoftware.it/tbosg

[dlzc]

Dear CV:

On Thursday, May 28, 2020 at 5:58:13 AM UTC-7, CV wrote:
...

> I tried to use fact() - factorial of a number and
> see that this function is limited to 21 max
> (from CT contrib: ctmath2.c).

WHY use a function, when you only need to store an array of 30 "numbers"?

...

> Anyway, as it is limited to 21, does anyone know
> a way to calculate factorial() without any limit
> or at least use a greater limit like 30 (without
> rounding issues)?

http://mathforum.org/library/drmath/sets/select/dm_factorial_list.html
... just treat them as strings, if required.

30! is 33 digits long, so handle the lowest order 18 digits as one number, and the highest order digits as another number.

> I know: the problem relies on the double
> precision math that harbour/C uses.

Break it into parts.

> I tried myself with my own function (with a
> for-loop or recursively) and obviously I get
> the same wrong result. If I try with greater
> numbers than 21 I get rounding problems (the
> proper upper limit should be 20 for the CT
> source file, if someone can fix it).

I once had a "rational" data type (in Quick Basic) with numerators and denominators for a given number, when solving arrays. The multiplication and division operations were then easy, but addition and subtraction took a lot of time. Would be similar overhead here.

I won't touch the gamma function. There are higher precision C libraries that have all the needed precision, but they'd have to be "wrapped" for use in (x)Harbour.

David A. Smith

[CV]

Hi Enrico

Thank you for the test with BCC.

So it seems to be an issue with the Pelles C compiler that xharbour.com uses from long time ago and never updated.

I sent to my friend the executable compiled with harbour+mingw and everything is fine now.

Regards
Claudio Voskian

[CV]

David

I can modify everything and write a dictionary with numbers, but the real issue is the error in an outdated C compiler (pelles), that is my main concern (or at least one that upsets me).

Regards
Claudio Voskian

[dlzc]

Dear CV:

On Thursday, May 28, 2020 at 8:36:25 AM UTC-7, CV wrote:
> David
>
> ... the real issue is the error in an outdated

> C compiler (pelles), that is my main concern
> (or at least one that upsets me).

http://www.smorgasbordet.com/pellesc/
... it appears not to have been updated since 2018.

It is not the basis for the commercial xHarbour, and is just one option that Mel links to on his download site. Might need to involve the developers at PellesC about this issue. Not something we can fix here.

David A. Smith

[meds...@gmail.com]

Hi David and Claudio:

Below is my FACTOR.PRG proggie compiled in both PellesC 9.0 and BCC 7.4, and the results of running both

Hope this helps.

-Mel

*********************************
Results from BCC 7.4 below
Factorial for nFac = 1 is: 1
Factorial for nFac = 2 is: 2
Factorial for nFac = 3 is: 6
Factorial for nFac = 4 is: 24
Factorial for nFac = 5 is: 120
Factorial for nFac = 6 is: 720
Factorial for nFac = 7 is: 5,040
Factorial for nFac = 8 is: 40,320
Factorial for nFac = 9 is: 362,880
Factorial for nFac = 10 is: 3,628,800
Factorial for nFac = 11 is: 39,916,800
Factorial for nFac = 12 is: 479,001,600
Factorial for nFac = 13 is: 6,227,020,800
Factorial for nFac = 14 is: 87,178,291,200
Factorial for nFac = 15 is: 1,307,674,368,000
Factorial for nFac = 16 is: 20,922,789,888,000
Factorial for nFac = 17 is: 355,687,428,096,000
Factorial for nFac = 18 is: 6,402,373,705,728,000
Factorial for nFac = 19 is: 121,645,100,408,832,000
Factorial for nFac = 20 is: 2,432,902,008,176,640,000
Factorial for nFac = 21 is: 51,090,942,171,709,440,000
Factorial for nFac = 22 is: -1
Factorial for nFac = 23 is: -1
Factorial for nFac = 24 is: -1
Factorial for nFac = 25 is: -1

Results from Pellesc 9.0 below

Factorial for nFac = 1 is: 1
Factorial for nFac = 2 is: 2
Factorial for nFac = 3 is: 6
Factorial for nFac = 4 is: 24
Factorial for nFac = 5 is: 120
Factorial for nFac = 6 is: 720
Factorial for nFac = 7 is: 5,040
Factorial for nFac = 8 is: 40,320
Factorial for nFac = 9 is: 362,880
Factorial for nFac = 10 is: 3,628,800
Factorial for nFac = 11 is: 39,916,800
Factorial for nFac = 12 is: 479,001,600
Factorial for nFac = 13 is: 6,227,020,800
Factorial for nFac = 14 is: 87,178,291,200
Factorial for nFac = 15 is: 1,307,674,368,000
Factorial for nFac = 16 is: 20,922,789,888,000
Factorial for nFac = 17 is: 355,687,428,096,000
Factorial for nFac = 18 is: 6,402,373,705,728,000
Factorial for nFac = 19 is: 121,645,100,408,832,000
Factorial for nFac = 20 is: 2,432,902,008,176,640,000
Factorial for nFac = 21 is: 51,090,942,171,709,430,000
Factorial for nFac = 22 is: -1
Factorial for nFac = 23 is: -1
Factorial for nFac = 24 is: -1
Factorial for nFac = 25 is: -1

The proggie that creatd the results above:

#define CRLF chr(13)+chr(10)

function Main()
// this program is built with PellesC64 and BCC 7.40 under xHarbour
local I, nFac
local cFactFile := "c:\cgi\FACTFILE.TXT"
local cFact,cFacts
cls
? "Computing Factorial for ever increasing numbers ..."

cFact := cFacts := ""

FOR I = 1 to 25
nFac := FACT(I)
cFact := " Factorial for nFac = "+TRANSFORM(I,"99")+" is: "+transform(nfac,"@Z 999,999,999,999,999,999,999,999")
? cFact
cFacts := cFacts + CRLF + cFACT
NEXT
cFacts := cFacts + CRLF

memowrit(cFactFile,cFacts)

? "Press Any Key to Finish."

inkey(0)

return NIL

***********************

[Ella Stern]

What about applying an algorithm like this?

public class MyTest
{
public static void main(String[] args)
{
double myResult = myFactorial(30.0);
System.out.println(myResult);
}

public static double myFactorial(double nr) {
double tmp = 1.0;
for (int i = 1; i <= nr; i++) {
tmp *= i * 0.1;
}
tmp *= Math.pow(10, nr);
return tmp;

[CV]

Hi MEL!

Thank you for the tests!!
It seems that PELLES C V9 is also wrong
No matter the version, it is the same bad result from the outdated Pelles C that xharbour.com use in their distro.

Users of xharbour.com have been using this Pelles C version unknown from the very beggining (In my case, from 2007), so I can't imagine how many operations like these have been done WRONGLY!

I beg someone in charge at xh.COM fixes it or use a different C compiler (I asked them many times without any luck; it seems no one cares about).

Regards
Claudio Voskian

[dlzc]

On Thursday, June 4, 2020 at 1:22:47 AM UTC-7, Ella Stern wrote:
> What about applying an algorithm like this?
>
> public class MyTest
> {
> public static void main(String[] args)
> {
> double myResult = myFactorial(30.0);
> System.out.println(myResult);
> }
>
> public static double myFactorial(double nr) {
> double tmp = 1.0;
> for (int i = 1; i <= nr; i++) {
> tmp *= i * 0.1;
> }
> tmp *= Math.pow(10, nr);
> return tmp;
> }
>
> }

He said in his first post that he tried writing his own loop, and got exactly the same result.

"I tried myself with my own function (with a for-loop or recursively) and obviously I get the same wrong result."

... so it is a "roundoff error" at extreme range of 64 bit integer math, unique to Pelles C.

Still not sure what he'd want a gamma() function for! ;-)

David A. Smith

[Ella Stern]

When using "double" instead of "long" there is possible to store huge numbers, but the precision is of 15 significant digits.

I'm not sure that "standard C" does have solution for getting a bigger precision. Java has Math.BigDecimal (in their Math library), C# has Decimal type (huge numbers stored as strings and undergoing to custom operations).

IMHO writing custom functions designed and tailored to the specific calculations is the way to go...

[meds...@gmail.com]

Hi Ella:

But the current problem is that Claudio Voskian is annoyed with PellesC because it is not the same as other C compilers such that it cannot do the following multiplication correctly:

21 * 2432902008176640000

This calc is performed by hand correctly, and correctly by BCC 7.4, MinGw 10.1, etc, etc.

But, Pellesc and PellesC64 both give a faulty result.

I have posted a Bug Report on the Pelles Forum, and am awaiting a response from Pelles.

-Mel

[meds...@gmail.com]

Hi All:

Pelle Orinius has responded to my bug report and has requested others to assist in helping too.

In the meantime he has provided a slightly better factorial function for others to test. He said he tried it using MSVC and got 'better results'.

It is displayed below.

-Mel


***** a newer Factorial function from Pelle *****
#include <stdio.h>

double tail_recursive_factorial_aux(double n, double acc) {
return n < 1 ? acc : tail_recursive_factorial_aux(n - 1, acc * n);
}
double tail_recursive_factorial(double n) {
return tail_recursive_factorial_aux(n, 1);
}

int main(void) {
for (double n = 1; n < 25; n++) {
union {
double d;
unsigned char b[8];
} u;
u.d = tail_recursive_factorial(n);
printf("%3.0f! = %40.15f - 0x%02hhX:0x%02hhX:0x%02hhX:0x%02hhX:0x%02hhX:0x%02hhX:0x%02hhX:0x%02hhX\n", n, u.d, u.b[7], u.b[6], u.b[5], u.b[4], u.b[3], u.b[2], u.b[1], u.b[0]);
}
}
***** End of factorial function *****

[CV]

Mel and everyone else

The fact() function is anecdotic, it was a mean to show a bug. A big one.

The main problem is the arithmetic algorithm used inside Pelles C (as it was not fixed in a decade or so): I can't imagine how many operations have been done erroneously.

It is enough for me *to not use it anymore* until fixed.

Regards
Claudio Voskian

[Drako]

El jueves, 28 de mayo de 2020, 7:58:13 (UTC-5), CV escribió:
> Hi
>
> I tried to use fact() - factorial of a number and see that this function is limited to 21 max (from CT contrib: ctmath2.c).
>

> Regards
> Claudio Voskian

I think the solution would be to implement The Poor Man's algorithm:
http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
Good enough to compute the factorial up to n=10000 in a few seconds.

[dlzc]

Dear Drako:

On Friday, June 5, 2020 at 12:01:52 PM UTC-7, Drako wrote:
...

> I think the solution would be to implement
> The Poor Man's algorithm:
> http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
> Good enough to compute the factorial up to
> n=10000 in a few seconds.

Sorry, no. The PellesC library will also break this. And the specific problem occurs when 21! is calculated, as described. Going into floats will be even less accurate.

[Darkwood]

On Thursday, May 28, 2020 at 5:58:13 AM UTC-7, CV wrote:

Try using alternative facts... :)

[Drako]

I thought that the problem was to solve the Factorial calculus, the Poor's man algorithm does that with light math using strings so it can go beyond integer math capabilities of FPU.

[dlzc]

Dear Drako:

What is shown on *that* link, involves direct multiplication. Which break Pelles C at 21!, and 64-bit integers (any compiler) above 21!. I did not look any further.

You probably meant to have someone click the link to this:
http://www.luschny.de/math/factorial/csharp/FactorialPoorMans.cs.html

But Claudio was mostly trying to highlight that the math functions behind Pelles C, which is used in xHarbour.com commercial products, is suspect. The fact() function was just a obvious / simple tip of an iceberg of indeterminate size.

David A. Smith

[CV]

Hi

To make things worst (the iceberg emerges):

x := FACT(20)
? x // 2432902008176640000.00 = Correct!
x := INT(x)
? x // 2432902008176640512 = INCORRECT!

Here the INT() function is working incorrectly and fact() is part of xharbour libs.
Another reason for not using this compiler anymore.

Using xharbour.com with pelles C.

Regards
Claudio Voskian