Wow! 31 is a sufficient maximum for the full Problem 485; 29 isn't. (Trial and error, not mathematical analysis!!!)
It takes about 27 seconds with 31 as maximum prime compared to the 45 seconds using all primes <: %: n . Cheers, Mike PS I played around with that "sliding window" article I cited earlier. Not surprisingly, it crawls like a snail with an effort to program it in J. Too loopy. Perhaps I could tacitise it, but not easy and hardly worth the effort. On 25/10/2014 17:51, 'Pascal Jasmin' via Programming wrote:
combined our approaches with a simple modification to dsieve add 2 lines to top of function: ( %: y) dsieve y : change the while line to: while. p <: x do. for 1e5 S 4e8, you do not need to look at primes higher than 13, and so this runs in 9 seconds timespacex '100000 ([: +/ >./\) 13 dsieve 40000000' 9.74692 2.51659e9 Smaller ranges (y) increase the required prime relative to max u (x)... 1e4 S 4e6 requires primes up to 23 10000 ([: +/ >./\) 23 dsieve 4000000 932689064 There is probably a mathematical relationship but a simple conjecture would be that if including the next prime in the sieve loop does not increase the S total, then no other prime above that would either. So there could be a built in loop short circuit. ________________________________ From: Mike Day <[email protected]> To: [email protected] Sent: Saturday, October 25, 2014 6:29 AM Subject: Re: [Jprogramming] An easy Euler Project one (485) 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
