[QUOTE=CRGreathouse;347486]
Write programs in your favorite text editor and then just read the file into gp with \r filename. You can even edit the gprc.txt or .gprc file to automatically read this file when you start your session to make this more convenient.
[/QUOTE]

Yeah, I'm doing it quite often actually, but instead of creting the file and reading using PARI, I am just creating the code and pasting. Only the functions I'm making PARI read.

[QUOTE=CRGreathouse;347486]
See [url]http://math.crg4.com/pari/highlight.html[/url] for syntax highlighting if your text editor does not have this already. If you're on a friend's computer (or at school, work, etc. and can't install programs) you can get syntax highlighting online at pastebin.com.[/QUOTE]

I saw your site few days ago, it's very good from your part do that, but my favorite text editor isn't in your list of highlight syntax.

[QUOTE=CRGreathouse;347482]
You could write a function to do that.
[code]fracpart(x)=x - floor(x)[/code]
or even
[code]fracpart(x)=x - x\1[/code]
where \ means "divide and round down".
[/QUOTE]

I just saw PARI have the function "frac(x)", I feel stupid for asking that.
Using "frac(x)" I could do that:

[CODE]
a=1.2345;
f=frac(a)*10^4;
c=ceil(f)
[/CODE]

Using this logic I could kinda make what I was needing, but it won't be automatic, I'll have to set the value "10^n" manually. :/

[QUOTE=ismillo;347524]I saw your site few days ago, it's very good from your part do that, but my favorite text editor isn't in your list of highlight syntax.[/QUOTE]

What do you use?

In addition to the editors I mentioned explicitly, anything GtkSourceView-based or GeSHi-based should work as well, though you may need an update or download for it. Alternately you can use a different editor for GP than for your other work -- I did this for Ruby, where I found an editor that was particularly good for it even though I didn't prefer it overall.

Or just use your existing one without highlighting, it's not that big of a deal unless you work with it a lot.

[QUOTE=ismillo;347524]Using this logic I could kinda make what I was needing, but it won't be automatic, I'll have to set the value "10^n" manually. :/[/QUOTE]

So write a function that does it!

I don't quite understand what you need but this seems very easy to automate.

[QUOTE=CRGreathouse;347554]What do you use?

In addition to the editors I mentioned explicitly, anything GtkSourceView-based or GeSHi-based should work as well, though you may need an update or download for it. Alternately you can use a different editor for GP than for your other work -- I did this for Ruby, where I found an editor that was particularly good for it even though I didn't prefer it overall.

Or just use your existing one without highlighting, it's not that big of a deal unless you work with it a lot.

[/QUOTE]

I do use Sublime Text.
The truth is, my problem isn't highlight syntax, my problem is the lack of freedom in PARI interface. I do like very much see what I'm writing and testing them while I writing.


[QUOTE=CRGreathouse;347554]

So write a function that does it!

I don't quite understand what you need but this seems very easy to automate.
[/QUOTE]

I was trying to get the fractional part of the number in put 'em in an array. Just like "digits(n)", but instead of the number, only the fractional part.

Something like this:

[CODE]
gp>a=1.23456;
gp>c=digits(ceil(frac(a)*((#Str(a)-#Str(a))+10^10)))
%1 = [2,3,4,5,6,0,0,0,0,0,0]
[/CODE]

But when I get "length" of the number, returns the precision series.

It works awesome with numbers like [TEX]\pi[/TEX] and [TEX]\sqrt[]{2}[/TEX].

Somehow if I could remove the 0s from the last digits, but the 0s must be in a large sequence, otherwise it would break the number.

[QUOTE=ismillo;347555]I do use Sublime Text.[/QUOTE]

I don't own a copy. You might be able to make your own syntax highlighting file by following their instructions here:
[url]http://docs.sublimetext.info/en/latest/extensibility/syntaxdefs.html[/url]
if you like.

[QUOTE=ismillo;347555]my problem is the lack of freedom in PARI interface. I do like very much see what I'm writing and testing them while I writing.[/QUOTE]

Sorry to hear that. But I don't really understand -- I can see what I'm writing in gp and I very often iterate small functions directly in gp. Is that not working for you, or do I not understand what you want to see?

gp uses [url=http://en.wikipedia.org/wiki/GNU_readline]readline[/url] to provide these facilities and it works very well for me. It lets you edit previous lines, has tab completion (type "isps" and hit TAB to get "ispseudoprime()"), and so forth.

[QUOTE=ismillo;347555]I was trying to get the fractional part of the number in put 'em in an array. Just like "digits(n)", but instead of the number, only the fractional part.[/QUOTE]

Does this do what you want?
[code]fracpart(x)={
x-=x/1;
my(t=precision(x),s=log(x)\log(10),v=vector(t));
x/=10^s;
for(i=1,t,
v[i]=x\1;
x=10*(x-v[i])
);
v
};[/code]
I'm sure it could be made more efficient but the basic idea is sound.

[QUOTE=CRGreathouse;347561]I don't own a copy. You might be able to make your own syntax highlighting file by following their instructions here:
[url]http://docs.sublimetext.info/en/latest/extensibility/syntaxdefs.html[/url]
if you like.



Sorry to hear that. But I don't really understand -- I can see what I'm writing in gp and I very often iterate small functions directly in gp. Is that not working for you, or do I not understand what you want to see?

gp uses [url=http://en.wikipedia.org/wiki/GNU_readline]readline[/url] to provide these facilities and it works very well for me. It lets you edit previous lines, has tab completion (type "isps" and hit TAB to get "ispseudoprime()"), and so forth.
[/QUOTE]

I will try to make gp highlight for Sublime with I succeed I'll post here. :D

It's basically this but, it's ok soon or later I'll get used to it. :D

[QUOTE=CRGreathouse;347561]

Does this do what you want?
[code]fracpart(x)={
x-=x/1;
my(t=precision(x),s=log(x)\log(10),v=vector(t));
x/=10^s;
for(i=1,t,
v[i]=x\1;
x=10*(x-v[i])
);
v
};[/code]
I'm sure it could be made more efficient but the basic idea is sound.[/QUOTE]

Didn't work here, I get this error:

"*** log: domain error in log: argument = 0"

It's pointed to "s=log(x)".

Edit: I think the error is here "x-=x/1", it should be "x-=x\1;", I guess?
When I changed common slash to backslash it's gave no error, and the result is giving the digits normally, however the last digits is "0" or "9".

[QUOTE=ismillo;347564]I will try to make gp highlight for Sublime with I succeed I'll post here.[/QUOTE]

Cool, I look forward to it!

[QUOTE=CRGreathouse;347566]Cool, I look forward to it![/QUOTE]

reading a possibly old manual it seems there's so many different ways to get the same answer:


squaring x ( it seems)
[CODE]norm(x) \\ for integers at least
sqr(x)
x^2[/CODE]

multiply x by 2 (it seems)
[CODE]trace(x)
2*x
x<<1[/CODE]

how many that can you the same thing are there ?

[QUOTE=science_man_88;347727]reading a possibly old manual it seems there's so many different ways to get the same answer:


squaring x ( it seems)
[CODE]norm(x) \\ for integers at least
sqr(x)
x^2[/CODE]

multiply x by 2 (it seems)
[CODE]trace(x)
2*x
x<<1[/CODE]

how many that can you the same thing are there ?[/QUOTE]

Lots of ways. But be careful, these aren't quite equivalent. trace(I) == 0, for example.

So, about the highlight, I couldn't make it, since it's necessary JSON and I don't know nothing about it.

Now, about the fractional part to integer, I could make this:
[CODE]
r2i(a)={
my(b=Vec(Str(a)),c=b[#Str(ceil(a))+2..#b],d=eval(c));
forstep(x=#d,1,-1,if(d[x]!=0,h=x;break;));
e=vecextract(c,[1..h]);
x=eval(e);
v2i1d(x)
}
[/CODE]

v2i1d(x) is CRG function for vectors to integer with one digit each position.

[CODE]
v2i1d(a)=subst(Pol(a), 'x, 10)
[/CODE]

It has some problems, however, I think it will handle ok.

If the last part of your script is just removing trailing zeros, you can do that faster with

[code]removeTrailingZeros(n)={
n / 10^valuation(n, 10)
};[/code]

[QUOTE=CRGreathouse;348149]If the last part of your script is just removing trailing zeros, you can do that faster with

[code]removeTrailingZeros(n)={
n / 10^valuation(n, 10)
};[/code][/QUOTE]

This is awesome. Work like a charm. Thanks. :D

[CODE]
removeTrailingZeros(n)={
n / 10^valuation(n, 10)
};

v2i1d(a)=subst(Pol(a), 'x, 10);

r2i(a)={
my(b=Vec(Str(a)),c=b[#Str(ceil(a))+2..#b],d=eval(c));
removeTrailingZeros(v2i1d(d))
}.;
[/CODE]

factorisation
 

Hi, i just started learning pari and im trying to use store the prime factor of a value. im using a=factor(1599) but it stores only the results of the prime factor. can someone give me some guidance on how i could store any of the prime factor value to a variable?

Pari/gp command
 

hi, im new to pari/gp and i need some guidance on how to go about storing my factor value into a variable. i have for example 192384 and i need to get 1 of the prime factor of this value.

i tried storing it to a variable "a=factor(192384)" but it stores all the prime value to it. can anyone guide me how i could get 1 of the prime factor please.

I was trying to look up a list of mersenne factors somewhere else online and by Google search got to this code on rosettacode.org but there's some parts I don't understand can anyone help explain it ?


[CODE]factorMersenne(p)={
forstep(q=2*p+1,[COLOR="red"]sqrt(2)<<(p\2)[/COLOR],2*p,
[COLOR="Red"][1,0,0,0,0,0,1][q%8][/COLOR] && Mod(2, q)^p==1 && return(q)
);
1<<p-1
};[/CODE]

I don't understand why that's the limit, I get it as sqrt(2)^(p+1) so is that from some part of the p+1 test ? as to the other part that's red I tried q=10 by itself for this part of the code and I see it changes the answer from 2 to 0, but I can't even speculate on this part. edit: doh it's a call to the vector shown there should be an easier way.

sqrt(2)<<(p\2) = sqrt(2) * 2^((p-1)/2)which is a fancy way of saying sqrt(Mp). [a << b = a left shifted b bits]

[1,0,0,0,0,0,1][q%8] is a shorthand way of saying v=[1,0,0,0,0,0,1]; v[q%8]. Just a way to test if q%8 == 1 or 7

When you check for factors of a number, you only need to check up to square root of the number. The limit is an obfuscated way to write sqrt(2^p-1), or sqrt(1<<p-1), if you like, which might be a bit faster because of the internal pari representation. In fact, its speed depends on the number of decimals, and it could be much faster if you first set default(realprecision,5), get the root and set the realprecision back. The backslash is integer division. Anyhow, the speed of this calculus is not relevant, because it is only executed once in the beginning, and only takes few milliseconds.

The "vectored" part you can understand it better if you write:
v=[9,8,7,6];
a=v[2]

This will return 8. Now imagine you don't use v, but write directly the constant, as "a=[9,8,7,6][2]", this will still return 8. (indexing starts with 1)

The faster way to write "if x%8==1 or x%8==7" is what you see marked red in your text. There is no faster way, this involve a modulus and a selection (pointer). If you do "y=x%8, if y==1 or y==7", you may do the calculus only once, but still do the comparison two times at each loop.

So, this loop is very fast, but due to the fact that it tests all 2kp+1 candidates, it is still about 3-4 times slower than the "tf_420c" function which I posted here around, which tests only the candidates in the right classes (which stand a miniscule chance to be prime).

Even that function is much slower that other stuff like mfaktc :razz: because it does not do sieving, it does not run multithreaded, and pari is not compiled, but interpreted.

edit: Grrrr axn!

[QUOTE=LaurV;372033]When you check for factors of a number, you only need to check up to square root of the number. The limit is an obfuscated way to write sqrt(2^p-1), or sqrt(1<<p-1), if you like, which might be a bit faster because of the internal pari representation. In fact, its speed depends on the number of decimals, and it could be much faster if you first set default(realprecision,5), get the root and set the realprecision back. The backslash is integer division. Anyhow, the speed of this calculus is not relevant, because it is only executed once in the beginning, and only takes few milliseconds.

The "vectored" part you can understand it better if you write:
v=[9,8,7,6];
a=v[2]

This will return 8. Now imagine you don't use v, but write directly the constant, as "a=[9,8,7,6][2]", this will still return 8. (indexing starts with 1)

The faster way to write "if x%8==1 or x%8==7" is what you see marked red in your text. There is no faster way, this involve a modulus and a selection (pointer). If you do "y=x%8, if y==1 or y==7", you may do the calculus only once, but still do the comparison two times at each loop.

So, this loop is very fast, but due to the fact that it tests all 2kp+1 candidates, it is still about 3-4 times slower than the "tf_420c" function which I posted here around, which tests only the candidates in the right classes (which stand a miniscule chance to be prime).

Even that function is much slower that other stuff like mfaktc :razz: because it does not do sieving, it does not run multithreaded, and pari is not compiled, but interpreted.

edit: Grrrr axn![/QUOTE]

Do you know of anyone ever trying something like: [CODE]? a=[1,2,3,4,5];
forstep(x=1,100,v=eval(a),
print(x","a);
a=vector(#a,n,a[#a-n+1]+a[n%#a+1])
;v=eval(a)
)
1,[1, 2, 3, 4, 5]
2,[7, 7, 7, 7, 2]
16,[9, 14, 14, 9, 14]
39,[28, 23, 23, 28, 18]
67,[41, 51, 51, 41, 56][/CODE] ?

[QUOTE=vittoire;367091]hi, im new to pari/gp and i need some guidance on how to go about storing my factor value into a variable. i have for example 192384 and i need to get 1 of the prime factor of this value.

i tried storing it to a variable "a=factor(192384)" but it stores all the prime value to it. can anyone guide me how i could get 1 of the prime factor please.[/QUOTE]

where is it in the factorization matrix originally ? what column and what row ?

You could do factor(192384)[1,1] to get the smallest prime factor, factor(192384)[2,1] to get the second (if it exists), and so forth.

installed a newer version of pari
 

I just want people to know that if some of the codes I post elsewhere don't seem correct I now have both 2.5.3 and 2.7.0 and probably another. I love some of these new add-on's to existing functions etc.

yes I made another try at aligen that I would hope is at least more inclusive than previous versions:

[CODE]aligentry(n)={c=[];
forpart(x=n-1,
if(#x==#vecsort(x,,8) && numdiv(x[1]*x[#x])-2==#x && setminus(divisors(x[1]*x[#x]),Vec(x))==[1,x[1]*x[#x]],
c=concat(c,x[1]*x[#x])
)
,a=[2,n-1])
;c
}[/CODE]

science_man, I don't understand what that function is trying to do. But looking back through this thread using the search term "aligen" I found this exchange, and maybe a better answer to it will help you:

[QUOTE=science_man_88;284065]is there an easy way to create all the distinct partitions of a number n without the number 1 involved ( though I also have a way with 1 involved I think) ? because I'm trying to speed up my aligen function and it's based on partitions.[/QUOTE]

[QUOTE=CRGreathouse;291835]Here's some naive code:

[code]partNo1(n)={
select(v->vecmin(v)>1,partitions(n))
};[/QUOTE]

With newer versions of gp, there's now a faster way to do this:
[code]partNo1(n)={
partitions(n, [2, n])
};[/code]

If you wanted those partitions with only distinct parts, you can do
[code]distinctPartNo1(n)={
select(v->#v==#Set(v), partitions(n, [2, n]))
};[/code]
or better yet
[code]distinctPartNo1(n)={
my(d=sqrtint(2*n+2));
if(d*(d+1)/2 > 2*n+2, d--);
select(v->#v==#Set(v), partitions(n, [2, n], d))
};[/code]
which is faster and requires less memory.

[QUOTE=CRGreathouse;372745]science_man, I don't understand what that function is trying to do. But looking back through this thread using the search term "aligen" I found this exchange, and maybe a better answer to it will help you:





With newer versions of gp, there's now a faster way to do this:
[code]partNo1(n)={
partitions(n, [2, n])
};[/code]

If you wanted those partitions with only distinct parts, you can do
[code]distinctPartNo1(n)={
select(v->#v==#Set(v), partitions(n, [2, n]))
};[/code]
or better yet
[code]distinctPartNo1(n)={
my(d=sqrtint(2*n+2));
if(d*(d+1)/2 > 2*n+2, d--);
select(v->#v==#Set(v), partitions(n, [2, n], d))
};[/code]
which is faster and requires less memory.[/QUOTE]

then I'm confused why does forpart exists ? it goes over the partitions of it's first value, it then checks if that partition has all distinct members, then for those with the proper number of elements to be the proper divisors without 1 of it's elements at opposite ends, then and only then does it check it those are the true divisors (because that could take a long time I wanted to cut them down first). the "a=" part is the same as the end of your partitions call.

[QUOTE=science_man_88;372747]then I'm confused why does forpart exists ? it goes over the partitions of it's first value, it then checks if that partition has all distinct members, then for those with the proper number of elements to be the proper divisors without 1 of it's elements at opposite ends, then and only then does it check it those are the true divisors (because that could take a long time I wanted to cut them down first). the "a=" part is the same as the end of your partitions call.[/QUOTE]

You can certainly use forpart -- and it should be easy to adapt my code to that use.

[QUOTE=CRGreathouse;372786]You can certainly use forpart -- and it should be easy to adapt my code to that use.[/QUOTE]

adding:

[CODE]floor(solve(y=1,n-1,sum(z=1,floor(y),z)-(n-1)))[/CODE] into my code cut the time in half and looks to still be as accurate to the older one. is that the type of fix you were talking about ?

Something like that, yes. Actually I'm a bit surprised solve() doesn't choke on that, wouldn't it be better to use something like

solve(y=1,n-1, y*(y+1)/2-n+1)\1

to remove the flat parts. I guess it doesn't really matter -- this is outside the hot path so speed is nearly irrelevant.

[QUOTE=CRGreathouse;372793]Something like that, yes. Actually I'm a bit surprised solve() doesn't choke on that, wouldn't it be better to use something like

solve(y=1,n-1, y*(y+1)/2-n+1)\1

to remove the flat parts. I guess it doesn't really matter -- this is outside the hot path so speed is nearly irrelevant.[/QUOTE]

yeah yours becomes faster by 19,715 ms by aligentry(130) I guess a quicker way is to check if every member of x is divisible by the primes in x. if not it can't be a list of divisors. parvector and select may be able to do this eliminating some of the checks.

Your original code:
[code]aligentry0(n)={c=[];
forpart(x=n-1,
if(#x==#vecsort(x,,8) && numdiv(x[1]*x[#x])-2==#x && setminus(divisors(x[1]*x[#x]),Vec(x))==[1,x[1]*x[#x]],
c=concat(c,x[1]*x[#x])
)
,a=[2,n-1])
;c
};[/code]

Cleaned up with indentation, lexical scoping ([b]always[/b] use "my" or "local" for variables), Set, a list instead of concatenation (much faster), etc.:

[code]aligentry1(n)={
my(c=List());
forpart(x=n-1,
if(#x==#Set(x) && numdiv(x[1]*x[#x])-2==#x && setminus(divisors(x[1]*x[#x]),Vec(x))==[1,x[1]*x[#x]],
listput(c, x[1]*x[#x])
)
,
[2,n-2]
,
solve(y=1,n-1, y*(y+1)/2-n+1)\1
);
Vec(c)
};[/code]

Additional thought: you're looking for a partition x1 > ... > xn equal to the proper divisors (> 1) of x1 * xk. So a faster method would be to loop through the primes p which might be the least member, look for partitions of n-p with all parts at least q = nextprime(p+1) and at most q*(n-q)\nextprime(q+1), and then add on p 'by hand'. This trims off a lot of partitions to check, I think.

one thing I realized was that x=vecsort(x,,8) eliminated more than #x=#vecsort(x), because it also asked setisset(x) effectively which allows use to only use ones that are in order already.

[QUOTE=CRGreathouse;372876]

Additional thought: you're looking for a partition x1 > ... > xn equal to the proper divisors (> 1) of x1 * xk. So a faster method would be to loop through the primes p which might be the least member, look for partitions of n-p with all parts at least q = nextprime(p+1) and at most q*(n-q)\nextprime(q+1), and then add on p 'by hand'. This trims off a lot of partitions to check, I think.[/QUOTE]

I did have difficulty trying to get yours to work, actually the reason I'm doing n-1 instead of n is because 1 is already eliminated that way so if prime p gets eliminated I'd look for partitions of n-(p+1) I might just make an outer parforprime ( a parallel forprime) loop to do this after figuring out a limit , and technically if n-1 is prime [TEX](n-1)^2[/TEX] gets missed by your code.
edit: doh, what am I thinking I just realized something before I posted and forgot it already. [2,n-3] will catch all for primes above it as wells so unless I go down through the primes both the parforprime and others save none. edit2: nevermind I realized that I'm not looking up to n-1 again.

[QUOTE=science_man_88;372886]one thing I realized was that x=vecsort(x,,8) eliminated more than #x=#vecsort(x), because it also asked setisset(x) effectively which allows use to only use ones that are in order already.[/QUOTE]

#x == #vecsort(x) is true for any vector x, since vecsort doesn't remove any elements unless you give it flag 8.

Set(x) is basically a better version of vecsort(x,,8) if your version is new enough that Set([1])==[1] .

[QUOTE=science_man_88;372894]I did have difficulty trying to get yours to work, actually the reason I'm doing n-1 instead of n is because 1 is already eliminated that way so if prime p gets eliminated I'd look for partitions of n-(p+1)[/QUOTE]

Sounds good, that's what I figured.

[QUOTE=science_man_88;372894]parforprime ( a parallel forprime)[/QUOTE]

You should test this before using to make sure your version can actually make reasonable use of it. Try something like

parforprime(p=2,7, runs_for_10_seconds())

to see if it runs in ~10 seconds or ~40 seconds (assuming 4 cores, adjust to suit otherwise).

parforprime still didn't make sense and seems to be slower for simple things. on another note:

[CODE]p=11;a=binary(2^((p-1)/2-1)-1);s=14;for(x=1,#a,s=(Mod(s,(2^p-1))^2-2)^2-2);s[/CODE]

here's another way to do the lucas lehmer test I was trying to make it close to:

[url]http://www.mersenne.org/various/math.php#trial_factoring[/url] , anyways off to make my other code better I guess.

[QUOTE=science_man_88;372972]parforprime still didn't make sense and seems to be slower for simple things.[/QUOTE]

It's definitely slower for simple things -- lots of overhead. I wouldn't consider doing less than a million cycles with each iteration.

[QUOTE=CRGreathouse;373004]It's definitely slower for simple things -- lots of overhead. I wouldn't consider doing less than a million cycles with each iteration.[/QUOTE]

well there's that, but also once p=2 is done it covers everything else because every value in the limits [amin,amax] will fall between 2,n-3 the only thing this misses is when n-1 is prime. So doing it first doesn't make sense but even doing it in reverse order doesn't either without a check on if it's already in c because as p decreases the range [p,n-(p+1)] is larger than [nextprime(p),n-(nextprime(p)+1)], for some reason I'm getting your code to be slower than mine overall.

Which codes are you comparing? I think we've both posted several.

[QUOTE=CRGreathouse;373027]Which codes are you comparing? I think we've both posted several.[/QUOTE]

I guess I made alterations that made my code faster:

[url]http://mersenneforum.org/showpost.php?p=372876&postcount=2453[/url]

if I take out the solve part you have that my original didn't yours takes 378 ms if I add your solve in my code takes 129 ms so it's really only the solve I can find as faster. if I add in the ( I typed it wrong before) improvement of using the check x==vecsort(x,,8) I get it down to ~111 ms for mine. and even with your solve back into it plus an equivalent to my check by using Vec(x)==Set(x) yours doesn't go below 200 ms.

Since you didn't answer my question I won't bother doing my own analysis.

[QUOTE=CRGreathouse;373038]Since you didn't answer my question I won't bother doing my own analysis.[/QUOTE]

I had a link to the codes I altered from, and how altering them changed things.

polynomial form
 

I know that:
2x+1 single mersennes when x is the previous single mersenne
2x^2+4x+1 double mersennes when x is the previous double mersenne

I know the relation of the exponents of a triple mersenne as they are the double mersennes but does anyone know of a general form other than (2^y+1)*(2^x-1) for triple, quadruple, quintuple, etc. mersennes ? edit: sorry wrong thread for some reason I thought this was the theory of mersenne primes thread

does anyone know why they haven't merged Pol and Polrev ? I could see a newer version with {n=1} being a switch where n=1 makes it Pol and -1 makes it Polrev, though I guess that's because I'm trying to make my own version that does something like both.

[QUOTE=science_man_88;373699]does anyone know why they haven't merged Pol and Polrev ?[/QUOTE]

Why would they?

[QUOTE=CRGreathouse;373708]Why would they?[/QUOTE]

I was just thinking that it would reduce the number of functions, but I can't make anything that quick my best attempt involved parsum:

[CODE]Pols(t,v='x,{n=1})=if(n==1,parsum(w=1,#t-1,t[w]*v^(#t-w),t[#t]),parsum(w=2,#t,t[w]*v^(#t-(#t-w+1)),t[1]))[/CODE]

edit: which is about 1200 times as slow for a vector of length 10000

[code]parmatrix(m,n,{X},{Y},{expr=0})={
a=matrix(m,n,x,y,0);
parfor(x=1,n,
parvector(m,y,
a[x,y]=eval(expr)
)
)
;a
}[/code]

I was trying to make a parallel function for matrices ( in case it might ever be useful)

You forgot my().

[QUOTE=CRGreathouse;374024]You forgot my().[/QUOTE]

It doesn't fully work as written even with my, what am I missing?

I would do
[code]parmatrix(m,n,expr=(a,b)->0)={
my(a=matrix(m,n));
parfor(x=1,m,
for(y=1,n,
a[x,y]=expr(x,y)
)
);
a
};[/code]

with the understanding that the third argument needs to be given as a closure:

[code]f(x,y)=x^2+sigma(y);
parmatrix(5, 5, f) \\ user function called by name
parmatrix(4, 6, gcd) \\ built-in function called by name
parmatrix(4, 4, (x,y) -> x+y) \\ anonymous function[/code]

You can't define a function in gp with a true expr argument (you'd need to use the PARI library directly in C).

[QUOTE=CRGreathouse;374081]I would do
[code]parmatrix(m,n,expr=(a,b)->0)={
my(a=matrix(m,n));
parfor(x=1,m,
for(y=1,n,
a[x,y]=expr(x,y)
)
);
a
};[/code]

with the understanding that the third argument needs to be given as a closure:

[code]f(x,y)=x^2+sigma(y);
parmatrix(5, 5, f) \\ user function called by name
parmatrix(4, 6, gcd) \\ built-in function called by name
parmatrix(4, 4, (x,y) -> x+y) \\ anonymous function[/code]

You can't define a function in gp with a true expr argument (you'd need to use the PARI library directly in C).[/QUOTE]

Okay, Thanks for the code and hint. One further question what function won't return all 0's from the matrix default of 0's ?

[QUOTE=science_man_88;374086]One further question what function won't return all 0's from the matrix default of 0's ?[/QUOTE]

I don't know what you mean. I gave three examples of functions that don't return all 0s.

[QUOTE=CRGreathouse;374090]I don't know what you mean. I gave three examples of functions that don't return all 0s.[/QUOTE]

I tried every one of those and got back 0's:

[CODE](11:48) gp > parmatrix(m,n,expr=(a,b)->0)={
my(a=matrix(m,n));
parfor(x=1,m,
for(y=1,n,
a[x,y]=expr(x,y)
)
);
a
};
(12:12) gp > f(x,y)=x^2+sigma(y);
(12:13) gp > parmatrix(5, 5, f)
%67 =
[0 0 0 0 0]

[0 0 0 0 0]

[0 0 0 0 0]

[0 0 0 0 0]

[0 0 0 0 0]

(12:13) gp > parmatrix(4, 6, gcd)
%68 =
[0 0 0 0 0 0]

[0 0 0 0 0 0]

[0 0 0 0 0 0]

[0 0 0 0 0 0]

(12:13) gp > parmatrix(4, 4, (x,y) -> x+y)
%69 =
[0 0 0 0]

[0 0 0 0]

[0 0 0 0]

[0 0 0 0][/CODE]

Would you paste in your startup message, like
[code] GP/PARI CALCULATOR Version 2.8.0 (development 16370-23a2ae3)
i686 running mingw (ix86/GMP-5.1.3 kernel) 32-bit version
compiled: Apr 14 2014, gcc version 4.8.2 (GCC)
threading engine: single
(readline v6.2 enabled, extended help enabled)

Copyright (C) 2000-2014 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes
WITHOUT ANY WARRANTY WHATSOEVER.[/code]
? Also try this not-actually-parallel version:
[code]parmatrix(m,n,expr=(a,b)->0)={
my(a=matrix(m,n));
for(x=1,m,
for(y=1,n,
a[x,y]=expr(x,y)
)
);
a
};[/code]

[QUOTE=CRGreathouse;374095]Would you paste in your startup message, like
[code] GP/PARI CALCULATOR Version 2.8.0 (development 16370-23a2ae3)
i686 running mingw (ix86/GMP-5.1.3 kernel) 32-bit version
compiled: Apr 14 2014, gcc version 4.8.2 (GCC)
threading engine: single
(readline v6.2 enabled, extended help enabled)

Copyright (C) 2000-2014 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes
WITHOUT ANY WARRANTY WHATSOEVER.[/code]
? Also try this not-actually-parallel version:
[code]parmatrix(m,n,expr=(a,b)->0)={
my(a=matrix(m,n));
for(x=1,m,
for(y=1,n,
a[x,y]=expr(x,y)
)
);
a
};[/code][/QUOTE]

[QUOTE]Reading GPRC: gprc.txt ...Done.

GP/PARI CALCULATOR Version 2.7.0 (released)
i686 running mingw (ix86/GMP-5.1.3 kernel) 32-bit version
compiled: Mar 23 2014, gcc version 4.6.3 (GCC)
threading engine: single
(readline v6.2 enabled, extended help enabled)

Copyright (C) 2000-2014 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY WHATSOEVER.[/QUOTE]

[CODE](12:44) gp > parmatrix(m,n,expr=(a,b)->0)={
my(a=matrix(m,n));
for(x=1,m,
for(y=1,n,
a[x,y]=expr(x,y)
)
);
a
};
(14:03) gp > parmatrix(10,10,f(x,y)=x^2+sigma(y))
%7 =
[ 2 4 5 8 7 13 9 16 14 19]

[ 5 7 8 11 10 16 12 19 17 22]

[ 10 12 13 16 15 21 17 24 22 27]

[ 17 19 20 23 22 28 24 31 29 34]

[ 26 28 29 32 31 37 33 40 38 43]

[ 37 39 40 43 42 48 44 51 49 54]

[ 50 52 53 56 55 61 57 64 62 67]

[ 65 67 68 71 70 76 72 79 77 82]

[ 82 84 85 88 87 93 89 96 94 99]

[101 103 104 107 106 112 108 115 113 118]

(14:04) gp > parmatrix(m,n,expr=(a,b)->0)=my(a=matrix(m,n));parfor(x=1,m,for(y=1,n,a[x,y]=expr(x,y)));a;
(14:04) gp > parmatrix(10,10,f(x,y)=x^2+sigma(y))
%9 =
[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

(14:04) gp > parmatrix(m,n,expr=(a,b)->0)=my(a=matrix(m,n));for(x=1,m,parvector(n,y,a[x,y]=expr(x,y)));a;
(14:05) gp > parmatrix(10,10,f(x,y)=x^2+sigma(y))
%11 =
[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]

[0 0 0 0 0 0 0 0 0 0]
[/CODE]

it appears the parallel loops are the problem.

Actually, on second thought, this might be a problem with the threading model -- the other threads can't change the global matrix a. Hmm. I don't know of an easy solution, maybe there is something that can be done with inline().

The 'right' solution would require working in PARI. You'd generate each column in parallel, then combine them with something like mkmat.

[QUOTE=CRGreathouse;374097]Actually, on second thought, this might be a problem with the threading model -- the other threads can't change the global matrix a. Hmm. I don't know of an easy solution, maybe there is something that can be done with inline().

The 'right' solution would require working in PARI. You'd generate each column in parallel, then combine them with something like mkmat.[/QUOTE]

[CODE]parmatrix(m,n,expr=(a,b)->0)={
my(a=matrix(m,n));
for(x=1,m,
a[x,]=parvector(n,y,eval(expr(x,y)))
);a;[/CODE]

one step closer it appears to fall flat with both parallel though. I even tried to turn it into a vector and do two parvector(); edit: if you want to alter the reading to column's I could make a flag that if 1 could go by column and if not go by row.

I tried the vector one again but with only one parallel function ( though I wonder if breaking parvector into pareval and vector could work):

[CODE]parmatrix1(m,n,expr(a,b)->0)
=my(a=Vec(matrix(m,n)));
for(x=1,n,
a[x]=parvector(m,y,eval(expr(x,y)))
);
a=Mat(a)[/CODE]

I'm thinking of making a check to see which is greater and working it that way the problem is this rotates the values of m and n if n is greater than m because doing the greater of the two in parvector got me down to around 8200 ms for one of the things I tried in expr ( admittedly in a 1000 by 1100 matrix) but I'd also have to change the order in which m and n are declared in a to not mess things up and this rotates the matrix.

I've been reading up on set's etc. ( partly brought on my the fog comment by RDS) and I'm wondering why forpart() doesn't allow the variable to go over the partitions of a set (aka the elements of the power set of the set x=[1,2,3,4,5,6,7,8,9,10,11,12] for example). I'm guessing there's probably already a Set based function like this ?

I notice also that you have "threading model: single" which means that you can't execute code in parallel anyway.

[QUOTE=science_man_88;374200]I've been reading up on set's etc. ( partly brought on my the fog comment by RDS) and I'm wondering why forpart() doesn't allow the variable to go over the partitions of a set (aka the elements of the power set of the set x=[1,2,3,4,5,6,7,8,9,10,11,12] for example). I'm guessing there's probably already a Set based function like this ?[/QUOTE]

If you want to loop over partitions, use
[code]for(i=0,n!-1, p = numtopart(n, i); print(p))[/code]
where you replace prime(p) with useful code and n with the number of items in the partition.

If you want to loop over subsets of a set, use
[code]for(i=0,2^#v-1, p = vecextract(v,i); print(p))[/code]
where v is the set and print(p) is replaced by useful code.

[QUOTE=CRGreathouse;374375]
If you want to loop over subsets of a set, use
[code]for(i=0,2^#v-1, p = vecextract(v,i); print(p))[/code]
where v is the set and print(p) is replaced by useful code.[/QUOTE]

once I saw them in order I thought it would cut a check out of aligentry0 but for n=73 it increases the number of vectors to checks to >19 * numbpart(73)^3 when the whole before took numbpart(73)+14.

maybe make a version that doesn't crash gp.exe
 

I made a recent edit on the OEIS but I know a further formula from it that leads to the Fermat quotients however my codes ( based on quick ways of setting things up that I've read in the OEIS) crashes at b(1152), I was hoping someone could suggest an edit that wouldn't crash it:

[CODE]a(n)=if(n,4*a(n-1)+1,0)[/CODE]// calculates A002450
[CODE]b(n)=p=prime(n);(3*a((p-1)/2))/p[/CODE] // calculates the Fermat quotients assuming p is calculable ( yes I know that form of finding it has a limit)

[QUOTE=science_man_88;384315]I was hoping someone could suggest an edit that wouldn't crash it:

[CODE]a(n)=if(n,4*a(n-1)+1,0)[/CODE]// calculates A002450
[CODE]b(n)=p=prime(n);(3*a((p-1)/2))/p[/CODE] // calculates the Fermat quotients assuming p is calculable ( yes I know that form of finding it has a limit)[/QUOTE]

You want to know why a(4655) fails. Does (4^4655 - 1)/3 fail also?

[QUOTE=CRGreathouse;384330] Does (4^4655 - 1)/3 fail also?[/QUOTE]

no. edit: but for me it fails at a(4651)

[QUOTE=science_man_88;384344]no. edit: but for me it fails at a(4651)[/QUOTE]

So the problem is with your function a(). You'll need to rewrite it to not fill the stack.

try at mersenne semiprimes
 

[CODE] = ()->for(k=0,2^199,if((isprime(floor(log(6*k+2)/log(2)))||(issquare(floor(log(6*k+2)/log(2)))&&isprimepower(floor(log(6*k+2)/log(2)))==2
&&isprime(2^floor(sqrt(log(6*k+2)/log(2)))-1)))&&#factor(6*k+1)[,1]==2&&sum(x=1,#factor(6*k+1)[,2],factor(6*k+1)[x,2])==2,print1(floor(log(6*k+2)/log(2))","));k=4*k)[/CODE]

I tried different versions like this before just going back through my pari history but I can't get much improvement if at all in fact some things I thought could be improvements don't appear to be like pre-calculating [CODE]floor(log(6*k+2)/log(2))[/CODE] seems to not really help much and I have been debating if I should use issquarefree(x) or core(x)==x because one seems faster above a certain number once I take the limit away, I even tried k as the exponent and a forstep and nothing seems to improve it even though I know there are checks I could do to improve it more like the squarefree scenario ( of course I can limit mod 120 etc first to lower the amount of squarefree checks when checking the prime part of the OR) any other ideas that might speed it up

Since your function is an unreadable single line I won't dissect it, but I notice that you're using floor(log(6*k+2)/log(2)) a lot. Let's define t = floor(log(6*k+2)/log(2)), which you'll notice stays the same for long streches of k and then goes up by 1. Why don't you loop over values of t directly, and then re-create the associated values of k if t passes your tests?

[QUOTE=CRGreathouse;385523]Since your function is an unreadable single line I won't dissect it, but I notice that you're using floor(log(6*k+2)/log(2)) a lot. Let's define t = floor(log(6*k+2)/log(2)), which you'll notice stays the same for long streches of k and then goes up by 1. Why don't you loop over values of t directly, and then re-create the associated values of k if t passes your tests?[/QUOTE]

each loop through the for loop k is multiplied by 4 and then adds one more. I'm looping through k values that produce odd exponent mersenne numbers I took care of it in my newest version it still doesn't trim a lot of time like i thought. I replaced it with p=logint(6*k+2,2) for each k used and more but I realize I used %39 to call it so in my new session I can't call it back to the quickest form I got to:

[CODE]try() ={
n=0;
m=vector(120,n,gcd(n-1,120)==1);
for(k=0,2^199,
p=logint(6*k+2,2);
if(
(
(
m[(p%120)+1] && core(p)==p && isprime(p)
)
||
(
issquare(p) && isprimepower(p,&n)==2 && isprime(2^n-1)
)
)
&& #factor(6*k+1)[,1]==2 && sum(x=1,2,factor(6*k+1)[x,2])==2,
print1(p",")
)
; k=4*k
)}[/CODE]

is the best I could do at rewriting my best implementation.

OK, that gives me a lot more to go on.

First, you should [i]always[/i] declare scope for your variables. In about 100% of cases this should be done with my(), with some residual cases using local() instead. Otherwise you overwrite global variables which is BAD and will cause unintended side-effects down the road.
[code]try()={
my(m=vector(120,n,gcd(n-1,120)==1), q);
for(k=0,2^199,
my(p=logint(6*k+2,2));
if(
(
(m[(p%120)+1] && core(p)==p && isprime(p))
||
(issquare(p) && isprimepower(p,&q)==2 && isprime(2^q-1))
) && #factor(6*k+1)[,1]==2 && sum(x=1,2,factor(6*k+1)[x,2])==2,
print1(p",")
);
k=4*k
)
};[/code]

Second you should split out your handling of the test of p from the rest of the program. This is just good general form. A first attempt would be a straight split:
[code]m=vector(120,n,gcd(n-1,120)==1); \\ global

is(p,k)={
my(q);
(
(m[(p%120)+1] && core(p)==p && isprime(p))
||
(issquare(p) && isprimepower(p,&q)==2 && isprime(2^q-1))
) && #factor(6*k+1)[,1]==2 && sum(x=1,2,factor(6*k+1)[x,2])==2
};

try()={
for(k=0,2^199,
my(p=logint(6*k+2,2));
if(is(p, k),
print1(p",")
);
k=4*k
)
};[/code]

Let's work on is() first. There's no need for m, it would be better to do the gcd each time. Much faster than core(n)==n is the function issquarefree(n), but you don't even need this at all because primes are always squarefree. That gives:
[code]is(p,k)={
my(q);
(
(gcd(p,120)==1 && isprime(p))
||
(issquare(p, &q) && isprime(q) && isprime(2^q-1))
) && #factor(6*k+1)[,1]==2 && sum(x=1,2,factor(6*k+1)[x,2])==2
};

try()={
for(k=0,2^199,
my(p=logint(6*k+2,2));
if(is(p, k),
print1(p",")
);
k=4*k
)
};[/code]

Next I notice that you're factoring 6k+1 three times, which seems extreme. One way to handle this would be to store the results of factoring, but in this case there's an easier solution: just look at the exponent vector. Numbers relatively prime to 120 are just those relatively prime to 30, but we can simplify further since that path is just looking for primes so we can just test that p > 5. In fact that test can be pulled outside (since the only prime square <= 5 is 4 which doesn't pass). Oh, and we need a definition of logint, so here we go:
[code]logint(N,b)=my(k);while(N>=b,N\=b;k++);k;

is(p,k)={
my(q);
p>5 &&
(isprime(p) || (issquare(p, &q) && isprime(q) && isprime(2^q-1)))
&& factor(6*k+1)[,2]==[1,1]~
};

try()={
for(k=0,2^199,
my(p=logint(6*k+2,2));
if(is(p, k),
print1(p",")
);
k=4*k
)
};[/code]
Now let's look at try(). I don't like abusing for() loops like this, a while() loop would make the transformation k -> 4k+1 clearer. But I think that it would be even clearer if you just wrote (4^n - 1)/3, and since the calculations inside is() take far longer than this calculation I'd just do that directly. (k <= 2^199 corresponds to n <= 100 here.) That also makes it clear that p = 2n+1 which was not at all clear before! Also we can remove our new function logint() since it isn't used.
[code]is(p,k)={
my(n);
p>5 &&
(isprime(p) || (issquare(p, &n) && isprime(n) && isprime(2^n-1)))
&& factor(6*k+1)[,2]==[1,1]~
};

try()={
for(n=0,100,
my(k=(4^n-1)/3, p=2*n+1);
if(is(p, k),
print1(p",")
)
)
};[/code]
Now the code is much clearer, and you can see that you're just looking for semiprimes of the form 2^m - 1 subject to some conditions on m -- probably necessary conditions for semiprimality. But we can go further! You're ultimately testing 2^p - 1 for primality, and if p is the square of a prime q then the factorization of 2^p - 1 must be 2^q-1 times another prime, so we don't need to factor 2^p-1 in generality in that case. This also lets us clean up the code to remove k, which is no longer needed.
[code]is(p)={
my(q);
if(p < 7, return(0));
if(issquare(p, &q) && isprime(q) && isprime(2^q-1)),
return(isprime((2^p-1)/(2^q-1)))
);
isprime(p) && factor(2^p-1)[,2]==[1,1]~
};

try()={
for(n=0,100,
my(p=2*n+1);
if(is(p),
print1(p",")
)
)
};[/code]
Evidently we can clean up the loop in try() to go directly over the odd numbers, no need for n. Also, since we're requiring it to be at least 7, let's just start the loop there and remove that condition from is().
[code]is(p)={
my(q);
if(issquare(p, &q) && isprime(q) && isprime(2^q-1),
return(isprime((2^p-1)/(2^q-1)))
);
isprime(p) && factor(2^p-1)[,2]==[1,1]~
};

try()={
forstep(p=7,201,2,
if(is(p), print1(p","))
)
};[/code]
At this point I'm done simplifying, but you could improve this further. For example, you could replace the test isprime(2^q-1) with a function LL(q) which uses the Lucas-Lehmer test to see if 2^q-1 is prime, and you could write trial division code that checks for small prime factors of 2^p-1 (using only the ones which could be factors, of course) -- if you find one you don't have to complete the factorization, you can just test the cofactor for primality.

Edit: Another good step would be to replace try() with try(lim=201) and using forstep(p=7,lim,2 so you can test different regions without editing the function.

sometimes I think you know too much about PARI in the sense that my all inbuilt function using codes are half as fast as the ones you make. I was doing the other testing to decrease the number of checks but I guess my checks took so much time that the decrease didn't matter. I thought about using k as the exponent but could never get it to work. my code caught 9 but not 4 and yours appears to forget about 9. I'll see what I can implement of your further suggestions.

I get
[code]9,11,23,37,41,49,59,67,83,97,101,103,109,131,137,139,149,167,197,199,[/code]
in 7 seconds. This could be improved, no doubt.

Edit: I did have to remove a superfluous paren from my code, try it now to see if it works.

[QUOTE=CRGreathouse;385545]I get
[code]9,11,23,37,41,49,59,67,83,97,101,103,109,131,137,139,149,167,197,199,[/code]
in 7 seconds. This could be improved, no doubt.

Edit: I did have to remove a superfluous paren from my code, try it now to see if it works.[/QUOTE]

[CODE](16:51) gp > is(p)=my(q);if(issquare(p,&q)&&isprime(q)&&isprime(2^q-1),(isprime((2^p-1)/(2^q-1))));isprime(p)&&factor(2^p-1)[,2]==[1,1]~;
(17:57) gp > try(lim=201)={forstep(p=7,lim,2,if(is(p), print1(p",")))};
(17:57) gp > try()
[B]11,23,37,41,59,67,83,97,101,103,109,131,137,139,149,167,197,199,[/B]
(17:58) gp > ##
*** last result computed in 3,970 ms.[/CODE]

oh and I tried to make the LL check but it made it take longer probably because I suck I was thinking of just using your M code ( if I can find it) or a vector of known to be mersenne prime exponents because even with a complete list I could use the code in theory to check up to 924309391636849 or a little higher.

The ultimate optimization would be something like
[url=http://factordb.com/index.php?query=2^n-1&use=n&n=1&VP=on&OD=on&FF=on&PRP=on&U=on&C=on&perpage=200&format=1&sent=Show]2^p-1[/url]
[url=http://factordb.com/index.php?query=2^%28n^2%29-1&use=n&n=1&VP=on&OD=on&FF=on&PRP=on&U=on&C=on&perpage=200&format=1&sent=Show]2^(p^2)-1[/url]

[QUOTE=CRGreathouse;385549]The ultimate optimization would be something like
[url=http://factordb.com/index.php?query=2^n-1&use=n&n=1&VP=on&OD=on&FF=on&PRP=on&U=on&C=on&perpage=200&format=1&sent=Show]2^p-1[/url]
[url=http://factordb.com/index.php?query=2^%28n^2%29-1&use=n&n=1&VP=on&OD=on&FF=on&PRP=on&U=on&C=on&perpage=200&format=1&sent=Show]2^(p^2)-1[/url][/QUOTE]

I'm trying now to optimize to check if the division is prime from the formula I got for when the exponent's a certain multiple plus something and the reduced binomial theorem on [URL="http://en.wikipedia.org/wiki/Binomial_theorem#Statement_of_the_theorem"]wikipedia [/URL] for when you have a variable plus 1 to a power I get that the result of division is:
[CODE]my(n=q);sum(k=0,n-1,binomial(n,k+1)*y^k)[/CODE]

but I'm trying to think the best way to use this to test primality.

PARI Routine for Finding Cycles
 

I've been running a compiled version of the following routine on several machines in my quest for a 4-cycle sequence. I am running a compiled version, more due to the ease of calling it via a bash script, than because of any speed difference. There is only a slight difference. I have tried many methods of speeding this up further, to include handling the counter outside the main function and different types of loops, etc. As can be seen, I am skipping many numbers in my search, based on information from another member on what they searched. All of this is in a different thread. On my fastest machine, a 1M block (in 9*10^12) is covered in about 44 seconds. The routine does return 2-cycles and found the one known 4-cycle that was within its search pattern.

Finally, on to my questions:

1. Is there room for any appreciable speed improvement in the routine?
2. Is the sigma() function something I could use with a GPU?
3. If 2, then how would I go about (learning and) modifying it to make use of a GPU?
[code]
findcyc(a,e)=
{
if(a%2==1,a+=1);
print("version 20140530b, running ", a, " through ", e);
gettime();
forstep(b=a,e,2,
if(b%4==0 || b%3==0,
i=b;
for(d=1,4,i=sigma(i)-i;
if(i<b,
break(),
if(i==b,
write("../../Desktop/cyclesFound", b, " ", d);
break()
)
)
);
if(b%1000000==0,
print(b,"/",e," ",gettime()/1000+0.," s/M");
write("../../Desktop/cycleStatus", b," ",e)
)
)
)
}
[/code]Thanks for any/all assistance...

Essentially all the work is done in sigma(), which spends essentially all of its time factoring the number. The best thing you could do would be to set up a sieve to factor the first instance of i more-or-less 'for free'. Next it would probably be worthwhile to write a customized factor function that could cut out early if the result would give a sigma value that's too small. This has value, since those cases would tend to be harder to factor than normal.

Normally I wouldn't suggest it, but since you brought up GPUs I thought I'd mention another option. Perl's ntheory module has a faster sigma than Pari for 64-bit numbers (it gets complicated for bigints). It took my computer 7.6s seconds per 1M block for this at 9 * 10^12, vs. Pari 2.8.0 at 19.2 seconds. That's only a 2x improvement, but maybe worth it if there isn't a good Pari workaround. Also note that on my computer, Pari 2.8.0 ran almost 2x faster than Pari 2.5.5, so make sure you're using a recent version.
[code]
use ntheory qw/:all/;
use Time::HiRes qw(gettimeofday tv_interval);
sub findcyc {
my($a, $e) = @_;
$a++ if $a & 1;
print "running $a through $e\n";
my $start = [gettimeofday];
for (my $b = $a; $b <= $e; $b += 2) {
if ($b % 4 == 0 || $b % 3 == 0) {
my $i = $b;
for my $d (1 .. 4) {
$i = divisor_sum($i) - $i;
if ($i <= $b) {
print "found $b $d\n" if $i == $b;
last;
}
}
}
print " $b/$e ", tv_interval($start), "\n" unless $b % 1000000;
}
}
[/code]

[QUOTE=EdH;385628]I've been running a compiled version of the following routine on several machines in my quest for a 4-cycle sequence. I am running a compiled version, more due to the ease of calling it via a bash script, than because of any speed difference. There is only a slight difference. I have tried many methods of speeding this up further, to include handling the counter outside the main function and different types of loops, etc. As can be seen, I am skipping many numbers in my search, based on information from another member on what they searched. All of this is in a different thread. On my fastest machine, a 1M block (in 9*10^12) is covered in about 44 seconds. The routine does return 2-cycles and found the one known 4-cycle that was within its search pattern.

Finally, on to my questions:

1. Is there room for any appreciable speed improvement in the routine?
2. Is the sigma() function something I could use with a GPU?
3. If 2, then how would I go about (learning and) modifying it to make use of a GPU?
[code]
findcyc(a,e)=
{
if(a%2==1,a+=1);
print("version 20140530b, running ", a, " through ", e);
gettime();
forstep(b=a,e,2,
if(b%4==0 || b%3==0,
i=b;
for(d=1,4,i=sigma(i)-i;
if(i<b,
break(),
if(i==b,
write("../../Desktop/cyclesFound", b, " ", d);
break()
)
)
);
if(b%1000000==0,
print(b,"/",e," ",gettime()/1000+0.," s/M");
write("../../Desktop/cycleStatus", b," ",e)
)
)
)
}
[/code]Thanks for any/all assistance...[/QUOTE]

I've thought of a few things that might help without a gpu, but I can't implement them efficiently if you know what a%12 is you could make a vector of steps to only hit those that are 0 mod 4 or 0 mod 3 with values mod 12 that are 0 mod 2 eliminating the need for the if statement, but you could keep it and try a vector of mod 12 values like [1,0,0,0,1,0,1,0,1,0,0,0][(b%12)+1] which eliminates all the 1 that aren't also 0 mod 2 ( learned this earlier in the thread). I was thinking if you added efficient checks to see if i can be the start of a 3-cycle but that may have too much overhead. Laurv/CRG might know more than me since all my gp are thread engine:single ( mostly because I use the Self-installing distributions for Windows I think) but in newer versions maybe parfor's or something could speed it up by checking multiple candidates at once. or using select/parselect maybe you could select all the ones that fit your conditions quicker.

Thanks for all the help. I will definitely look into all of this as time permits. I did try other methods of choosing the candidates, even trying the selection outside the Pari routine, but none (so far) gained any speed. Many methods slowed the program down.

I had not considered adding a sieve on the front, but I did think about a customized function, but assumed I couldn't do better than what was available. I will have to study that some more.

The Perl info is something I'll have to visit, as well. (I'll have to knock some rust off my Perl programming.) As to the version of Pari, I'll have to see, since every machine has its own install and I compiled the code on each individual computer to be sure any specifics were met, some may not be running the newest.

Of course, the reason for bringing up GPUs is that I have two systems with old NVidia boards (48 cores, each) that are functional with CUDA. All my iterations so far are running as single threads, manually invoked across all available CPU cores per machine. If I could add in some of the NVidia cores, I should gain a little on those machines. But, I have a lot of study to get where I can use them and I'm growing really slow in all my computer "playings" due to other "life interests."

Again, thanks everyone for all the help...

Short Followup
 

Alas, all the several machines I checked (Fedora, Ubuntu, Debian) have Pari 2.5.x installed and nothing newer available via their repositories. Perhaps, rather than opt for individually, manually reinstalling Pari, I should look into the Perl option mentioned. I've also thought of trying one machine with a newer version of Pari and see if the compiled application will run across the others...

Last Followup Here (for now, at least)
 

Well, some preliminary results on one of my slower machines has proven danaj quite correct in the suggestion to take a look at Perl. Some basic tests have shown quite a difference in Perl's favor. My PARI (2.5.5) script takes 116 seconds to cover a 1M block (9*10^12 range), while the Perl (5.18.2) version completes the same range in 33 seconds.

I will leave the PARI thread, for now, in pursuit of more Perl testing.

Thanks danaj!

[QUOTE=EdH;385768]Well, some preliminary results on one of my slower machines has proven danaj quite correct in the suggestion to take a look at Perl. Some basic tests have shown quite a difference in Perl's favor. My PARI (2.5.5) script takes 116 seconds to cover a 1M block (9*10^12 range), while the Perl (5.18.2) version completes the same range in 33 seconds.

I will leave the PARI thread, for now, in pursuit of more Perl testing.

Thanks danaj![/QUOTE]

using my thoughts from thinking back to the asm book I have and CRG's post about my recent code outside of the writes ( which I think we could use listput and write only once at the end) I have broken your code into two codes that get back 0 ms for each individually so 2 instances of pari could in theory do this section of 1 million in 0 ms ( at least for version 2.7.2 but I realized that 2.5.3 was faster for some of my codes that I came up with:



admittedly it also semi-messes up the printing


[CODE]findcyc4(a,e) = {
if(e%4!=0,
e+=4-(e%4)
);
print("version 20140530b, running ",a," through ",e);
gettime();
forstep(b=e,a,-4,
i=b;
d=1;
i=sigma(i)-i;
if(cmp(b,i),
break(),
if(cmp(b,i)==0,
break()
)
);
d+=1;
i=sigma(i)-i;
if(cmp(b,i),
break(),
if(cmp(b,i)==0,
break()
)
);
d+=1;
i=sigma(i)-i;
if(cmp(b,i),
break(),
if(cmp(b,i)==0,
break()
)
);
d+=1;
i=sigma(i)-i;
if(cmp(b,i),
break(),
if(cmp(b,i)==0,
break()
)
);
if(b%1000000==0,
print(b,"/",e," ",gettime()/1000+0.," s/M")
)
)
}[/CODE]

[CODE]findcyc6(a,e) = {
if(e%12!=6,
e+=6-(e%12)
);
print("version 20140530b, running ",a," through ",e);
gettime();
forstep(b=e,a,-12,
i=b;
d=1;
i=sigma(i)-i;
if(cmp(b,i),
break(),
if(cmp(b,i)==0,
break()
)
);
d+=1;
i=sigma(i)-i;
if(cmp(b,i),
break(),
if(cmp(b,i)==0,
break()
)
);
d+=1;
i=sigma(i)-i;
if(cmp(b,i),
break(),
if(cmp(b,i)==0,
break()
)
);
d+=1;
i=sigma(i)-i;
if(cmp(b,i),
break(),
if(cmp(b,i)==0,
break()
)
);
if(b%1000000==0,
print(b,"/",e," ",gettime()/1000+0.," s/M")
)
)
}[/CODE]

maybe someone else can find some more improvements that can take the time down even more.edit: just realized cmp may not be a defined function in my oldest PARI ( which I should uninstall) maybe check that first.

[QUOTE=science_man_88;385842]I have broken your code into two codes that get back 0 ms for each individually so 2 instances of pari could in theory do this section of 1 million in 0 ms[/QUOTE]

:sad:

[QUOTE=axn;385843]:sad:[/QUOTE]

maybe I could just put the second forstep loop in the first code and check e again so maybe I can recombine them.edit: though my ultimate thought if It could ever easily be done (with a slight modification to my repeat script to write to file I think it could be ) was to unroll the forstep loops as we can predict how many times each would be done (e-a)/abs(<respective step>) and yes I forgot the my() for local variables.edit:[B][U] DOH[/U][/B] I realized why I got 0ms I did break() not next(). now I'm back up to just under 18 seconds.

[QUOTE=science_man_88;385845]maybe I could just put the second forstep loop in the first code and check e again so maybe I can recombine them.edit: though my ultimate thought if It could ever easily be done (with a slight modification to my repeat script to write to file I think it could be ) was to unroll the forstep loops as we can predict how many times each would be done (e-a)/abs(<respective step>) and yes I forgot the my() for local variables.edit:[B][U] DOH[/U][/B] I realized why I got 0ms I did break() not next(). now I'm back up to just under 18 seconds.[/QUOTE]
I'm going to play with Perl a bit more as time permits. The basic conversion is running around 10 seconds per million on some of my machines. I am interested in what you come up with, though. I could always swap back later...

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...

[QUOTE=EdH;386269]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...[/QUOTE]

I broke the code down into two codes one handles those that are 0 mod 4 the other 0 mod 6 (Okay technically 6 mod 12 at the time, the combo of 0 mod 2 and 0 mod 3)

but right now the 0 mod 4 code has 3 times the numbers to test until I get the 0 mod 12 numbers over to the 0 mod 6 code and remove it from the need to be checked by the 0 mod 4 code. I also pointed out that to be abundant a number may have a divisor that is either perfect of abundant so I've made a script to test that before going on to test sigma(b)-b ( still ironing things out) only the ones that fail that script need to be tested to show they are abundant before proceeding with d=2 , I'm also thinking of further breaking down the findcyc4 script to avoid an if elseif like construct to find e mod 12 ( because once I move the other check over to the findcyc6 script inputs that are 0 mod 12 don't need to be considered on the findcyc4 script).[CODE]isaorp(x)=my(t=0);forstep(y=#x-1,2,-1,if(!isprimepower(x[y]) && vecsum(divisors(x[y]))>=x[y],return(1),t+=1));if(t==#x-2,return(0))[/CODE]

okay so with some adjustment ( saw a code CRG posted in a sequence in the OEIS that helped as a new start):

[CODE]isaorp(n) = fordiv(n,d,if(5<d && d<sqrtint(n) && !isprimepower(d) && sigma(d)>=2*d,return(1)))[/CODE]

which along with other code alterations brought my attempts at both altering the code from scratch and altering my own versions of the code down from about 3 minutes down to under a minute again. don't know if I'd call isaorp an efficient sigma replacement but I did the math and if you assumed no overlap it can eliminate the need to check after the first sigma in potentially over 20% of cases I still end up having to do sigma at the moment because the way of doing it twice in one go slows down the code because instead of just calling sigma twice in a row it uses it three times the way I originally thought about skipping it. still haven't broken down findcyc4 to the two cases there in theory you could use isaorp before at least the first few sigma because it can can still eliminate potentially 20% of results needing sigma on the larger number itself.

Note that fordiv does a full factorization, then creates the full divisor list, then loops calling the code block. So breaking early doesn't save much time. Measuring, calling fordiv(n,d,break()) takes almost 2x more time than sigma(n). I don't see how this would save time unless you did the checks before calling fordiv.

Your function says the number 100001024 is not abundant/perfect. Isn't it an abundant number?