[Jprogramming] An easy Euler Project one (485)

[Mike Day]

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 <rauldmiller at gmail.com>
>> To: Programming forum <programming at jsoftware.com>
>> 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