Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #13362 > unrolled thread

Reverse factorial

Back to article view | Back to comp.lang.forth

Contents

  Reverse factorial "Ed" <invalid@nospam.com> - 2012-06-29 11:51 +1000
    Re: Reverse factorial mhx@iae.nl (Marcel Hendrix) - 2012-06-30 02:05 +0200
      Re: Reverse factorial "Ed" <invalid@nospam.com> - 2012-06-30 12:26 +1000
      Re: Reverse factorial "Ed" <invalid@nospam.com> - 2012-07-01 23:17 +1000
        Re: Reverse factorial mhx@iae.nl (Marcel Hendrix) - 2012-07-01 16:51 +0200
          Re: Reverse factorial "Ed" <invalid@nospam.com> - 2012-07-02 05:40 +1000
            Re: Reverse factorial Paul Rubin <no.email@nospam.invalid> - 2012-07-02 22:08 -0700
              Re: Reverse factorial "Ed" <invalid@nospam.com> - 2012-07-03 17:02 +1000
              Re: Reverse factorial mhx@iae.nl (Marcel Hendrix) - 2012-07-03 21:14 +0200
                Re: Reverse factorial Paul Rubin <no.email@nospam.invalid> - 2012-07-03 12:44 -0700
                  Re: Reverse factorial mhx@iae.nl (Marcel Hendrix) - 2012-07-03 22:37 +0200
                    Re: Reverse factorial "Ed" <invalid@nospam.com> - 2012-07-05 23:14 +1000

#13362 — Reverse factorial

Here's a reverse factorial routine to wile away your free time.

As the routine is highly f/p intensive it can double as a benchmark
for testing relative performance between systems and implementation
details such as separate/common/hardware f/p stacks.

--

\ N!.f
\
\ Forth version adapted from ...
\ http://www.xpl0.org/guess.html
\
\ N!.xpl         2-Mar-2010
\ Use 32-bit XPL0 (xpx.bat) for speed and accuracy
\ include c:\cxpl\codes;
\ real    N, A;
\ [N:= 0.0;  A:= 0.0;
\ repeat  N:= N + 1.0;
\         A:= A + Log(N);
\ until A >= 1E9-1.0;
\ RlOut(0, N);
\ ];
\

15 set-precision

1e9 fconstant #digits

: N! ( -- )
  0.0e ( N)
  0.0e ( A)
  begin
    fswap 1.0e f+ fswap
    fover flog f+
    fdup [ #digits 1.0e f- ] fliteral f< 0=
  until
  ( N A) fdrop f.
;

: run ( -- )
  cr cr ." Smallest number with a factorial of at least "
  #digits fs. ." digits is ... " cr n!
;

run

--

DX-Forth 3.98  2012-06-26

Forth-94
80387 15-digit floating point (common stack)

include n!

Smallest number with a factorial of at least 1.E9 digits is ...
130202809.  ok

[toc] | [next] | [standalone]

#13397

"Ed" <invalid@nospam.com> writes Re: Reverse factorial

> Here's a reverse factorial routine to wile away your free time.

> As the routine is highly f/p intensive it can double as a benchmark
> for testing relative performance between systems and implementation
> details such as separate/common/hardware f/p stacks.

Not this one.

[..]
> DX-Forth 3.98  2012-06-26

> Forth-94
> 80387 15-digit floating point (common stack)

> include n!

> Smallest number with a factorial of at least 1.E9 digits is ...
> 130202809.  ok

[..]

FORTH> runs
Smallest number with a factorial of at least 10^6 digits is 205023   \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^7 digits is 1723508   \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^8 digits is 14842907   \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^9 digits is 130202809   \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^10 digits is 1158787578   \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^11 digits is 10433891464   \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^12 digits is 94851898541   \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^13 digits is 869200494600   \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^14 digits is 8019346203786   \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^15 digits is 74419210652836   \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^16 digits is 694100859679692   \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^17 digits is 6502464891216880   \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^18 digits is 61154108320430276   \ 0.000 seconds elapsed. ok

-marcel

-- ----------------------------------------------------------
: LOG(N!) ( n -- ) ( F: -- log[n!] )
	S>F FLOCAL n
	PI*2 n F* 				FLOG F2/
	n 1e FEXP F/ 		 		FLOG n F* F+
	12e n F* 1/F F1+     
	288e n FSQR F* 1/F 		 	F+
	139e 512840e n F* n F* n F*  F/ 	F-
	571e 2488320e n FSQR FSQR F* F/  	F-  FLOG F+ ;

: FIND-n! ( F: #digits -- n )
	FLOCAL #digits
	#digits F>S 1 LOCALS| low high |
	BEGIN  
	  low high + 2/ LOG(N!) #digits 
	  F>= IF  low high + 2/  TO high  
	    ELSE  low high + 2/  TO low
	   ENDIF
	  low 1+ high >=
	UNTIL
	high ;

: run ( #digits -- )
	CR ." Smallest number with a factorial of at least 10^"
	DUP . ." digits is " timer-reset 10e S>F F** FIND-n! . ."   \ " .elapsed ;

: runs ( -- )  #19 6 DO  I run  LOOP ;

[toc] | [prev] | [next] | [standalone]

#13398

Marcel Hendrix wrote:
> "Ed" <invalid@nospam.com> writes Re: Reverse factorial
>
> > Here's a reverse factorial routine to wile away your free time.
>
> > As the routine is highly f/p intensive it can double as a benchmark
> > for testing relative performance between systems and implementation
> > details such as separate/common/hardware f/p stacks.
>
> Not this one.
> ...

Clearly you need a slower machine :)

The lack of variables and bare minimum of non-float instructions
made the routine suitable for testing my forth system with and
without a separate f/p stack.  The common stack version won
by a factor of two.  I've yet to figure out why :)

[toc] | [prev] | [next] | [standalone]

#13417

Marcel Hendrix wrote:
>
> [..]
>
> FORTH> runs
> Smallest number with a factorial of at least 10^6 digits is 205023   \ 0.000 seconds elapsed.
> Smallest number with a factorial of at least 10^7 digits is 1723508   \ 0.000 seconds
> elapsed.
> Smallest number with a factorial of at least 10^8 digits is 14842907   \ 0.000 seconds
> elapsed.
> Smallest number with a factorial of at least 10^9 digits is 130202809   \ 0.000 seconds
> elapsed.
> Smallest number with a factorial of at least 10^10 digits is 1158787578   \ 0.000 seconds
> elapsed.
 > ...

I get these results using your algorithm

Smallest number with a factorial of at least 10^0 digits: 1.
Smallest number with a factorial of at least 10^1 digits: 10.
Smallest number with a factorial of at least 10^2 digits: 70.
Smallest number with a factorial of at least 10^3 digits: 450.
Smallest number with a factorial of at least 10^4 digits: 3249.
Smallest number with a factorial of at least 10^5 digits: 25206.
Smallest number with a factorial of at least 10^6 digits: 205023.

While the original (slow) routine gives

Smallest number with a factorial of at least 1.E0 digits: 1.
Smallest number with a factorial of at least 1.E1 digits: 13.
Smallest number with a factorial of at least 1.E2 digits: 70.
Smallest number with a factorial of at least 1.E3 digits: 450.
Smallest number with a factorial of at least 1.E4 digits: 3249.
Smallest number with a factorial of at least 1.E5 digits: 25206.
Smallest number with a factorial of at least 1.E6 digits: 205022.

[toc] | [prev] | [next] | [standalone]

#13418

"Ed" <invalid@nospam.com> writes Re: Reverse factorial

> Marcel Hendrix wrote:
[..]
> I get these results using your algorithm

> Smallest number with a factorial of at least 10^0 digits: 1.
> Smallest number with a factorial of at least 10^1 digits: 10.
[..]
> Smallest number with a factorial of at least 10^6 digits: 205023.

> While the original (slow) routine gives

> Smallest number with a factorial of at least 1.E0 digits: 1.
> Smallest number with a factorial of at least 1.E1 digits: 13.
[..]
> Smallest number with a factorial of at least 1.E6 digits: 205022.
[..]

According to Wolfram Alpha, 

log10(13!) = 9.79428031638948
log10(14!) = 10.940408352067719

Therefore the result for 10^1 digits should be 14, not 13. My program
is wrong, but not because log(n!) is faulty (see below).

According to Wolfram Alpha, 

log10(205022!) = 999999.08681421082059069559626773741428638233100481270
log10(205023!) = 1.00000439861679486551512255600216645545545382220e6

Therefore 205023 as the "Smallest number with a factorial of at least 10^6 
digits" is correct.

The fix to my program is:

: FIND-n! ( F: #digits -- n )
...
	#digits F>S  3 2 */  1 LOCALS| low high |
...
	high ;

The result for 10^0 is wrong for both programs (0! = 1).
As should be obvious, I didn't test my code for 10^0 and 10^1 digits. 

-marcel

[toc] | [prev] | [next] | [standalone]

#13426

Marcel Hendrix wrote:
> ...
> According to Wolfram Alpha,
>
> log10(13!) = 9.79428031638948
> log10(14!) = 10.940408352067719
>
> Therefore the result for 10^1 digits should be 14, not 13. My program
> is wrong, but not because log(n!) is faulty (see below).

14!  produces an integer with 11 digits - more than we need
13!  produces an integer with 10 digits - better
12!  produces an integer with  9 digits - too low

> According to Wolfram Alpha,
>
> log10(205022!) = 999999.08681421082059069559626773741428638233100481270
> log10(205023!) = 1.00000439861679486551512255600216645545545382220e6
>
> Therefore 205023 as the "Smallest number with a factorial of at least 10^6
> digits" is correct.

Given the previous paragraph, I wouldn't be too sure.

> The fix to my program is:
>
> : FIND-n! ( F: #digits -- n )
> ...
> #digits F>S  3 2 */  1 LOCALS| low high |
> ...
> high ;

Introducing a bias to overcome the error in the approximation
algorithm at small values.  Commonly known as fudging the
result.

> The result for 10^0 is wrong for both programs (0! = 1).

10^0 digits = 1 digit which satisfies 1!

0! = 1 is a convention supported by mathematical logic
(the same logic that gave you infinity along with some
infinities being greater than others :)

0! = 1 would not be directly computable and a program
would need to support it as a special case.

[toc] | [prev] | [next] | [standalone]

#13462

"Ed" <invalid@nospam.com> writes:
>> log10(205022!) = 999999.08681421082059069559626773741428638233100481270
>> log10(205023!) = 1.00000439861679486551512255600216645545545382220e6
>>
>> Therefore 205023 as the "Smallest number with a factorial of at least 10^6
>> digits" is correct.

Of course 10**999999.08 has 1e6 digits just like 10^9 has 10 digits.

I calculated 205022! exactly with GHC (uses GMP underneath).  This took
about 8 seconds including the output conversion (which used around 300MB
of memory, doesn't seem good).  The first digits are:

 1221277091680869925277501633213930559467939034445151838141073953...

I think for inverse factorial of large numbers, it probably is pretty
accurate to use a numerical rootfinder on the first few terms of the
Stirling series:

http://en.wikipedia.org/wiki/Stirling%27s_approximation#Speed_of_convergence_and_error_estimates 

Even just the first (1/12n) correction term makes Stirling's
approximation pretty accurate.

[toc] | [prev] | [next] | [standalone]

#13463

Paul Rubin wrote:
> "Ed" <invalid@nospam.com> writes:

That was Marcel.  Maths is beyond me.

> >> log10(205022!) = 999999.08681421082059069559626773741428638233100481270
> >> log10(205023!) = 1.00000439861679486551512255600216645545545382220e6
> >>
> >> Therefore 205023 as the "Smallest number with a factorial of at least 10^6
> >> digits" is correct.
>
> Of course 10**999999.08 has 1e6 digits just like 10^9 has 10 digits.
>
> I calculated 205022! exactly with GHC (uses GMP underneath).  This took
> about 8 seconds including the output conversion (which used around 300MB
> of memory, doesn't seem good).  The first digits are:
>
>  1221277091680869925277501633213930559467939034445151838141073953...
>
> I think for inverse factorial of large numbers, it probably is pretty
> accurate to use a numerical rootfinder on the first few terms of the
> Stirling series:
>
> http://en.wikipedia.org/wiki/Stirling%27s_approximation#Speed_of_convergence_and_error_estimates
>
> Even just the first (1/12n) correction term makes Stirling's
> approximation pretty accurate.

[toc] | [prev] | [next] | [standalone]

#13482

Paul Rubin <no.email@nospam.invalid> writes Re: Reverse factorial
[..]
>>> log10(205022!) = 999999.08681421082059069559626773741428638233100481270
>>> log10(205023!) = 1.00000439861679486551512255600216645545545382220e6

>>> Therefore 205023 as the "Smallest number with a factorial of at least 10^6
>>> digits" is correct.

> Of course 10**999999.08 has 1e6 digits just like 10^9 has 10 digits.

It is beyond my semantic capabability to understand what "Smallest 
number with a factorial of at least 10^x digits" means exactly. Nail the
constraint in terms of log(n!) and the program will probably be correct.

[..]
> I think for inverse factorial of large numbers, it probably is pretty
> accurate to use a numerical rootfinder on the first few terms of the
> Stirling series:

> http://en.wikipedia.org/wiki/Stirling%27s_approximation#Speed_of_convergence_and_error_estimates 

> Even just the first (1/12n) correction term makes Stirling's
> approximation pretty accurate.

Hmm, apparently that algorithm doesn't jump out of the Forth for you :-)

The upperbound of log(n!) is not strictly n (for low n where I was not 
intending to use the words), and so the binary search walked against the
bound before finding the root.

-marcel

[toc] | [prev] | [next] | [standalone]

#13486

mhx@iae.nl (Marcel Hendrix) writes:
> It is beyond my semantic capabability to understand what "Smallest 
> number with a factorial of at least 10^x digits" means exactly. 

n = minimum { k | length(decimal representation of k!) >= 10^x }

For d>0, the length of the decimal representation of d is  
1 + floor(log10 d).

[toc] | [prev] | [next] | [standalone]

#13488

Paul Rubin <no.email@nospam.invalid> writes Re: Reverse factorial

> mhx@iae.nl (Marcel Hendrix) writes:
>> It is beyond my semantic capabability to understand what "Smallest 
>> number with a factorial of at least 10^x digits" means exactly. 

> n = minimum { k | length(decimal representation of k!) >= 10^x }

> For d>0, the length of the decimal representation of d is  
> 1 + floor(log10 d).

Thanks!

FORTH> runs
Smallest number with a factorial of at least 10^00 digits is  0.0000000000000000000e+0000  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^01 digits is  1.3000000000000000000e+0001  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^02 digits is  7.0000000000000000000e+0001  \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^03 digits is  4.5000000000000000000e+0002  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^04 digits is  3.2490000000000000000e+0003  \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^05 digits is  2.5206000000000000000e+0004  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^06 digits is  2.0502200000000000000e+0005  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^07 digits is  1.7235080000000000000e+0006  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^08 digits is  1.4842907000000000000e+0007  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^09 digits is  1.3020280900000000000e+0008  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^10 digits is  1.1587875780000000000e+0009  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^11 digits is  1.0433891464000000000e+0010  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^12 digits is  9.4851898541000000000e+0010  \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^13 digits is  8.6920049460000000000e+0011  \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^14 digits is  8.0193462037860000000e+0012  \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^15 digits is  7.4419210652836000000e+0013  \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^16 digits is  6.9410085967969200000e+0014  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^17 digits is  6.5024648912168800000e+0015  \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^18 digits is  6.1154108320430276000e+0016  \ 0.001 seconds elapsed.
Smallest number with a factorial of at least 10^19 digits is  5.7713453304452275000e+0017  \ 0.000 seconds elapsed.
Smallest number with a factorial of at least 10^20 digits is  5.4635317748670943960e+0018  \ 0.000 seconds elapsed. ok

-marcel

[toc] | [prev] | [next] | [standalone]

#13604

Marcel Hendrix wrote:
> Paul Rubin <no.email@nospam.invalid> writes Re: Reverse factorial
>
> > mhx@iae.nl (Marcel Hendrix) writes:
> >> It is beyond my semantic capabability to understand what "Smallest
> >> number with a factorial of at least 10^x digits" means exactly.
>
> > n = minimum { k | length(decimal representation of k!) >= 10^x }
>
> > For d>0, the length of the decimal representation of d is
> > 1 + floor(log10 d).
>
> Thanks!
>
> FORTH> runs
> Smallest number with a factorial of at least 10^00 digits is  0.0000000000000000000e+0000  \
> 0.000 seconds elapsed.
> Smallest number with a factorial of at least 10^01 digits is 1.3000000000000000000e+0001  \
> 0.000 seconds elapsed.
> Smallest number with a factorial of at least 10^02 digits is  7.0000000000000000000e+0001  \
> 0.001 seconds elapsed.
> Smallest number with a factorial of at least 10^03 digits is  4.5000000000000000000e+0002  \
> 0.000 seconds elapsed.
> Smallest number with a factorial of at least 10^04 digits is 3.2490000000000000000e+0003  \
> 0.001 seconds elapsed.
> Smallest number with a factorial of at least 10^05 digits is  2.5206000000000000000e+0004  \
> 0.000 seconds elapsed.
> Smallest number with a factorial of at least 10^06 digits is  2.0502200000000000000e+0005  \
> 0.000 seconds elapsed.
> ...

What program generated these?  Was the result for 10^00 calculated
or plugged-in?

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.forth

csiph-web