Text Only | Text with Attachments

QB64.org Forum

Active Forums => Programs => Topic started by: STxAxTIC on October 26, 2019, 10:40:41 pm

Title: Math challenge: Evaluate big factorials like 1000!
Post by: STxAxTIC on October 26, 2019, 10:40:41 pm
I think most people know that the factorial of a number (strictly integer for now) is that number multiplied by every integer less than it, all the way down to 1.

N! = N * (N-1) * (N-2) * ... * (N-(N-1))

so for example:

10! = 10*9*8*7*6*5*4*3*2*1

etc.

Whoever can write the tightest code to evaluate 100! and 1000! gets the prize of their answer being highlighted green. The answer should be in regular power-of-ten scientific notation with reasonable precision, not infinite. Note that using Stirling's approximation is in-exact and is not allowed. (See the screenshot.) And once again, you need not reach for big-number libraries.
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: Ashish on October 27, 2019, 01:37:39 am
I think this is it -
Code: QB64: [Select]
  1. 'Code to evaluate 100! or 1000! in scientific notation
  2. 'By Ashish for STxAxTIC
  3. FOR i = 1 TO 1000
  4.     v# = v# + (LOG(i) / LOG(10.#))
  5.     IF i<101 THEN k# = k# + (LOG(i)/LOG(10.#))
  7. PRINT "1000! = "; 10 ^ (v# - INT(v#)); " * 10^"; INT(v#)
  8. PRINT "100! = "; 10 ^ (k# - INT(k#)); " * 10^"; INT(k#)
  9.  

Is I am the winner?
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: STxAxTIC on October 27, 2019, 05:57:55 am
This is why everyone here loves you Ashish, perfect job.

For anyone else who hasn't mastered the question, the proof of why his answer works is attached below (it requires the results of the previous challenge, in theory.)
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: jack on October 27, 2019, 07:41:33 am
in this particular case you could have dispensed with logarithms and used _FLOAT
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: STxAxTIC on October 27, 2019, 10:05:38 am
Lol, fine. Pretend the question was 10!^(10!)... Nice observation though jack... you mean we could have just done brutal multiplication right? I dont know the limits of float off hand.
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: jack on October 27, 2019, 10:18:32 am
sorry that I spoiled the fun
the _FLOAT range should be 3.4E-4932 to 1.1E+4932
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: STxAxTIC on October 27, 2019, 10:24:11 am
Actually float starts failing somewhere before 175-factorial, it'll never make to 1000 unless this code is too tight:

Code: QB64: [Select]
  2. a = 1
  3. FOR k = 1 TO 165
  4.     a = a * k

Hat remains tipped in Ashish's direction.
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: jack on October 27, 2019, 10:57:51 am
@STxAxTIC
you are right, somehow _FLOAT is not an extended 80-bit float anymore or at least don't behave that way
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: jack on October 27, 2019, 11:10:51 am
here _FLOAT gives the extended result
Code: [Select]
f## = 1## / 3##
f## = (f## * 1000000000000000##) - 333333333333333##
PRINT f##

f# = 1# / 3#
f# = (f# * 1000000000000000#) - 333333333333333#
PRINT f#
but obviously _FLOAT is broken
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: STxAxTIC on October 27, 2019, 11:17:21 am
I think I know what a major issue is, which may or may not have led to the observed failure using _FLOAT.

Suppose we're halfway through a factorial calculation, and the intermediate result happens to contain more digits than the precision of whatever type we're using. For instance the temporary answer might be: 1.2345678XXXXXXXX * 10^50

...But we know the digits go beyond that final 8 but have been chopped off. Note though (and this is the important part) that there's no way to multiply a less precise number, even something as innocent as 3, into an an answer that's already imprecise at the single-digit order. Depending on what's after the last digit, you will even start losing track of the order of magnitude. For factorials this is effect especially bad - the calculation goes wrong early and just compounds wildly from there.

Anyway, thanks for turning this stone over.
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: SMcNeill on October 27, 2019, 11:33:11 am
You could probably multiply in groups and not have such a precision issue.

1000! = (1000*999*998*997*996*995*994*993*992*991) * 990!
990! = (990*989*988*987*986*985*984*983*982*981) * 980!

The other solution is just to swap over to string math.  It’s slow, but precise and reliable.
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: STxAxTIC on October 27, 2019, 11:38:26 am
I don't see offhand how doing the multiplication in groups would be better without resorting to logs or other manicuring tricks, I'll think about it.

String math does get right to the answer eventually, but in real life that's still too slow (for numbers that matter in useful science). In the biz we use Stirling's approximation for this stuff (see top post).

EDIT

Forgot to mention that doing it backwards is probably better than forwards... As in start with 1*2*3*4, rather than N*(N-1)*(N-2).
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: jack on February 04, 2020, 03:57:29 am
here's a solution without logarithms
Code: [Select]
DIM a AS _FLOAT
a = 1
n = 1000
IF n < 1755 THEN
    FOR k = 1 TO n
        a = a * k
    NEXT
    e = 0
    WHILE a > 1F308
        a = a / 10##
        e = e + 1
    WEND

    f$ = STR$(a)
    d = INSTR(f$, "D")
    d$ = MID$(f$, d + 1)
    e = e + VAL(MID$(f$, d + 1))
    f$ = LEFT$(f$, LEN(f$) - LEN(d$) + 1) + LTRIM$(STR$(e))

    PRINT "factorial of "; n; " = "; f$
ELSE
    PRINT n; " is too large"
END IF
it works up to n = 1754
I noticed that the 32-bit version of QB64 is more accurate,  the 64-bit version doesn't preserve the precision of a _Float but rather that of a double
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: Richard Frost on February 04, 2020, 07:01:40 am
Is there a prize for the weirdest effort?  I can only get to 170!, but I do use recursion, which is deadly short.  The
code for a single factorial could be 5 lines.  So I jazzed up my pathetic numerical skill with a font guaranteed  to make
people reach for their glasses.
Code: QB64: [Select]
  1. DEFINT A-Z
  2. COMMON SHARED n##, xx, yy, y
  3. DIM SHARED p(12, 4)
  4. DATA "0",7,5,5,5,7,"1",2,6,2,2,7,"2",7,1,7,4,7,"3",7,1,7,1,7,"4",5,5,7,1,1,"5",7,4,7,1,7
  5. DATA "6",7,4,7,5,7,"7",7,1,1,1,1,"8",7,5,7,5,7,"9",7,5,7,1,7,".",0,0,0,0,2,"+",0,2,7,2,0
  6. xx = 40: yy = 58: y = yy: SCREEN 12
  7. FOR i = 0 TO 11 '                                     load font
  8.     READ g$: FOR j = 0 TO 4: READ z: p(i, j) = z * 4096: NEXT
  10. FOR what = 1 TO 171
  11.     n## = 1: fact what: t$ = STR$(what) + STR$(n##)
  12.     IF what < 171 THEN tinyprint t$ ELSE LOCATE 25, 50: PRINT t$
  15. SUB fact (t)
  16.     IF t = 1 THEN EXIT SUB
  17.     n## = n## * t: fact t - 1 '                       look Mom, recursion!
  19. SUB tinyprint (t$) STATIC
  20.     x = xx
  21.     FOR i = 1 TO LEN(t$)
  22.         d$ = MID$(t$, i, 1): d = VAL(d$)
  23.         eh = INSTR(".+ ", d$): IF eh THEN d = eh + 9
  24.         FOR j = 0 TO 4
  25.             LINE (x, y + j)-(x + 4, y + j), 15, , p(d, j)
  26.         NEXT j
  27.         x = x + 4
  28.     NEXT i
  29.     y = y + 8: IF y > 400 THEN y = yy: xx = xx + 120
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: SMcNeill on February 04, 2020, 07:41:29 am
Quote from: jack on February 04, 2020, 03:57:29 am
I noticed that the 32-bit version of QB64 is more accurate,  the 64-bit version doesn't preserve the precision of a _Float but rather that of a double

As I explained elsewhere, this type of statement is very OS/compiler/flag specific.  _FLOATs translate to long doubles and, as such, play under their rules for how they’re implemented in C.

 
Quote from: SMcNeill on February 04, 2020, 07:34:40 am
_FLOAT isn’t “broke” on 64-bit versions.  You need to take a moment to read up on the long double variable type to try and understand the difference in behavior: https://en.wikipedia.org/wiki/Long_double

The precision results you see with _FLOAT depends heavily upon your system, the compiler you’re using, and sometimes even which option flags you have set with your compiler.

It’s not anything inside QB64 which is “broken”.  It’s simply the nature of a long double data type (_FLOAT).

_FLOATs are *required* to be as precise as doubles, but often offer more precision than that.  You can’t expect anything more from your programs than that, unless you want to write your own math routines or go to string math.
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: jack on February 04, 2020, 08:27:23 am
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
Quote
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.
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: david_uwi on February 04, 2020, 01:52:09 pm
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
  21. FOR j = m TO mx
  22.     PRINT USING "#"; f(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.

Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: bplus on February 04, 2020, 02:03:30 pm
Oh man that is nice, that first one stays accurate doesn't it?
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: david_uwi on February 04, 2020, 02:16:34 pm
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.
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: bplus on February 04, 2020, 03:14:31 pm
Well I may now have to take another look at string math, after seeing this and Petr's recent post.
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: SMcNeill on February 04, 2020, 04:28:42 pm
Quote from: jack on February 04, 2020, 08:27:23 am
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#

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. 

Quote
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.
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: STxAxTIC on February 04, 2020, 06:19:39 pm
Glad you found us david! I remember code similar to this from *.net in 2014, maybe sometime in June. Nice and fast!
Title: Re: Math challenge: Evaluate big factorials like 1000!
Post by: david_uwi on February 04, 2020, 06:34:45 pm
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.
Text Only | Text with Attachments