Thanks, Joe. Eugene's discussion is very elegant.
I could have preset d for low factors 2 3 5 ... Not sure if I'll try it for this problem. The sieve approach may also be used to derive prime factors and their powers for a list of integers, especially if it runs from 1 to a limiting n, by fairly trivial modifications to my code. The q vector already holds the powers, and each p can be saved somewhere. The structure of the resulting data, storing several p's and q's for each integer with or without that integer explicitly is the major problem unless n is huge; Roger's q: variants indicate some of the possibilities. I've used this for some Project Euler problems, not always with success (!), and frequently in Fortran which doesn't have built-in primes or factoring, unlike J or Pari GP. Cheers, Mike On 23/10/2014 16:18, Joe Bogner wrote:
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
