Math challenge: Evaluate big factorials like 1000!

[jack]

how can you say it's not broke if it does not work?
clearly _Float claims to be 80-bit float which it isn't, and if it isn't then why bother having it if it's the same as double?
I know MS C treats long double the same as double, but what does that have to do with _float?
btw, this same behavior is observed with QB64 under linux, now I know linux supports 80-bit long double
and if you had taken the trouble to look at my code you would have noticed that the calculations were done in quasi 80-bit float, otherwise why would it give a correct answer?
one thing is certain, it's stupid for me to waste my time on this forum and with QB64
from the QB64 Wiki

QB64 always allocates 32 bytes to store this value.
It is safe to assume this value is at least as precise as DOUBLE.
Under the current implementation it is stored in a 10-byte floating point variable.
_FLOAT variables can also use the ## variable name type suffix.
Values returned may be expressed using exponential or scientific notation using E for SINGLE or D for DOUBLE precision.
According to IEEE-754 this can store a value of up to 1.1897E+4932 compared to a DOUBLE which goes up to 1.7976E+308.

[david_uwi]

This was a challenge on the qbasic forum several years (decades?)ago. Here is my solution (you may need to make the window bigger if you want to do really big numbers (if the number of digits is > 32768 change mx).

Code: QB64: [Select]

1. DEFINT A-S
2. mx = 32000
3. WIDTH 140
4. DIM f(mx) AS INTEGER
5. INPUT "input number"; n
6. f(mx) = 1
7. m = mx
8. FOR i = 2 TO n
9.     FOR j = mx TO m STEP -1
10.         f(j) = f(j) * i + ic
11.         ic = f(j) \ 10
12.         f(j) = f(j) MOD 10
13.     NEXT j
14.     WHILE ic <> 0
15.         m = m - 1
16.         f(m) = f(m) + ic
17.         ic = f(m) \ 10
18.         f(m) = f(m) MOD 10
19.     WEND
20. NEXT i
21. FOR j = m TO mx
22.     PRINT USING "#"; f(j);
23. NEXT j
24. PRINT " ="; n; "!"
25.  

I've got another solution using an approximation N! approx= N ln(N)-N +1/2 ln(2*pi*N)

Code: QB64: [Select]

1. DEFDBL A-Z
2. n = 50
3. pi = 3.141593
4. a = n * LOG(n) - n + 0.5 * LOG(2 * pi * n)
5. b = a / LOG(10)
6. c = b - INT(b)
7. PRINT 10 ^ c; "x10^"; INT(b)

It's nowhere near as accurate.

[bplus]

Oh man that is nice, that first one stays accurate doesn't it?

[david_uwi]

Thanks it is kinda neat. It is totally accurate and may get a little slow for VERY big numbers. Notice how many zeros there are at the end due to multiplying by multiples of 10 (and 2&5s) and once you get a zero at the end it stays.

[bplus]

Well I may now have to take another look at string math, after seeing this and Petr's recent post.

[SMcNeill]

how can you say it's not broke if it does not work?
clearly _Float claims to be 80-bit float which it isn't, and if it isn't then why bother having it if it's the same as double?
I know MS C treats long double the same as double, but what does that have to do with _float?
btw, this same behavior is observed with QB64 under linux, now I know linux supports 80-bit long double
and if you had taken the trouble to look at my code you would have noticed that the calculations were done in quasi 80-bit float, otherwise why would it give a correct answer?
one thing is certain, it's stupid for me to waste my time on this forum and with QB64
from the QB64 Wiki

Except it *does* work.  I'll share you a small program and a screenshot of it working:

Code: QB64: [Select]

1. f## = 8## / 9##
2. PRINT "  123456789012345678901"
3. PRINT USING "#.#####################"; f##
4. PRINT f##
5.  
6. f# = 8## / 9##
7.  
8. PRINT USING "#.#####################"; f#
9. PRINT f#

Now, as you can tell by looking at the screenshot, PRINT is only printing 16 digits of the value (as I've mentioned elsewhere before).  When we format our output by using PRINT USING, you can compare and see that DOUBLE loses precision after the 16th decimal place, whereas FLOAT holds the precision to the 19th decimal place. 

As you can also see from the screenshot, this is QB64x64 version which I'm running the program on, so obviously it's not an issue inside QB64x64 itself.

So why are you seeing different results with your code?? 

It's as I told you before:  FLOAT translates into a long double variable type in C, and how that long double operates depends on various factors such as your PC's CPU/math coprocessor, OS, compiler, and flags set. 

In C and related programming languages, long double refers to a floating-point data type that is often more precise than double-precision though the language standard only requires it to be at least as precise as double.

Specifications say an apple pie must be a pie that was baked with apple filling.  Some places offer apple pie with whip cream as their standard apple pie...   The cook in your PC is providing you an apple pie and you're exclaiming, "This restaurant is broken!  There's no whip cream on my apple pie!!"

There doesn't have to be whip cream on an apple pie for it to be an apple pie, and _FLOATS are only required to be as precise as a double.
 Anything beyond that is extra.

As for the wiki reference, you need to remember that a majority of it was written back in the SDL days of QB64 when there was only a 32-bit version of QB64 available.  It's only recently that a 64-bit version was packaged and offered officially for the user base, and not everything inside the wiki may cover all the minute differences in the two versions.  A note may need to be made somewhere in the wiki about how precision can be affected by how long doubles are handled in 64-bit versions, but there's nothing "broken" which can be fixed.  At least, at least not on my machine there isn't -- as you can see from the screenshot I provide below.

All I can say, since you appear to not be satisfied with the results you're achieving, is for you to look into the source itself and see if you can correct the issue as it occurs for your system.  The fix for you may be as simple as just using a different set of compiler options with your system, and you can find a list of those switches here: https://linux.die.net/man/1/g++

Some of the flags which seem like they may be relevant to your issue, seems to me to possibly be:

-mabi=ibmlongdouble: Change the current ABI to use IBM extended precision long double. This is a PowerPC 32-bit SYSV ABI option.

-mabi=ieeelongdouble: Change the current ABI to use IEEE extended precision long double. This is a PowerPC 32-bit Linux ABI option.

-mlong-double-64//-mlong-double-128: These switches control the size of "long double" type. A size of 64bit makes the "long double" type equivalent to the "double" type. This is the default.

-mhard-quad-float: Generate output containing quad-word (long double) floating point instructions.

-msoft-quad-float: Generate output containing library calls for quad-word (long double) floating point instructions.

-mno-align-double: Control whether GCC aligns "double", "long double", and "long long" variables on a two word boundary or a one word boundary. Aligning "double" variables on a two word boundary will produce code that runs somewhat faster on a Pentium at the expense of more memory.
On x86-64, -malign-double is enabled by default.

-mfpmath=387: Use the standard 387 floating point coprocessor present majority of chips and emulated otherwise. Code compiled with this option will run almost everywhere. The temporary results are computed in 80bit precision instead of precision specified by the type resulting in slightly different results compared to most of other chips.

[STxAxTIC]

Glad you found us david! I remember code similar to this from *.net in 2014, maybe sometime in June. Nice and fast!

[david_uwi]

My file says July 2014. I didn't realize it was so recent after all we lost Mac in 2008 - what a great programmer he was.