I'm glad the Perl module is working well. Pari/GP has been quite an inspiration, and of course has many more features and some things are much faster in Pari (e.g. partitions, partitions iterator). Pari, like Perl 6 and Python, has a seamless transition to >64-bit integers that Perl struggles with, and of course arbitrary precision floating point.

Fundamental changes that can cut down on the number of calls to sigma or reduce the argument would be the biggest win, albeit one that would apply to both systems. CRG mentions a few ideas.

I thought about doing something like the mod 12 setup to cut down on the loop overhead, since that's not an area Perl does well in, but at ~10^13 enough of the time is spent in sigma/divisor_sum that it didn't seem worth a huge amount. In fact the best way to solve that (for me) is to just add a new function to XS.xs since that cuts out even more overhead (no input validation on each call to sigma, no conversion of variables, etc.), something like:
[code]
void findcyc(IN UV a, IN UV e)
CODE:
{
UV b, d, i;
if (a & 1) a++;
for (b = a; b <= e; b += 2) {
if (b % 4 == 0 || b % 3 == 0) {
for (i = b, d = 1; d <= 4 && i > b; d++) {
i = divisor_sum(i,1) - i;
if (i == b) printf("found %lu %lu\n", b, d);
}
}
}
}[/code]Now it's all in C and the overhead is negligible. You could do the same with Pari. This works very well at 10^7, but does very little at 9*10^12 -- the overhead truly is in the factoring at this point. I found this method was nice when I was doing tabulations of pseudoprimes since the work per call there was quite small even with larger inputs, so any overhead reduction was welcome in tasks that would run for months.

parallel for loops (and parforprime) are a neat new feature in Pari. Doing that in Perl is possible but gets a bit ugly. This also just works for a single machine. My tactic for things like this is just running lots of copies at different ranges across the machines, and do manual book-keeping to track who is running what.

Thanks again for the help. The approach you mention is what I've been doing for the last few months. I have a manually invoked set of procedures that assign blocks and set up each machine to run one iteration per core/thread. With the change to Perl, I've had to increase the assignment block.:smile:

A question, though: What is the ",1" within the divisor_sum call. I couldn't find a reference to it and have been using just divisor_sum(i). Of course, I'm also just running in Perl with no C extensions, since I'm not familiar with that aspect (yet).

All of the changes I've tried also resulted in no noticeable speedup from my original conversion.

This portion of the thread should probably be moved out of the "PARI commands" thread. Of course, I'd hope to find it, if it did move.:smile:

P.S. The most work was getting ntheory into my Perl on all my machines. I'm not really familiar with modules and used cpan, but sometimes I had to go through the installation of Bundle::CPAN twice to get it to install correctly so ntheory would install.

Via PM, someone reminded me of the ,1 addition to the divisor_sum routine. I now remember seeing it in some documentation after all. It has also dawned on me that the reason I was using PARI, and now Perl, was to use their function to do the sigma/divisor_sum, so I didn't have to implement it on my own. However, as CRG pointed out, I should go ahead and do some of the work ahead of time and really, what I should be doing is studying the sigma/divisor_sum algorithm and just implementing a complete specialized C/C++ program. Other than the function itself, the rest of the coding is quite simple.

Now, off to find the algorithm...

Edit: ah, yes, I had seen the ,1 at the ntheory document on cpan...

The Perl module does all its real work in C (or C+GMP for bigints). It has pure Perl code for everything, but that's quite slow compared to the C code or Pari's C code. So it is possible to extract it all from Perl. I have various standalone portions (e.g. the LMO, Lehmer, etc. prime count code), but I haven't done it with factoring. I know someone had mentioned they were building the whole thing as a C library. Honestly I should turn it all into a plain C library, then have the Perl module use it.

[quote]Now, off to find the algorithm...[/quote]In Pari you can look for sumdiv in basemath/arith2.c, factoru in basemath/ifactor1.c. For ntheory, it's all in factor.c.

Thank you, everyone, for all the help you've provided. I'm quite confused at the moment - I guess the term would be overwhelmed. Unfortunately, I'll be busy with other things today and won't be able to even look at any of this.

I thought for a couple moments of simply grabbing factor.c and writing the trivial loops around it, but when I started looking at the includes, and their includes, etc., it cycles back to perl.h, so I'm not sure I can make that route work either. I suppose I need to learn how to do XS, next. I tried to look for the pari file, but didn't see it. This is probably because most of my machines are using repository pari/gp, although all have pari-gp2c installed. I also couldn't find, on line, any efficient algorithms for the sigma/divisor_sum type functions. I'm kind of leaning toward the Perl-XS side, but with the pari times being reported, I could easily swing back.

On the bright side, I have crossed into 10^13 sooner than expected and discovered a Y2K style flaw in my bash script that keeps track of progress across all the machines. It shouldn't take much to remedy, but I don't have the time right now.

Thanks, again, for all the help. I need to go study a bunch...

[QUOTE=EdH;386050] I also couldn't find, on line, any efficient algorithms for the sigma/divisor_sum type functions. I'm kind of leaning toward the Perl-XS side, but with the pari times being reported, I could easily swing back.
[/QUOTE]

I have a possibility in theory using sigma(b/4) ( for findcyc4) and setminus() ( for general use) but that might be just as slow though every second b for the findcyc4 can be divided by eight and add a new step. the pari times cmp wasn't in 2.5.3 so that could be a defeating thing ( though one can replace !cmp() with !bitxor() in that case). machines may be a major factor because for Danaj one thing takes over 3 minutes but on mine less than one. which means on his the codes I have posted ( which really haven't changed, except replacing break() with next()) would take about 51 seconds. also any parallel coding would need to be optimized to the machine if it was done. the reason I say sigma(b/4) for findcyc4 is because you can literally get the divisors of 4*x by finding whats new with divisors(x)*2 ( where setminus can come in) and then divisors(x)*4, this way the most you have to find divisors for is x.

I just thought I'd show everyone the versions of is and try in my code file:
[CODE]is(p)={
if(setintersect([p],vector(#Mevec,n,MeVec[n]^2))==[p],
my(q=sqrtint(p));
isprime(sum(k=0,q-1,binomial(q,k+1)*(2^q-1)^k))
);
(isprime(p) && factorint(2^p-1)[,2]==[1,1]~)
}[/CODE]

[CODE]try(lim=199) ={
forprime(p=7,lim,
if(setintersect([p],MeVec)==[p],
if(is(p^2),
print1(p",")
),
if(is(p),
print1(p",")
)
)
);
}[/CODE]

[code]setintersect([p],vector(#Mevec,n,MeVec[n]^2))==[p][/code]
is slow (and just looks wrong). Why would you intersect a set with a singleton? Either pre-compute
[code]sqSet = Set(vector(#Mevec,n,MeVec[n]^2))[/code]
and use
[code]setsearch(sqSet, p)[/code]
or compute terms one-by-one and exit early, like
[code]canFind(p)=for(i=1,#Mevec,if(Mevec[i]^2==p, return(1))); 0[/code]

Taking the first method you have
[code]sqSet = Set(vector(#Mevec,n,MeVec[n]^2));
is(p)={
if(setsearch(sqSet, p),
my(q=sqrtint(p));
isprime(sum(k=0,q-1,binomial(q,k+1)*(2^q-1)^k)) \\ ?
);
isprime(p) && factor(2^p-1)[,2]==[1,1]~
};[/code]
but this is still very strange, since you compute a big number, sum(k=0,q-1,binomial(q,k+1)*(2^q-1)^k), test whether it is prime, [i]and then throw away the result[/i]. Is that really what you wanted?

[QUOTE=CRGreathouse;386168]
Taking the first method you have
[code]sqSet = Set(vector(#Mevec,n,MeVec[n]^2));
is(p)={
if(setsearch(sqSet, p),
my(q=sqrtint(p));
isprime(sum(k=0,q-1,binomial(q,k+1)*(2^q-1)^k)) \\ ?
);
isprime(p) && factor(2^p-1)[,2]==[1,1]~
};[/code]
but this is still very strange, since you compute a big number, sum(k=0,q-1,binomial(q,k+1)*(2^q-1)^k), test whether it is prime, [i]and then throw away the result[/i]. Is that really what you wanted?[/QUOTE]

I see what you mean for the setsearch, but for the sum result it isn't used again because neither was what it replaced and finding out if the sum was prime seemed to take less time than the integer division it once held if I remember correct I've been trying to speed up EdH's code recently I think the factorint actually slowed it down though.

I just realized ( partly on my own and partly from danaj) my codes for EdH's problem has a major flaw. cmp(b,i) doesn't do the check it replaced it checks if b>i OR b<i or i<b OR i>b doh. though in other developments I still see ways of speeding up the codes because every multiple of a perfect or abundant number is abundant so findcyc 6 if I can get it to take the 0 mod 12 from findcyc4 we can eliminate the first round of checks if we can prove it divides by a perfect or abundant number.

[QUOTE=science_man_88;386245]I just realized ( partly on my own and partly from danaj) my codes for EdH's problem has a major flaw. cmp(b,i) doesn't do the check it replaced it checks if b>i OR b<i or i<b OR i>b doh. though in other developments I still see ways of speeding up the codes because every multiple of a perfect or abundant number is abundant so findcyc 6 if I can get it to take the 0 mod 12 from findcyc4 we can eliminate the first round of checks if we can prove it divides by a perfect or abundant number.[/QUOTE]
I'm not following a lot of what you're posting, but do have a query:

How does 0 mod 12 equate to (0 mod 4 || 0 mod 3)? (Or, am I too far removed from your code?)

Won't 0 mod 12 miss numbers like 9950971671152, which is 0 mod 4, 2 mod 3 and 8 mod 12?

Thanks for all the effort you' re putting in, though...