Mike, I've been following along this thread and have learned from it. Thanks for sharing. I stumbled upon this article last night that is also relevant:
http://www.jsoftware.com/jwiki/Doc/Articles/Play103 It's much slower but has some good commentary on the topic On Thu, Oct 23, 2014 at 11:06 AM, Mike Day <[email protected]> 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/8438 - Release Date: 10/23/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
