Marcel Hendrix wrote: > INCLUDE ../bignum/factor.frt
If you have a 64 bit Forth (and I know you have ;-), you don't need bignums for that. All intermediate results have at most 20 digits, and therefore easily fit in a 128 bit number; the sum with at most 13 digits fits in a 64 bit number. This one should be faster (how much depends on how efficient your gs^mod actually is). In Gforth with a 2GHz Athlon64, it takes about 1.5ms to run 1000t. I guess iForth64 should do it in less than 1ms, even when the division for the modulus might be the limiting factor.
\ sum of the last 10 digits of n^n for n=1..1000
Cell 8 < [IF] .( Needs a 64 bit Forth!) cr true abort [THEN]
&10000000000 Constant mod#
: *mod ( a b -- c ) um* mod# um/mod drop ;
: **mod ( x n -- n ) >r 1 swap BEGIN r@ 1 and IF tuck *mod swap THEN r> 2/ dup WHILE >r dup *mod REPEAT 2drop ;
: e48 ( n -- g ) 0 swap 1+ 1 DO I dup **mod + LOOP mod# mod ;
: Euler48 ( -- ) CR ." The last ten digits of the series " ." 1^1 + 2^2 + 3^3 + ... + 1000^1000 are " 1000 e48 . cr ;
Bernd Paysan <bernd.pay...@gmx.de> writes: >In Gforth with a 2GHz Athlon64, it takes about >1.5ms to run 1000t. I guess iForth64 should do it in less than 1ms, even >when the division for the modulus might be the limiting factor.
If you timed it with the development version of gforth-fast, I doubt it. I guesstimate that the program spends more then 1ms on UM/MOD on your machine.
>&10000000000 Constant mod#
>: *mod ( a b -- c ) um* mod# um/mod drop ;
Wasn't there a way to do integer division (and consequently compute the modulus) by multiplication with the inverse?
Anton Ertl wrote: > Bernd Paysan <bernd.pay...@gmx.de> writes: >>In Gforth with a 2GHz Athlon64, it takes about >>1.5ms to run 1000t. I guess iForth64 should do it in less than 1ms, even >>when the division for the modulus might be the limiting factor.
> If you timed it with the development version of gforth-fast, I doubt > it. I guesstimate that the program spends more then 1ms on UM/MOD on > your machine.
I had it run with the normal gforth. Gforth-fast on the 2GHz Athlon64: 955µs, on a 2.6GHz Core 2 Quad 650µs. So actually iForth64 must be under one 1ms, too ;-).
>>&10000000000 Constant mod#
>>: *mod ( a b -- c ) um* mod# um/mod drop ;
> Wasn't there a way to do integer division (and consequently compute > the modulus) by multiplication with the inverse?
Won't help much, since we need to divide a 128 bit number, and that requires too many multiplications to get the result (the inverse also must be a 128 bit number to be accurate enough, so we have four multiplications to divide down, and then another multiplication to restore the modulus from the fractional part). It's probably possible to use a 64 bit number and the right amount of bitshift, but for a Euler problem solution, this is too much work ;-).
I've tried something else: Using fractional numbers, i.e. scaling by 1e-10 and having the decimal point on bit 64. The problem is the same: If you use just a 64 bit inverse of 1e10, you don't have much accuracy left. Using 128 bits requires too many multiplications, and I still didn't get it accurate enough.
> Find the last ten digits of the series, 1^1 + 2^2 + 3^3 + ... + 1000^1000. >*)
>: INIT ( -- ) S" 10000000000" s2 >GIANT S" 0" s3 >GIANT ; >: TERM ( n -- ) DUP 0 <# #S #> s1 >GIANT s1 SWAP s2 GS^MOD s3 s1 GG+ ; >: 1000T ( -- g ) INIT #1001 1 DO I TERM LOOP s3 s2 GGMOD s3 .GIANT ;
>: Euler48 ( -- ) > CR ." The last ten digits of the series 1^1 + 2^2 + 3^3 + ... + 1000^1000 are " 1000T ;
>: .ABOUT CR ." Euler48 -- Find the last ten digits of the series, " > CR ." 1^1 + 2^2 + 3^3 + ... + 1000^1000." ;
> .ABOUT
>\ FORTH> euler48 >\ The last ten digits of the series 1^1 + 2^2 + 3^3 + ... + 1000^1000 are xxxxxxxxxx >\ 0.011 seconds elapsed. ok
Please don't publish spoilers without spoiler in the subject line.
Groetjes Albert
-- -- Albert van der Horst, UTRECHT,THE NETHERLANDS Economic growth -- like all pyramid schemes -- ultimately falters. albert@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst
In article <5jfbf5-14d....@annette.mikron.de>, Bernd Paysan <bernd.pay...@gmx.de> wrote: <SNIP>
>I've tried something else: Using fractional numbers, i.e. scaling by 1e-10 >and having the decimal point on bit 64. The problem is the same: If you use >just a 64 bit inverse of 1e10, you don't have much accuracy left. Using 128 >bits requires too many multiplications, and I still didn't get it accurate >enough.
This too must work: If we are prepared to do N^2 operations instead of N log N, for an exponentiation, we can calculate 877^877 without needing more precision than double on a 32 bit Forth. Repeatedly multiply with 877. This remains under 1000*10^10. Then calculate the last 10 digits by first calculating the last 5, then the last 5 of the remainder. (Using FM/MOD) Combine, rinse, repeat. (OK this is a very limited form of multiple precision.)
(Oeps, we are now into 1000^3 operations. This may become an overnighty.)
Some of the problems give me the idea I'm cheating when using a 64 bit Forth, in the sense that you don't need to address the real challenge.
-- -- Albert van der Horst, UTRECHT,THE NETHERLANDS Economic growth -- like all pyramid schemes -- ultimately falters. albert@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst
William James <w_a_x_...@yahoo.com> wrote Re: Euler problem #48
>> \ FORTH> euler48 >> \ The last ten digits of the series 1^1 + 2^2 + 3^3 + ... + 1000^1000 are xxxxxxxxxx >> \ 0.011 seconds elapsed. ok > In PARI/gp. 0 ms. It yields more than 10 digits. > sum(i=1,1000,lift( Mod(i,10^10)^i ) )
You mean, it calculates the wrong answer, faster?
Sum = Mod( sum(i=1,1000,lift( Mod(i,10^10)^i ) ), 10^10 )
Methink that for i<1001 Mod(i,10^10) always equals i. What do I miss here?
>Fixes that, surely.
>-marcel
-- -- Albert van der Horst, UTRECHT,THE NETHERLANDS Economic growth -- like all pyramid schemes -- ultimately falters. albert@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst
