OK, Jon.
Ihad been fiddling with the mu function; my previous solution reported
here used an explicit sieve for each numerator. The mu function avoids
needing a sieve.
mu(1) is 1
mu(n) is 0 if n has one or more repeated prime factors
otherwise
mu(n) is (-1)^k where k is the number of (non-repeated) prime factors
So if we build a sieve on factors of 36, we only need to consider 2 and
3 and mark off 1s at 2 4 6 ..., and at 3 6 ...
If we added (rather than or-ed) the ones for six as well as for 2 and 3
the sieve would be
(0) 0 1 1 1 0 3 0 1 1 1 0 3 ...
But if we add mu(2 3 6) = _1 _1 1 instead, the sieve would correctly be
(0) 0 1 1 1 0 1 0 1 1 1 0 1 ...
But knowing these mus, we can just count up the numbers of ones and
minus ones they would generate in the sieve.
So we can easily count the number of non-coprimes and subtract it from
the size of the interval.
ppmu =: 3 : 0 NB.
q =. >:i.y NB.
t =. 0x NB.
lim =. y <. 2 NB. interval limits
for_i. }.q}.~(-<.-:y) do. NB.
NB. get interval limits, hi-1 and lo; half-open interval is (lo,hi-1]
lim =. y <. lim + 3 2 NB.
NB. simple factors * mu; eg 36 ==> _2 _3 6, but not 4 or 9 as they're
not square-free
muf =. }.,@>@(*/&.>/)@(<@(1 ,])/.~) -~.q:i NB.
ti =. (-/lim) + +/ -/ (**<.&| ) lim %/mufNB.count coprimes to i
t =. t + ti NB.
end. NB.
t NB.
)
timer'ppmu 12000'
+--------+-------+
|0.170166|7295372|
+--------+-------+
NB. Your f can be
((1 = +.) *. (1r3&< *. 1r2&>)@:%~) >:@:i.
NB. and I needed to add a rank here:
g =:+/@: f"0
timer'h >:i.12000'
+-------+-------+
|18.2639|7295372|
+-------+-------+
Part of the problem with this approach is that f applies +. to the whole
interval [0,i-1] although you only require that part in (i%3,i%2)
So if
hh =: +/@:((] ([ (1 +/@:= +.) ] #~ (1r3 < ] * 1r2 > ])@:%~) i.)"0)@:>:@i.
timer'hh 12000'
+-------+-------+
|4.62402|7295372|
+-------+-------+
Times for ppmu are approximately O(n^2)
If the problem were much larger, it would be necessary to avoid applying
all those q: and instead work up another sieve on the numerators,
probably looping over all primes up to n.
BTW, the Euler problems are accessible again, though for reference only
at the moment.
Mike
On 21/06/2014 17:18, Jon Hough wrote:
Not that it matters, but since I asked the question, I will show my solution.
coprime =. 1&=@:+.
NB. get the integers coprime to y, which when put into a fraction with y as
denominator lie between 1/3 and 1/2.
f=.(coprime *. (((1r3&<)*.(1r2&>))@:%~)) (>:@:i.)
NB. get the number
g =.+/@: f
h =. +/@:g
h >: i.12000
gives the answer:
7295372
This seems to be the same as others have. But since Project Euler is no more,
we'll never know if it is correct.Note: I haven't tested the speed but it seems
to be a bit clunky, taking about 3~4 seconds on my PC.
Regards,Jon
Date: Tue, 17 Jun 2014 19:35:04 +0100
From: [email protected]
To: [email protected]
Subject: Re: [Jprogramming] Project Euler 73
In private correspondence, "Pascal Jasmin" suggests I air these
functions in the forum.
Please note that the following is a spoiler for an Euler Project
solution method, so I'm placing listings some way down.
On my Samsung ultranotebook, that's not ultrafast: (sorry for any
line-wraps)
NB. Yesterday's modified version "pppp" of my old function "pp" used for
a correct Euler submission.
timespacex 'pppp 12000'
0.584461 278400
NB. Today's slightly more complicated version "ppppp"
timespacex 'ppppp 12000'
0.177553 268928
NB. Here's the slower, simpler one, "pppp"
NB. Sorry - it's got two nested loops!!! I'll put an NB. on each line
to help with line-wrapping problems
Both these methods take the 1/3 and 2/3 limits on fractions as given.
pppp =: 3 : 0 NB.
q =. >:i.y NB.
t =. 0 NB.
for_i. }.q}.~(-<.-:y) do. NB. all numerators from 2 to n%2
qp =. y#0 NB. initialise sieve
for_qi. ~. q: i do. NB. loop over nub of prime factors of
numerator
qp =. qp +. y$1{.~-qi NB. set sieve on all products of prime
factors of current numerator
end. NB.
t =. t + +/-. qp#~ (q < 3 * i) * q > 2 * i NB. sum over all
non-products within (1/2, 1/3)
end. NB.
t NB.
)
As I said yesterday, I think it's pretty trivial.
This following version, "ppppp", speeds up even more at the expense of a
bit more complexity, by being more frugal with the sieve.
ppppp =: 3 : 0 NB.
q =. >:i.y NB.
t =. 0 NB.
for_i. }.q}.~(-<.-:y) do. NB.
low =. >: 2 * i NB.
hi =. y <. 3 * i NB.
qp =. 0#~ nq =. 1 + hi - low NB. sieve in interval (1+2i, 3i-1)
for_qi. ~. q: i do. NB.
qp =. qp +. (low-1) }. (low+nq-1)$1{.~-qi NB. offsetting the
ones for prime products
end. NB.
t =. t + +/-. qp NB.
end. NB.
t NB.
)
pppp 12000
7295372
ppppp 12000
7295372
JVERSION
Engine: j701/2011-01-10/11:25
Library: 8.02.09
Qt IDE: 1.1.2/5.3.0
Platform: Win 64
Installer: J802 install
InstallPath: c:/d/j802
More recent Euler problems tend to have much larger N, and methods like
the above often only serve to explore the foothills (or perhaps
molehills) of those puzzles.
Have fun!
Mike
-----
No virus found in this message.
Checked by AVG - www.avg.com
Version: 2014.0.4592 / Virus Database: 3972/7694 - Release Date: 06/17/14
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
-----
No virus found in this message.
Checked by AVG - www.avg.com
Version: 2014.0.4592 / Virus Database: 3972/7717 - Release Date: 06/21/14
-----
No virus found in this message.
Checked by AVG - www.avg.com
Version: 2014.0.4592 / Virus Database: 3972/7718 - Release Date: 06/21/14
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
