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 ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
