Math challenge: Evaluate big factorials like 1000!

[STxAxTIC]

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.

[Ashish]

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.#))
6. NEXT
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?

[STxAxTIC]

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

[jack]

in this particular case you could have dispensed with logarithms and used _FLOAT

[STxAxTIC]

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.

[jack]

sorry that I spoiled the fun
the _FLOAT range should be 3.4E-4932 to 1.1E+4932

[STxAxTIC]

Actually float starts failing somewhere before 175-factorial, it'll never make to 1000 unless this code is too tight:

Code: QB64: [Select]

1. DIM a AS _FLOAT
2. a = 1
3. FOR k = 1 TO 165
4.     a = a * k
5. NEXT
6. PRINT a

Hat remains tipped in Ashish's direction.

[jack]

@STxAxTIC
you are right, somehow _FLOAT is not an extended 80-bit float anymore or at least don't behave that way

[jack]

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

[STxAxTIC]

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.

[SMcNeill]

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.

[STxAxTIC]

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

[jack]

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

[Richard Frost]

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
9. 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$
13. NEXT
14. SLEEP
15. SUB fact (t)
16.     IF t = 1 THEN EXIT SUB
17.     n## = n## * t: fact t - 1 '                       look Mom, recursion!
18. END SUB
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
30. END SUB

[SMcNeill]

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.

 

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