FWIW, I've finally cracked the minute for the full problem
using the improved sieve below. It takes about 45 sec on
my laptop.
3 points:
0) Problem 485 solvers' discussion reveals this discussion about
maxima in a sliding window:
http://richardhartersworld.com/cri/2001/slidingmin.html
I doubt if it outperforms J's elegant method, eg
10 >./\ v
but I might have a look.
1) Several solvers seemed to have used a sieve, others brute-
force, as I did originally. Some discussion of the windowing,
including (0) above.
2) I'm puzzled why this method of producing an arithmetic
progression of indices into an array
i =. <:+/\p#~<.n%p
is quicker than the more obvious:
i =. I.n$1{.~ - p
Here's the improved sieve; I expect there will be gratuitous
line-throws again:
dsieve =: 3 : 0
NB. SLIGHTLY FASTER WITH m THAN s IN PREVIOUS VN
m =. 1#~n =. y NB. running product array
d =. m NB. number of divisors array
NB. every number i is a product of pij^qij for some distinct primes, pij
NB. required di is the product of (qij+1)
p =. 2 NB. seed first prime candidate
while. p <: %:n do.
np =. (<.n%p)<. <. p^npl =. <: >.p^.n
pl =. <. p^i.npl+1
NB. get q = powers of prime p in multiples of p
NB. eg q for p=3 is 1 1 2 1 1 2 1 1 3 for n < 27^2
NB. cf 3 6 9 12 15 18 21 24 27
q =. (}: , >:@{:)@:($~ (p * #))^:(0>.npl-1)1
q =. q$~<.n%p NB. extend q for all n(%p)
NB. i =. I.n$1{.~ - p NB. indices of all multiples of p
NB. next line, alternative for getting index to amend s and d,
NB. is FASTER!
i =. <:+/\p#~<.n%p
m =. (<. (i{m) * q{ pl) i } m NB. embed p^q for all multiples of p
d =. ( (i{d) * q + 1) i } d NB. adjust d with contributions from p^q
NB. first three sets of m & d
NB. +-+------------------------------------------------------------+
NB. |2|1 2 1 4 1 2 1 8 1 2 1 4 1 2 1 16 1 2 1 4 1 2 1 8 1 2 1 4 1 2|
NB. | |1 2 1 3 1 2 1 4 1 2 1 3 1 2 1 5 1 2 1 3 1 2 1 4 1 2 1 3 1 2|
NB. +-+------------------------------------------------------------+
NB. +-+----------------------------------------------------------------+
NB. |3|1 2 3 4 1 6 1 8 9 2 1 12 1 2 3 16 1 18 1 4 3 2 1 24 1 2 27 4 1 6|
NB. | |1 2 2 3 1 4 1 4 3 2 1 6 1 2 2 5 1 6 1 3 2 2 1 8 1 2 4 3 1 4|
NB. +-+----------------------------------------------------------------+
NB.
+-+---------------------------------------------------------------------+
NB. |5|1 2 3 4 5 6 1 8 9 10 1 12 1 2 15 16 1 18 1 20 3 2 1 24 25 2 27 4
1 30|
NB. | |1 2 2 3 2 4 1 4 3 4 1 6 1 2 4 5 1 6 1 6 2 2 1 8 3 2 4 3 1 8|
NB.
+-+---------------------------------------------------------------------+
NB. smoutput p;m,:d
NB. wd'msgs'
p =. 4 p: p NB. next prime
NB. next line, alternative for next prime, is too slow
NB. p =. ((1+p+1 i.~p&}.) ) m
end.
NB. final m is 1+i.n except for primes >: %: n
NB. double d for every such prime
d * >:(<>:@i.@# ) m
)
Perhaps I'll stop here, who knows!
Mike
On 23/10/2014 16:06, Mike Day wrote:
In case anyone's interested, here's a faster, if slightly
more opaque, version of a sieve for number of divisors
for all positive integers up to a limit.
It's better commented though.
Sorry for the inevitable double line-throws.
dsieve =: 3 : 0
s =. <.>:i.n=.y NB. original numbers to be "factored"
d =. n#1 NB. number of divisors array
NB. every number i is a product of pij^qij for some distinct primes, pij
NB. required di is the product of (qij+1)
p =. 2 NB. seed first prime candidate
while. p <: %:n do.
np =. (<.n%p)<.<. p^npl =. <: >.p^.n
pl =. <. p^i.npl+1
NB. get q = powers of prime p in multiples of p
NB. eg q for p=3 is 1 1 2 1 1 2 1 1 3 for n < 27^2
NB. cf 3 6 9 12 15 18 21 24 27
q =. (p&((}:,>:@{:)@(]$~(*#)))^:(0>.npl-1))1
q =. q$~<.n%p NB. extend q for all n(%p)
i =. I.n$1{.~ - p NB. indices of all multiples of p
NB. next line, alternative for getting index to amend s and d,
NB. is slower
NB. i =. p&*&.>:i.<.n%p
s =. (<. (i{s) % q{ pl) i } s NB. remove p^q from multiples of p
d =. ( (i{d) * q + 1) i } d NB. adjust d with contributions from p^q
p =. 4 p: p NB. next prime
NB. next line, alternative for next prime, is too slow
NB. p =. (>:@i.&1 @: >&1 ) s
NB. after p=2, s = 1 1 3 1 5 3 7 1 9 5 11 3 13 7 15 1 ...
NB. after p=2, d = 1 2 1 3 1 2 1 4 1 2 1 3 1 2 1 5 ...
NB. after p=3, s = 1 1 1 1 5 1 7 1 1 5 11 1 13 7 5 1 ...
NB. after p=3, d = 1 2 2 3 1 4 1 4 3 2 1 6 1 2 2 5 ...
NB. etc...
NB. smoutput p;d,:s
NB. wd'msgs'
end.
NB. final s is all ones except for primes >: %: n
NB. double d for every such prime
d * 1 2{~ s>1 NB. or d * >: s > 1
)
Ssieve =: (+/@:(>./\ dsieve )~)/ NB. solves problem 485 in about 5
minutes.
Ssieve 1000 10
17176
Mike
On 22/10/2014 11:48, Mike Day wrote:
Having floated this boat, here's a sieve method for
obtaining all numbers of divisors for integers from
1 to n:
dsieve =: 3 : 0
s =. <.>:i.n=.y
d =. n#1
p =. 2
while. p < %:n do.
pp =. 1
np =. (<.n%p)<.<.p^npl =. <.p^.n
q =. np#0
pl =. p^i.npl+1
NB. get q = powers of prime p
NB. eg q for p=3 is 1 1 2 1 1 2 1 1 3 ...
while. pp<n do.
q =. q+np$1{.~-pp
pp =. <.pp*p
end.
NB. q =. +/@:((+/) ($(1{.~-) )"0 ]) pl NB. slower than the explicit loop
q =. q$~<.n%p
i =. I.mask =. n$1{.~ - p
s =. (<. (mask#s) % q{pl) i } s
d =. ((mask#d) * q + 1) i } d
p =. 4 p: p NB. or could search for next non-1 in }.s
end.
d * >:s>1
)
timer'#dsieve 1000000'
+-----+-------+
|1.446|1000000|
+-----+-------+
Pari GP does this in about 1.2 seconds, so a reasonable
comparison.
Thanks,
Mike
On 22/10/2014 00:15, 'Pascal Jasmin' via Programming wrote:
a complication with the approach is:
*/ 2 2 2 3 3 5 5 7 11
138600
d */ 2 2 2 3 3 5 5 7 11
144
and so for the ranges starting at 38600 to 63320 there is a greater
maximum, and it also appears that for some multiples of 138600, the
number of divisors exceeds (for part of the range) of multiples of
83160
a number that simplifies the process is:
*/ 2 2 2 3 5 7 11
9240
as it allows needing to examine with d only multiples of it.
----- Original Message -----
From: Raul Miller <[email protected]>
To: Programming forum <[email protected]>
Cc:
Sent: Tuesday, October 21, 2014 6:41 PM
Subject: Re: [Jprogramming] An easy Euler Project one (485)
See also http://oeis.org/A002182
Also Roger Hui's point about a sieve is probably a good one since this
problem only needs you to consider prime factors less than 10000.
Still,
those 1229 prime factors multiplied by the 1e8 limit means we would
need
something on the order of several terabytes of memory for intermediate
results if we were to store the entire sieve in the obvious way. So
some
extra work would be needed, and I'm not sure how that would pan out.
---
This email is free from viruses and malware because avast! Antivirus protection
is active.
http://www.avast.com
-----
No virus found in this message.
Checked by AVG - www.avg.com
Version: 2014.0.4765 / Virus Database: 4040/8449 - Release Date: 10/25/14
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
