Yes, I was more interested in the performance of "d" - as far as I
know, my
J method for the Euler problem doesn't need M. , as it only visits each
integer
once, at the expense of attempting to calculate all d values in one
array, but
never mind.
As for comparing numdiv() with J, it's tricky doing a straight comparison,
but this edited copy of a GP session might give a rough idea:
(NB. are my embedded comments)
<<
NB. define a user-function in Pari-GP to return a vector of d values
from n to n+k:
(09.02) nd(n,k)={vector(k,i,numdiv(n+i))}
NB. hit a memory limit:
(09:02) gp > #nd(1000000,1000000)
*** at top-level: #nd(1000000,1000000)
... NB. skip a bit
*** numdiv: the PARI stack overflows !
current stack size: 16000000 (15.259 Mbytes)
[hint] you can increase GP stack with allocatemem()
(09:02) gp > allocatemem()
*** Warning: new stack size = 32000000 (30.518 Mbytes).
NB. But Windows Task Manager shows gp.exe is using ca 17 MB.
(09:02) gp > #nd(1000000,1000000)
time = 1,231 ms.
%46 = 1000000
(09:02) gp > nd
%47 = (n,k)->vector(k,i,numdiv(n+i))
(09:02) gp > nd(n,k)={vector(k,i,numdiv(n+i))}
%48 = (n,k)->vector(k,i,numdiv(n+i))
(09:03) gp > #nd(1000000,1000000)
time = 1,227 ms. NB. ie 1.227 seconds
%49 = 1000000
(09:07) gp > #nd(1000000,100000)
time = 115 ms.
%50 = 100000
(09:08) gp > #nd(1000000,10000)
time = 12 ms.
NB. overhead of creating a vector with positive small integer values?
(09:23) gp > #vector(1000000,i,i+1000000)
time = 104 ms.
%53 = 1000000
>>
J time to evaluate powers of prime factors on every element in a vector;
presumably a J primitive for d would perform similarly. Roger will know
if there's a better comparator:
1 ts'#@{:@(__&q:)"0 q'[q=:1000000+i.1000000
4.85031 8.39565e6
1 ts'#@{:@(__&q:)"0 q'[q=:1000000+i.100000
0.46814 1.05562e6
1 ts'#@{:@(__&q:)"0 q'[q=:1000000+i.10000
0.0462873 138112
The following are closer to GP in time but with
a considerable memory hit:
1 ts'#q:"0 q'[q=:1000000+i.1000000
1.62014 8.05308e8
1 ts'#q:"0 q'[q=:1000000+i.100000
0.157885 8.38881e7
1 ts'#q:"0 q'[q=:1000000+i.10000
0.015095 1.04878e7
I would think Pari GP uses the same algorithm for
numdiv() for these fairly small (!!) integers,
namely product of 1+each prime's power. It's just
doing it in C binary whereas we're trying to use the
appropriate output from q: in interpreted J.
Pari GP of course uses very sophisticated methods
for "large" integer factoring, > 2^31 say, but
we're not in that area here, thank goodness.
JVERSION
Engine: j701/2011-01-10/11:25
Library: 8.03.05
Qt IDE: 1.2.4/5.3.2
Platform: Win 64
Installer: J803 install
InstallPath: c:/d/j803
My laptop has 6GB, at about 1.8 GHz.
Thanks again,
Mike
PS - I've forgotten how to access Roger's new prime stuff
in Dyalog APL...
On 22/10/2014 06:52, Roger Hui wrote:
I have only thought about d, the number of divisors, and not at all about
how it fits with the overall problem.
Another approach, if you are not applying d to a range of values, is to try
using M.
How much faster is Pari GP compare to J, on numdiv() versus J
implementations of d? I haven't researched finding number of divisors. If
there is a way to find it without factoring that may explain why Pari GP is
faster. Also, factoring in J isn't necessarily the best on large (>1e13,
say) arguments as I am far from an expert.
On Tue, Oct 21, 2014 at 4:33 PM, Mike Day <[email protected]>
wrote:
Yes, agreed - I started trying a sieve method earlier today.
So far it wasn't doing better than the prime factoring
approach using q: though it would no doubt benefit from
more careful programming.
As I've already got a solution I don't think I'll pursue
the sieve method further now. I would have tried
that, and also splitting the problem into smaller blocks.
Sure you're not tempted to add a primitive for d ? Why
is Pari GP quicker here? It's quicker even when my 6Gb
machine isn't thrashing for memory.
Pascal Jasmin's alternative for d is attractive. In tacit form
*/@:>:@:(#/.~)@:-.&0"1@:q:
is quicker and leaner than
*/@:(>:@{:)"2@(__&q:)
for >: i. 1000000 though slower albeit still leaner
for >: i. 10000000 ...
Cheers,
Mike
On 21/10/2014 21:47, Roger Hui wrote:
__ q: 30
2 3 5
1 1 1
The number of divisors are the possible vectors of exponents, which in
this
case are {0 1;0 1;0 1. For each prime, if the maximum exponent is e, then
the possible exponents for that prime is 0 1 2 3 ... e-1,e, i.e. the
number
of possible exponents is 1+e.
That's why I suggested using a sieve. Factoring is hard. But if you are
computing the number of divisors of a range, you can avoid factoring
altogether.
A few weeks ago I came across a neat (but inefficient) expression for the
divisors of n :
n=. 144
~. n +. 1+i.n
1 2 3 4 6 8 9 12 16 18 24 36 48 72 144
__ q: 144
2 3
4 2
*/ 1 + 4 2
15
# ~. n +. 1+i.n
15
On Tue, Oct 21, 2014 at 12:46 PM, 'Pascal Jasmin' via Programming <
[email protected]> wrote:
I don't completely understand number of divisors, but note:
*/@:(>:@{:)"2@(__&q:) 30
8
q: 30
2 3 5
I don't know why number of divisors of 39 would not be 3. Though I guess
the divisors are 1 2 3 5 6 10 15 30
So, your d function looks like it should be efficient as long as __ & q:
is comparable in speed to just q:
*/@:(>:@{:)"2@(__&q:) 105600
96
If not, there is this alternative basis for d:
*/ >: #/.~ q: 105600
96
----- Original Message -----
From: Mike Day <[email protected]>
To: [email protected]
Cc:
Sent: Tuesday, October 21, 2014 2:54 PM
Subject: Re: [Jprogramming] An easy Euler Project one (485)
No takers to the previous posting (below)?
0) Pari GP has a built-in function for number of divisors;
J doesn't, although it seems a nice candidate for an extension of p:
or q:
1) My fairly crude Pari GP function to solve this problem runs faster
than
the slightly more elegant although admittedly brute-force J verb(s).
The Pari GP solution for the actual Project Euler problem 485 took
just under 7 minutes, and is correct, so optimisation isn't really
required.
I've just received a J solution after 1 hour 20 minutes, also
correct.
Both run in interactive mode.
2) I did have a go at optimising the GP Pari, using a count array of
d-values in the current window (of size 100 000 in the actual problem);
that array can be considerably smaller than 100 000; I messed around
with forward and backward pointers but it was slower, and buggy.
3) Running a smaller problem under jpm suggests that the bottleneck
is the d (number of divisors) verb.
S 1e6 1e4 takes about 6.5 seconds of which 6 seconds is spent in "d".
4) So is there a better d verb out there - capable of running efficiently
on a vector argument? (Or a p: or q: extension as per 0 above?)
Thanks,
Mike
On 20/10/2014 15:26, Mike Day wrote:
NB. Maximum number of divisors
NB. Problem 485
NB. Published on Saturday, 18th October 2014, 04:00 pm
NB. Let d(n) be the number of divisors of n.
NB. Let M(n,k) be the maximum value of d(j) for n ≤ j ≤ n+k-1.
NB. Let S(u,k) be the sum of M(n,k) for 1 ≤ n ≤ u-k+1.
NB.
NB. You are given that S(1000,10)=17176.
NB.
NB. Find S(100 000 000,100 000).
NB. A brute-force J method:
d =: */@:(>:@{:)"2@(__&q:)
S =: +/@:(>./\ d@>:@i.)~/
timer'S 1000000 1000' NB. NOT the published problem!
+-----+---------+
|7.769|140671058|
+-----+---------+
timer'S 1000000 10000'
+-----+---------+
|7.147|175757800|
+-----+---------+
These are faster in a pedestrian Pari GP script,
but it does use the built-in function, "numdiv":
(08:10) gp > S(1000000,1000)
time = 3,302 ms.
%5 = 140671058
(15:13) gp > S(1000000,10000)
time = 2,889 ms.
%6 = 175757800
Curious.
Mike
PS The Pari script:
S(u,k)= {
rv=vector(k,i,numdiv(i)); \\ rolling vector of "d" values
t=vecmax(rv);curmx=t;
ipos=0;
for(n= k+1, u,
oldrvi=rv[ipos+1];
newrvi=numdiv(n);
rvi=rv[ipos+1]=newrvi;
if(curmx<rvi, curmx=rvi,
if(oldrvi==curmx,curmx=vecmax(rv));
);
t+=curmx;
ipos=(1+ipos)%k;
);
t;
}
---
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/8432 - Release Date: 10/22/14
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
