[SPOILER] Re: Euler problem #187
Anton Ertl
and it seems to be one of those that had not been solved by any of the
other registered Forth users (there is an advantage to us being so
few).
On a 3GHz Xeon 5160 takes 3.65s on gforth-fast and 3.32s on iForth. I
guess that the Sieve is running into the memory wall, that's why the
advantage of iForth is so small here.
As usual, you can download the program from
<http://www.complang.tuwien.ac.at/forth/programs/euler/187.fs>.
\ Problem 187:
\ A composite is a number containing at least two prime factors. For
\ example, 15 = 3 × 5; 9 = 3 × 3; 12 = 2 × 2 × 3.
\ There are ten composites below thirty containing precisely two, not
\ necessarily distinct, prime factors: 4, 6, 9, 10, 14, 15, 21, 22,
\ 25, 26.
\ How many composite integers, n < 10^8, have precisely two, not
\ necessarily distinct, prime factors?
\ Solution
\ I guess something smart can be done with prime-counting-only code,
\ but here I am going for a simple solution: Essentially, the numbers
\ we seek are of the form n=p1*p2, where p1<=p2, p1<sqrt(limit),
\ n<limit, p1 and p2 prime. I count them by generating the primes up
\ to sqrt(limit), then for each such prime p1 count the primes p2 with
\ p1<=p2<limit/p1. Since I'm lazy, I generate the primes up to
\ limit/2 and use them for counting.
\ As I said, I'm lazy, so I'll just adapt the Byte Sieve, as bad as it
\ may be.
\ If you use an ANS Forth system other than Gforth, get
\ <http://www.complang.tuwien.ac.at/forth/compat.zip> and uncomment this:
\ s" compat/loops.fs" included
\ s" compat/anslocal.fs" included
100000000 constant limit
limit 2/ constant plimit \ a little more than needed
plimit 2/ allocate throw constant flags
FLAGS plimit 2/ + CONSTANT EFLAG
: PRIMES ( -- n )
FLAGS plimit 2/ 1 FILL
3 EFLAG FLAGS DO
I C@ IF
DUP I + DUP EFLAG < IF
EFLAG SWAP DO
0 I C!
DUP +LOOP
ELSE
DROP THEN
THEN
2 + LOOP
DROP ;
: sqrt ( n1 -- n2 )
s>d d>f fsqrt f>d drop ;
: count-p2 { n1 -- n2 }
0 limit 1- n1 / 3 - 2/ 1+ n1 3 - 2/ +do
flags i + c@ + loop ;
: solution ( -- u )
1 \ counter, initialize with 1 to count 4=2*2
-1 limit sqrt 2/ -do
flags i + c@ if
i 2* 3 + count-p2 + then
1 -loop
2 count-p2 + ;
primes solution . cr
\ You can speed this up by a little by remembering what COUNT-P2
\ produced last time, and only actually counting the primes that were
\ not counted last time, but I am too lazy for that. Besides, the
\ speedup is only small, because PRIMES takes most of the execution
\ time already (3.14s out of 3.65s on a 3GHz Xeon 5160 with
\ gforth-fast); strange, I would have expected SOLUTION to take a
\ similar amount of time.
--
M. Anton Ertl http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
New standard: http://www.forth200x.org/forth200x.html
EuroForth 2008: http://www.complang.tuwien.ac.at/anton/euroforth/ef08.html
Marcel Hendrix
> That was another relatively easy one (or I would not have done it),
Your code, pasted into Gforth-fast, generates a number that is rejected by
the Euler solution checker?
-marcel
cac
> Your code, pasted into Gforth-fast, generates a number that is rejected by
> the Euler solution checker?
>
Picky, picky. I suspect that:
1 \ counter, initialize with 1 to count 4=2*2
sould be:
0 \ counter, initialize with 1 to count 4=2*2
-- Charles
Anton Ertl
If I run it on a 64-bit system, the result is correct. If I run it on
a 32-bit system, or if I run it on iForth, it is too high by 1. It's
not obvious where the bug is.
In addition, depending on which version you are using, you may see the
gcc-2.95/gforth-fast transcendentals bug; work around that by defining
: fsqrt 0 fsqrt drop ;
- anton
Marcel Hendrix
> m...@iae.nl (Marcel Hendrix) writes:
> >an...@mips.complang.tuwien.ac.at (Anton Ertl) writes Re: [SPOILER] Re: Euler problem #187
[..]
>>Your code, pasted into Gforth-fast, generates a number that is rejected by
>>the Euler solution checker?
> If I run it on a 64-bit system, the result is correct. If I run it on
> a 32-bit system, or if I run it on iForth, it is too high by 1. It's
> not obvious where the bug is.
You have strong evidence that the 64-bit system (and the algorithm)
are perfect? I can not read the code with these funny loops.
As there are only about 18 million numbers, one could write them to disk
and find the odd one out by comparing 64-bit and 32-bit output files.
Maybe you encountered a CPU/FPU bug in 64-bit mode :-)
-marcel
Anton Ertl
>an...@mips.complang.tuwien.ac.at (Anton Ertl) writes Re: [SPOILER] Re: Euler problem #187
>
>> m...@iae.nl (Marcel Hendrix) writes:
>> >an...@mips.complang.tuwien.ac.at (Anton Ertl) writes Re: [SPOILER] Re: Euler problem #187
>[..]
>>>Your code, pasted into Gforth-fast, generates a number that is rejected by
>>>the Euler solution checker?
>
>> If I run it on a 64-bit system, the result is correct. If I run it on
>> a 32-bit system, or if I run it on iForth, it is too high by 1. It's
>> not obvious where the bug is.
I found and fixed the bug:
: count-p2 { n1 -- n2 }
0 limit 1- n1 / 3 - 2/ 1+ n1 3 - 2/ 0 max +do
flags i + c@ + loop ;
Without the "0 MAX", for n1=2, this accessed element -1 of FLAGS; this
element corresponds to the number 1, but is obviously not initialized.
Apparently it happened that 32-bit Gforth and iForth had a 1 there and
64-bit Gforth had a 0 there. Strange, there could have been 254 other
values there.
>You have strong evidence that the 64-bit system (and the algorithm)
>are perfect?
Well, having it do the right thing for limit=30 and then getting the
result accepted by the Euler website was made me pretty confident, but
obviously that was overconfidence.
> I can not read the code with these funny loops.
Read up on them in
http://www.complang.tuwien.ac.at/forth/gforth/Docs-html/Counted-Loops.html
They are just improved versions of the regular counted loops. I guess
I would have made another error or two with the regular counted loops.
However, the reverse loop in SOLUTION can be replaced with a forward
loop here (which could be implemented with DO...LOOP), but if we want
to do the optimization I mentioned at the end of the program, we need
the reverse loop.
BTW, I now think that we can get a significant speedup by making
PRIMES stop when it reaches sqrt(plimit).
Marcel Hendrix
> I found and fixed the bug:
> : count-p2 { n1 -- n2 }
> 0 limit 1- n1 / 3 - 2/ 1+ n1 3 - 2/ 0 max +do
> flags i + c@ + loop ;
> Without the "0 MAX", for n1=2, this accessed element -1 of FLAGS; this
> element corresponds to the number 1, but is obviously not initialized.
> Apparently it happened that 32-bit Gforth and iForth had a 1 there and
> 64-bit Gforth had a 0 there. Strange, there could have been 254 other
> values there.
You are accessing the header of allocated memory here. With iForth for
Windows, there's debug info from C's heap manager there. Unloading and
reloading the program a few times will produce a different address,
therefore a different header, and therefore a different result
(because memory is not deallocated upon unloading in this program).
[..]
>> I can not read the code with these funny loops.
> Read up on them in
> http://www.complang.tuwien.ac.at/forth/gforth/Docs-html/Counted-Loops.html
> They are just improved versions of the regular counted loops. I guess
> I would have made another error or two with the regular counted loops.
They may be technically superior, but having to look up the 8 (?)
variants every time undoes any possible benefit for me. I really do
not understand what 'problem' the Gforth manual sees in
: test 0 3 do i . -1 +loop ; test 3 2 1 0 ok
-marcel
-- ----
FORTH> in
Euler187 -- How many composite integers, n < 10^8, have precisely two, not
-- necessarily distinct, prime factors?
ok
FORTH> flags 8 - 16 dump
$0E090028 00 00 00 00 20 00 09 0E--01 01 01 00 01 01 00 01
ok
FORTH> euler187
There are 17427272 composite integers, n < 10^8, having precisely two, not
necessarily distinct, prime factors.
0.250 seconds elapsed. ok
Anton Ertl
>an...@mips.complang.tuwien.ac.at (Anton Ertl) writes Re: [SPOILER] Re: Euler problem #187
>
>> Read up on them in
>
>> http://www.complang.tuwien.ac.at/forth/gforth/Docs-html/Counted-Loops.html
>
>> They are just improved versions of the regular counted loops. I guess
>> I would have made another error or two with the regular counted loops.
>
>They may be technically superior, but having to look up the 8 (?)
>variants every time undoes any possible benefit for me.
I think it's pretty mnemonic (see below), and it saves you remembering
and working around the traps that the Forth-94 counted loops lay out
for you.
Let's first look at the DOs: You can have DOs for up-counting loops
(+DO, U+DO), and DOs for down-counting loops (-DO, U-DO). The other
dimension is whether the limits are signed (+DO, -DO) or unsigned
(U+DO, U-DO). So this gives a table
signed unsigned
count up +DO U+DO
count down -DO U-DO
Next, let's look at the LOOPs: there is the usual LOOP and +LOOP for
counting up, and the new -LOOP for counting down.
I don't use counted loops often, but I don't need to look these words
up when I need them.
>I really do
>not understand what 'problem' the Gforth manual sees in
>
>: test 0 3 do i . -1 +loop ; test 3 2 1 0 ok
: test 3 0 do i . 1 +loop ; test 0 1 2 ok
If you don't see a problem, you are blind.
Here's a better solution for counting down:
: test 0 3 do i . 1 -loop ; test 3 2 1 ok
(Note that the compat/loops.fs version of -LOOP only works with -DO and
U-DO, not with DO or ?DO).
The other problems with your test loop construction is that it uses
DO, and the problem can be seen with:
: test1 do i . -1 +loop ;
0 0 test1
0 -3 test1
These can be avoided by using -DO.
There is only one excuse for using DO, and knowing that the parameters
are distinct is not one of them.
>FORTH> euler187
>There are 17427272 composite integers, n < 10^8, having precisely two, not
>necessarily distinct, prime factors.
>0.250 seconds elapsed. ok
It seems you have worked at getting the wrong answer way faster?
Marcel Hendrix
> an...@mips.complang.tuwien.ac.at (Anton Ertl) writes Re: [SPOILER] Re: Euler problem #187
>> I found and fixed the bug:
Never forget the old but nonetheless excellent work of others (in this
case Albert van der Horst).
-marcel
-- --------------------------------------------------------------------------------------
INCLUDE ../benchmar/benchpin.frt \ PI(N2) ( n1 -- n2 ) counts all primes below n1
(*
A composite is a number containing at least two prime factors. For
example, 15 = 3 * 5; 9 = 3 * 3; 12 = 2 * 2 * 3.
There are ten composites below thirty containing precisely two, not
necessarily distinct, prime factors: 4, 6, 9, 10, 14, 15, 21, 22,
25, 26.
How many composite integers, n < 10^8, have precisely two, not
necessarily distinct, prime factors?
Solution:
---------
This is PI(n/i) - pi(i) + 1 for all primes i not greater than sqrt(n).
*)
: solution ( -- n )
0 #10000 2 DO I Prime? IF #100000000 I / PI(N2) I PI(N2) - 1+ + ENDIF LOOP ;
: Euler187 ( -- )
CR ." There are " solution . ." composite integers, n < 10^8, having precisely two, not "
CR ." necessarily distinct, prime factors." ;
\ FORTH> euler187
\ There are 17xxxxxx composite integers, n < 10^8, having precisely two, not
\ necessarily distinct, prime factors.
\ 0.636 seconds elapsed. ok
Luca Masini
> Never forget the old but nonetheless excellent work of others (in this
> case Albert van der Horst).
[SNIPPED]
> -- --------------------------------------------------------------------------------------
> INCLUDE ../benchmar/benchpin.frt \ PI(N2) ( n1 -- n2 ) counts all primes below n1
[SNIPPED]
Is available the definition of PI(N2) ?
> Solution:
> ---------
> This is PI(n/i) - pi(i) + 1 for all primes i not greater than sqrt(n).
In J is a one-liner
+/ (_1&p:@:>.@(n%p:)-]) i._1&p:%:n [ n=:1e8
1742yyyy
where _1&p: is the number of prime less than the argument ( PI(N2) )
I'm trying to solve some problem using Forth now.
Luca.
Marcel Hendrix
> Marcel Hendrix wrote:
>> Never forget the old but nonetheless excellent work of others (in this
>> case Albert van der Horst).
> [SNIPPED]
>> -- --------------------------------------------------------------------------------------
>> INCLUDE ../benchmar/benchpin.frt \ PI(N2) ( n1 -- n2 ) counts all primes below n1
> [SNIPPED]
> Is available the definition of PI(N2) ?
I found something faintly resembling it:
ftp://ccreweb.org/software/kforth/examples/benchmarks/benchpin.4th
> > Solution:
> > ---------
> > This is PI(n/i) - pi(i) + 1 for all primes i not greater than sqrt(n).
> In J is a one-liner
> +/ (_1&p:@:>.@(n%p:)-]) i._1&p:%:n [ n=:1e81742yyyy
> where _1&p: is the number of prime less than the argument ( PI(N2) )
It takes some getting used too, I guess :-)
> I'm trying to solve some problem using Forth now.
-marcel
Albert van der Horst
>m...@iae.nl (Marcel Hendrix) writes Re: [SPOILER] Re: Euler problem #187
>
>> an...@mips.complang.tuwien.ac.at (Anton Ertl) writes Re: [SPOILER] Re: Euler problem #187
>
>>> I found and fixed the bug:
>
>Never forget the old but nonetheless excellent work of others (in this
>case Albert van der Horst).
>
>-marcel
>\ FORTH> euler187
>\ There are 17xxxxxx composite integers, n < 10^8, having precisely two, not
>\ necessarily distinct, prime factors.
>\ 0.636 seconds elapsed. ok
It is a benchmark intended to waste cycles in a recursive
program that actually calculates something.
Nice to see that it is up to production work.
I used it myself to find the 10,001 th prime, doing trial and error
by hand.
I solved a few of the harder problems.
The least significant non-zero digits of 10^12 ! (euler188) took me
halve a day.
The problem is that numbers ending in 5 together with an even number
give a zero. The bottom line is that you must account for all factors
5 and 2. On the other hand, you need not bother with the odd
numbers, they all cancel ... (which leaves 5*10^11 numbers to
multiply ).
The other one was the hyper exponentation:
1777 ^ 1777 ^ 1777 .. (1855 times). .. ^ 1777
The task is not to calculate all digits ;-) merely the last 5.
It turns out that some of the solutions were based on wrong
reasoning, but still give the right answer.
Both of these harder problems run in submillisecond time once the are
solved. The advantage of the harder problems is that the discussion is
not closed, such that I can show people some Forth code.
I have other things to do, but cannot help to tackle the colored
triangle, the capacitor and the triomino problem. The triomino problem
can be solved with techniques of the pentominos. A bitmap solution can
be spectacularly fast, and could show off Forth in the discussion.
The people at euler are responsive. On my suggestion they added the
following to the intended audience "professionals who want to keep
their programming and mathematics skills on the edge". Now Marcel,
Anton and myself need no longer be ashamed of working on those
problems.
I'm still at page 3 (of 4 dutch pages) with 14 solved.
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
Anton Ertl
>I have other things to do, but cannot help to tackle the colored
>triangle
I guess you mean Problem 189. I looked at this for half an hour or
so, and I think I have an approach to solve it, but am too lazy to
program it. Here's my approach:
Divide the triange into 4 equally sized subtriangles; each of the
subtriangles has 4 element triangles on each edge. There's a central
subtriangle and 3 corner subtriangles; they interface with the center
triangle along one edge, and not at all with each other.
Now, enumerate the possible colourings of one subtriangle, and for
each of the 81 colour combinations on one edge, remember how many
colourings have that edge.
Next, look at the interface of two subtriangles: for all compatible
combinations of edges, multiply the number of colourings represented
by the combinations and add up the products. That gives you the
number of colourings for the combination of two subtriangles.
Ideally, the three edges of the central subtriangle would be
independent enough that we can just apply the factor we determined
from adding the first corner triangle two times more to get the number
for the combination of all subtriangles.
But I suspect it won't be so easy. Then we have to deal with all the
egdes of the central subtriangle at once, i.e., 19683 potential
combinations. Should not be too bad, either.
cac
Hmmm. My approach was to divide the triangle into an upper half and a
lower half, and enumerate the boundary as you did. Where I get wedged
is the rotation and reflection issue. I need to figure out how many
non-distinct colourings there are, and subtract that from my answer;
or think about you solution for a while; I suspect it might be a much
better approach.
-- Charles
Anton Ertl
>Hmmm. My approach was to divide the triangle into an upper half and a
>lower half, and enumerate the boundary as you did. Where I get wedged
>is the rotation and reflection issue. I need to figure out how many
>non-distinct colourings there are, and subtract that from my answer;
As far as I understand the task, rotation and reflection is a
non-issue. If you enumerate all colourings for one orientation,
that's the answer.
If you had a way to create only one of, say, the two mirror images
(along some axis), then you would have to count it twice, unless the
colouring is symmetric. But with the usual enumerative approach we
get both mirror images anyway.
cac
> cac <c...@inbox.com> writes:
>> Hmmm. My approach was to divide the triangle into an upper half and a
>> lower half, and enumerate the boundary as you did. Where I get wedged
>> is the rotation and reflection issue. I need to figure out how many
>> non-distinct colourings there are, and subtract that from my answer;
>
> As far as I understand the task, rotation and reflection is a
> non-issue. If you enumerate all colourings for one orientation,
> that's the answer.
>
> If you had a way to create only one of, say, the two mirror images
> (along some axis), then you would have to count it twice, unless the
> colouring is symmetric. But with the usual enumerative approach we
> get both mirror images anyway.
From the problem statement:
A colouring C' which is obtained from a colouring C by rotation or
reflection is considered distinct from C unless the two are
identical.
How many distinct valid colourings are there for the above
configuration?
I read that as saying if a rotation or reflection is identical, it
is not distinct, and should not be counted. Am I confused as to
the interpretation?
For example, with a much smaller triangle, I generate the following
arrangements:
R R
R B G G B R
They are both generated by my algorithm, and add 2 to the enumeration.
Are they distinct by the problem definition?
-- Charles
Anton Ertl
>Anton Ertl wrote:
[Problem 189]
>> As far as I understand the task, rotation and reflection is a
>> non-issue. If you enumerate all colourings for one orientation,
>> that's the answer.
>>
>> If you had a way to create only one of, say, the two mirror images
>> (along some axis), then you would have to count it twice, unless the
>> colouring is symmetric. But with the usual enumerative approach we
>> get both mirror images anyway.
>
> From the problem statement:
>
> A colouring C' which is obtained from a colouring C by rotation or
>
> reflection is considered distinct from C unless the two are
> identical.
>
> How many distinct valid colourings are there for the above
> configuration?
>
>I read that as saying if a rotation or reflection is identical, it
>is not distinct, and should not be counted. Am I confused as to
>the interpretation?
>
>For example, with a much smaller triangle, I generate the following
>arrangements:
>
> R R
> R B G G B R
>
>They are both generated by my algorithm, and add 2 to the enumeration.
>Are they distinct by the problem definition?
As I understand it, they are.
cac
>> Are they distinct by the problem definition?
>
> As I understand it, they are.
>
I wrote the code some time ago, and I can't really say that I
can read it. I am going to rewrite it from scratch, and see if I
get the same answer, but I remain pretty sure that reflections
and rotations are relevant.
-- Charles
William James
wrote:
PARI/gp
{
default( primelimit, 50000000 );
top = 10^8;
count = 0;
forprime( pr=2, sqrtint( top - 1 ),
count += primepi( top \ pr ) - primepi( pr ) + 1 );
count
}
About 170ms on a slow, 800MHz laptop.
William James
wrote:
>
> The other one was the hyper exponentation:
> 1777 ^ 1777 ^ 1777 .. (1855 times). .. ^ 1777
PARI/gp:
hyper(n,e) = r=1; m=Mod(n,100 000 000);for(i=1,e,r=lift(m^r)); r
William James
Marcel Hendrix wrote:
> an...@mips.complang.tuwien.ac.at (Anton Ertl) writes Re: [SPOILER] Re: Euler problem #187
>
> > m...@iae.nl (Marcel Hendrix) writes:
> > >an...@mips.complang.tuwien.ac.at (Anton Ertl) writes Re: [SPOILER] Re: Euler problem #187
> [..]
> >>Your code, pasted into Gforth-fast, generates a number that is rejected by
> >>the Euler solution checker?
>
> > If I run it on a 64-bit system, the result is correct. If I run it on
> > a 32-bit system, or if I run it on iForth, it is too high by 1. It's
> > not obvious where the bug is.
>
> You have strong evidence that the 64-bit system (and the algorithm)
> are perfect? I cannot read the code with these funny loops.
Let it be known throughout the land of Zeedibul
That Hendrix finds Forth quite unreadable.
Bernd Paysan
>>
>> : PRIMES ( -- n )
>> FLAGS plimit 2/ 1 FILL
>> 3 EFLAG FLAGS DO
For speed, please run the outer loop only between FLAGS and
FLAGS+sqrt(plimit)/2. This is completely sufficient, as any non-prime is
already dividable by one number that's at most the square of plimit, so
after that, you'll never mark any new number as non-prime. This "small"
optimization improves the algorithm of sieve from O(n^1,5) to O(n), which
is a huge advantage if your n is 100 millions. The benchmark was never
optimized, since one idea of a benchmark is to run the same stupid thing,
and not change the rules.
> forprime( pr=2, sqrtint( top - 1 ),
Yes, that's the right end here.
--
Bernd Paysan
"If you want it done right, you have to do it yourself"
http://www.jwdt.com/~paysan/