mersenneforum.org

mersenneforum.org (https://www.mersenneforum.org/index.php)
-   PARI/GP (https://www.mersenneforum.org/forumdisplay.php?f=155)
-   -   PARI's commands (https://www.mersenneforum.org/showthread.php?t=13636)

science_man_88 2014-10-29 01:23

[QUOTE=danaj;386365]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 [B]much[/B] 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?[/QUOTE]

that's not what it's checking it's checking up to sqrtint(n) for abundant or perfect divisors of the number so the result of sigma of the full number doesn't need to be checked because it has properties that make it abundant. edit: sumdiv doesn't let you break out of it apparently so I used fordiv.

CRGreathouse 2014-10-29 01:50

[QUOTE=science_man_88;386169]finding out if the sum was prime seemed to take less time than the integer division it once held if I remember correct[/QUOTE]

No. Checking for primality in gp involves at least two GCD calculations, each of which take a number of integer divisions.

CRGreathouse 2014-10-29 01:56

[QUOTE=danaj;386365]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.[/QUOTE]

Right. If there was some good reason to do it, you should first factor and then use the factored result:

[code]f=factor(n);
fordiv(f, d,
if (blah(n) && foo(n) && bar(n) && sigma(d) >= 2*n, return(1))
);
sigma(f) ...[/code]

At least you don't have to factor it twice. But it would be better to code your own routine that uses early-abort factorization when the unfactored part is small enough, and the factor bound large enough, to prove that the number must be deficient, or else that the factored part alone suffices to show that the number is abundant. This should be orders of magnitude faster on average.

danaj 2014-10-29 02:53

[QUOTE=science_man_88;386366]that's not what it's checking it's checking up to sqrtint(n) for abundant or perfect divisors of the number so the result of sigma of the full number doesn't need to be checked because it has properties that make it abundant. edit: sumdiv doesn't let you break out of it apparently so I used fordiv.[/QUOTE]

But the number *is* abundant, and the function says it is not. We should continue looking at this number, but isaorp says we don't need to. Either the function is wrong or the name is wrong (because it is returning false when given input that "is abundant or perfect"). Assuming we are asking if n is abundant or perfect, then calling sigma(n) is *faster* than all the checks and early returns inside fordiv, because the setup for fordiv takes more time than all the work inside sigma. [By the way, it turns out that the ntheory library does the same thing for its divisor_sum and fordivisors loop, so no differences there. Not to mention the annoying restriction of no good way to early exit fordivisors. I have solutions for that but they're very inelegant in implementation.]

Based on EdH's program output, if the number is abundant, then we need the result of sigma(n). It seems what we really want is a fast "is this deficient" function, since that is the case where we get to exit early. Since the program limits the search to 4, I suppose a fast "is this perfect" could be used for the last call.

Charles's suggestion of inserting logic inside the factor routine would work since it could either return sigma(n) or truly exit early for deficient numbers (or both abundant and deficient for the last call).

Would it be possible to go from the other direction and see if we could skip the last sigma by constructing which abundant numbers have a sigma equal to 2x one of the perfect numbers? That may be completely infeasible.

As EdH mentioned, perhaps a new thread should be made for this, unless it is of enough interest regarding the Pari commands.

CRGreathouse 2014-10-29 05:29

I'm fine with a new thread.

I think the 'right' function would be one that takes n and k and returns -1, 0, or 1 based on whether sigma_{-1} (n) is less than, equal to, or greater than k, with early exit when possible.

science_man_88 2014-10-29 12:01

[QUOTE=danaj;386372]But the number *is* abundant, and the function says it is not. We should continue looking at this number, but isaorp says we don't need to. Either the function is wrong or the name is wrong (because it is returning false when given input that "is abundant or perfect"). Assuming we are asking if n is abundant or perfect, then calling sigma(n) is *faster* than all the checks and early returns inside fordiv, because the setup for fordiv takes more time than all the work inside sigma. [By the way, it turns out that the ntheory library does the same thing for its divisor_sum and fordivisors loop, so no differences there. Not to mention the annoying restriction of no good way to early exit fordivisors. I have solutions for that but they're very inelegant in implementation.]

Based on EdH's program output, if the number is abundant, then we need the result of sigma(n). It seems what we really want is a fast "is this deficient" function, since that is the case where we get to exit early. Since the program limits the search to 4, I suppose a fast "is this perfect" could be used for the last call.

Charles's suggestion of inserting logic inside the factor routine would work since it could either return sigma(n) or truly exit early for deficient numbers (or both abundant and deficient for the last call).

Would it be possible to go from the other direction and see if we could skip the last sigma by constructing which abundant numbers have a sigma equal to 2x one of the perfect numbers? That may be completely infeasible.

As EdH mentioned, perhaps a new thread should be made for this, unless it is of enough interest regarding the Pari commands.[/QUOTE]

maybe I should of named it isaorpdiv there are abundant numbers that have all divisors but themselves deficient it was a code to calculate them that I started with as they would be one of the lower divisors that could help exit the loop. to go from the other direction you would need to have all intermediate results in theory it's possible using a function like ali or aligen that I made earlier ( admittedly I may have to look them up.) but that doesn't really avoid a sigma and it involves partitions.

danaj 2014-10-31 07:07

I noticed this tonight. Works in 2.6.0, broken in 2.6.1 - 2.8 and tonight's git master. Looks like the commit that introduced it was 3bb8081166c9ac6915977323a4751669ae6c9ed0 in July 2013.

[code](23:05) gp > ispower(-167^10)
*** at top-level: ispower(-167^10)
*** ^----------------
*** ispower: domain error in mplog: argument <= 0
*** Break loop: type 'break' to go back to GP prompt
break> [/code]which ought to return 5.

LaurV 2014-10-31 12:46

I think is not broken, just fixed.
The powering has a higher priority than the unary minus.
It totally makes sense when you write "a^n-b^n", and it is arguable when you write just "-b^n", but to make it right, "(-b)^n" is more clear.

danaj 2014-10-31 16:24

I don't think that is it, as the result is supposed to be negative. -67^10, for example, works just fine. The different versions of Pari also print the same result for -167^10. To remove that effect, change the example to:[code]ispower(-16871927924929095158449)[/code]and one gets the same result (Pari 2.5.5 gives 5, past that commit it starts giving an mplog domain error). This should not produce an error. Note that -101^10 and smaller gives a correct result, but 103 and higher gives the error.

Lest we think this is a problem with evens, we can try -7079^5 which works and -8093^5 which does not. Or -(8093^5) if preferred.

LaurV 2014-10-31 17:06

[strike]Ah, ok, I understand now, but you confused me with that "5". It should return 10, and not 5. Why should it return 5? :surprised:[/strike]
scrap tat, I am stupid. I go to sleep. 0:08 AM here..

danaj 2014-10-31 17:14

Since it is negative, we have to get rid of even powers. So for instance, ispower(-4) returns 0 , ispower(-8) returns 3, ispower(-16) returns 0, ispower(-32) returns 5, etc. For -167^10: -27889^5 yields -167^10 but there isn't any way to raise to an even power and yield a negative result -- hence no 10. Similarly -167^20 yields 5, and -167^(3*5*7*2*2*2*2*2*2) yields 105 (3*5*7).

Edit: nvm, I see you got it. I noticed my module wasn't defined for negative powers so decided to add that since Pari had it and it made sense. It was interesting.

LaurV 2014-10-31 17:16

1 Attachment(s)
I told you I am not arguing with you anymore :razz:

Well, at least they had something to fix from the former versions, and broke it more... I just dug into my archives:
[ATTACH]11908[/ATTACH]

danaj 2014-10-31 17:29

[QUOTE=LaurV;386564]I told you I am not arguing with you anymore :razz:[/QUOTE] :( I'm sorry, I wasn't arguing but trying to write things out in detail to convince myself.

LaurV 2014-11-01 06:07

Haha, why are you so defensive? (don't need to reply, this is a rhetoric question).
That (as the first time too) was a show of respect for you from my side, and not a blame. It was a [URL="http://en.wikipedia.org/wiki/Compliment"]compliment[/URL]. I like to win my arguments, and telling you that I don't want to argue with you means that I evaluated the chances and I think I will lose, and I want to avoid that :razz: Like putting you in the same row with axn, Batalov, and few other guys here that I can count on the fingers of a single hand. It is even harder for me to tell that to a woman, because usually am a misogynist son of a bleach.

science_man_88 2014-12-03 23:15

another attempt at mass TF of mersenne numbers ( in parallel)
 
[CODE]{
my(a=vector(g=1<<12,n,[]),b=[0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1],p=parvector(primepi(1<<15)-1,n,prime(n+1)));
p=setminus(p,parselect(v->v%4==3&&ispseudoprime(v<<1+1)&&b[(v<<1+1)%24+1]&&isprime(v<<1+1)||setminus(MeVec,[v])!=MeVec,p));
p=parvector(2,n,parselect(v->v%3==n,p));
for(x=1,g,
t=x<<1;
m=x%3;
if(m==0,
a[x]=concat(a[x],o=parselect(v->ispseudoprime(l=v*t+1)&&b[(l)%24+1]&&isprime(l)&&!component(Mod(2,l)^v-1,2),concat(p[1],p[2]))),
m==1,
a[x]=concat(a[x],o=parselect(v->ispseudoprime(l=v*t+1)&&b[(l)%24+1]&&isprime(l)&&!component(Mod(2,l)^v-1,2),p[1])),
m==2,a[x]=concat(a[x],o=parselect(v->ispseudoprime(l=v*t+1)&&b[(l)%24+1]&&isprime(l)&&!component(Mod(2,l)^v-1,2),p[2]))
);
if(m!=0,
p[m]=setminus(p[m],o),
p[1]=setminus(p[1],setintersect(p[1],o));
p[2]=setminus(p[2],setintersect(p[2],o))
)
);
a
};[/CODE]

this is still slower than the non paralleled version I made recently:

[CODE]{
my(c=1<<12,a=vector(c,n,[]),b=[0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1]);
forprime(p=3,1<<15,
r=p<<1+1;
j=(p%3==2);
if((p%4==3 && ispseudoprime(r) && b[r%(#b)+1] && isprime(r))||setminus(MeVec[p])!=MeVec,
,
l=1<<((p-1)/2);
forstep(x=if(j,3,2),if(l<c,l,c),[1,2],
q=(x*p)<<1+1;
if(ispseudoprime(q) && b[q%#b+1] && isprime(q) && component(Mod(2,q)^p-1,2)==0,
a[x]=concat(a[x],p);
break()
)
)
)
);
a;
};[/CODE]

is there anything I can do to speed up either version ? they currently run around 7 seconds on my computer.

CRGreathouse 2014-12-03 23:52

I can't figure out what's going on -- I would have expected this to cause an error, with setminus(MeVec[p]) not having enough arguments and all. Why don't you explain what you're trying to do?

science_man_88 2014-12-04 00:11

[QUOTE=CRGreathouse;389054]I can't figure out what's going on -- I would have expected this to cause an error, with setminus(MeVec[p]) not having enough arguments and all. Why don't you explain what you're trying to do?[/QUOTE]

I used a text editor it should read MeVec,p

I'm trying to see if I can successfully use the parallel programming part to actually win out over the non parallel programming by a significant margin, it's doing TF on mersenne numbers 2^q-1 (q in p (a vector of primes) in the parallel version) which get narrowed through conditions known to eliminate exponents early. it then goes on to check the checks are ordered such as to eliminate further checks from being needed. most things that I use more than once get a variable name assigned to them so I don't continue to calculate them for each check. I guess I know a few more things that may decrease the time but one thing I have in the non parallel version seems to slow down the parallel version ( making sure x keeps the number below the sqrt of the mersenne number to be tested) before the loop I divide p into groups mod 6 so I don't test unnecessary x for each v in the parselects. I have tried other things but not in this current code like [][x%3+1] around what is now in the if statement ( laid out as a switch case) but then a lot of things got doubled partly because I didn't have the setminus against o I think. but it's basically eliminating candidates that don't work then testing potentially viable candidates.

CRGreathouse 2014-12-04 00:55

[QUOTE=science_man_88;389055]I'm trying to see if I can successfully use the parallel programming part to actually win out over the non parallel programming by a significant margin, it's doing TF on mersenne numbers 2^q-1[/QUOTE]

OK, so first things first: this should be its own function, which you can call for a single value of q. The function shouldn't know or care about the rest of the vector. Something like this:

[code]merFac(q, lim=1e4)={
forstep(p=2*q+1,lim,2*q,
if (Mod(2, p)^q==1, return(p))
);
0 \\ Couldn't find a factor
};[/code]

Then you can run this function on all the values of your vector:
[code]parapply(merFac, p)[/code]

science_man_88 2014-12-04 01:53

[QUOTE=CRGreathouse;389061]OK, so first things first: this should be its own function, which you can call for a single value of q. The function shouldn't know or care about the rest of the vector. Something like this:

[code]merFac(q, lim=1e4)={
forstep(p=2*q+1,lim,2*q,
if (Mod(2, p)^q==1, return(p))
);
0 \\ Couldn't find a factor
};[/code]

Then you can run this function on all the values of your vector:
[code]parapply(merFac, p)[/code][/QUOTE]

though it doesn't quite do what my code did ( group p based on k) this alone is too much faster to ignore I can do a test up to 2^26 edit:(with the first 2^15 primes greater than 3) just slightly more than what my code could do to 2^12 edit2:( for the primes up to 2^15 greater than 3)

CRGreathouse 2014-12-04 05:51

You may consider, not just for code you post here but for all code you write, (1) naming the functions and (2) using addhelp() to document what they are supposed to do. That will help keep it straight what you're trying to do, which is especially useful as the code gets complex.

And when possible, break the code down into small chunks (as I have) rather than using giant functions. Sometimes complexity is unavoidable, but when you can get around it by chunking things into smaller pieces it's often worthwhile. Another reason this can be helpful: you can see which parts of your program are slow, so you don't waste time optimizing things which are already fast enough.

science_man_88 2015-02-15 14:35

I have a question how do you get PARI/GP to install on android because my tablet is setup to not allow install from unknown sources.

edit never mind apparently it wasn't the setting I thought it was to disable it. this can be deleted.

science_man_88 2015-04-28 21:43

tried to make parallel versions of at least parts of previous functions in theory:

[CODE]
parvecsearch(v,x)={
b=0;
parfor(y=1,#v,
b=b+1;
(v[y]==x),
z,
if(z,
return(1),
if(b==#v,
return(0)
)
)
)[/CODE]
a code attempting part of vecextract in parallel that doesn't work as a function but works without it being function it seems
[CODE]a=[1,2,2,3,4,5];b=[2,4,5];parfor(x=1,#b,b[x]=a[b[x]]);b[/CODE]

the other two versions of this I scripted were with parapply and parvector in place of parfor but i bet this was a waste and already well used/thrown out.

science_man_88 2015-05-14 17:45

[QUOTE=science_man_88;401152]tried to make parallel versions of at least parts of previous functions in theory:

[CODE]
parvecsearch(v,x)={
b=0;
parfor(y=1,#v,
b=b+1;
(v[y]==x),
z,
if(z,
return(1),
if(b==#v,
return(0)
)
)
)[/CODE]
a code attempting part of vecextract in parallel that doesn't work as a function but works without it being function it seems
[CODE]a=[1,2,2,3,4,5];b=[2,4,5];parfor(x=1,#b,b[x]=a[b[x]]);b[/CODE]

the other two versions of this I scripted were with parapply and parvector in place of parfor but i bet this was a waste and already well used/thrown out.[/QUOTE]

I realize now I forgot two simple versions of both, i was coming back to post my attempts at the sieve of Sundaram etc:

for vecsearch I realize it can be done like:

[CODE]a=[2,3,4];#parselect(v->v==3,a)>0[/CODE]

for vecextract it's as simple as:

[CODE]a=[1,3,6,9,12,15,18,21,24,27];b=[5,8,10];c=parapply(v->a[v],b)[/CODE]

and for the sieve of sundaram if you allow specifying the start points it can be done in parallel:

[CODE]a=[1..4004000];parfor(i=1,2000,forstep(x=(sqr(i)+i)<<1,#a,i<<1+1,a[x]=0));parapply(v->v<<1+1,parselect(w->w!=0,vecsort(a,,8)))[/CODE] and yes I know of at least one more optimization by pre-calculating (sqr(i)+i)<<1 outside the loop it's used it.

science_man_88 2015-05-21 14:47

[CODE]parfactor(Y,{lim=Y})=a=[];parforprime(x=1,lim,if(!(Y%x),x,1),y,if(y==1,,D=Y;z=0;until((D%y)!=0,D=D/y;z+=1);a=concat(a,[[y,z]])));Mat(Col(a))[/CODE]

no doubt someone knows how to speed this up (probably with a lookup table but I don't know an efficient way) edit: sorry just thoughts on how to parallel most commands is what comes in my head I guess I should just use parapply.

science_man_88 2015-05-28 16:48

still haven't figured out a vecsort completely but:

parvecmin:

[CODE]parselect(v->!parselect(w->w<v,x),x)[/CODE]

parvecmax:

[CODE]parselect(v->!parselect(w->w>v,x),x)[/CODE]

parsumdiv:

[CODE]a=divisors(y);parsum(x=1,numdiv(y),a[x])[/CODE]

parsetunion:

[CODE]Set(concat(concat(parselect(v->!parselect(w->w==v,y),x),parselect(v->!parselect(w->w==v,x),y)),parselect(v->parselect(w->w==v,y),x)))[/CODE]

parsetminus:

[CODE]Set(parselect(v->!parselect(w->w==v,y),x))[/CODE]

parsetintersect:

[CODE]Set(parselect(v->parselect(w->w==v,y),x))[/CODE]

attempted parcontentvec:

[CODE]x=vecsort(x,,8);a=x[1];b=parvector(#x-1,n,x[n+1]);parapply(p->gcd(a,p),b)[1][/CODE]

parVecrev:

[CODE]parvector(#x,n,x[#x-n+1])[/CODE]

science_man_88 2015-06-25 22:21

[CODE]my(p=15,
a=p-2,
c=1<<p-1,
d=1<<a,
e=[0..d>>1]<<1,
f=parvector(#e,n,
g=parapply(r->[component(Mod(binomial(d,r),c),2),component(Mod(1<<(d-r+1),c),2),component(Mod(3,c)^(r/2),2)],Vec(e[n]));
g=g[1];
prod(h=1,3,
g[h]
)
)
);
vecsum(f)%c[/CODE]

I've had this idea before and even made code once I think but I'd like to make it more parallel some thoughts I had were:

since d-r is used in binomial coefficients, and the other two components that make it into g[1] maybe that and component and mod and some shifts can all be done in parallel. I say shifts even for g because 3^k can be rewritten as 1.5^k*2^k in theory and 1.5 can be rewritten in fraction form but I see that never really ends, also I don't know why I didn't just shift for /2 as I know it will be divisible by two I designed it that way. all this really means is a way to calculate LL residues using properties from :

[URL="https://en.wikipedia.org/wiki/Lucas%E2%80%93Lehmer_primality_test"]Lucas Lehmer test wiki[/URL]

and

[URL="https://en.wikipedia.org/wiki/Binomial_theorem"]binomial theorem wiki[/URL]

the general form of [TEX]S_i [/TEX] and the fact that the powers that make it work would eliminate half the terms and double the rest, is what I'm currently going on.

anyone want to help me speed it up ? It's still crappy speed wise right now and yes I've been using composite exponents to test timings because I wanted to test how far up it could go before I got impatient and it really wasn't far. edit: doh I guess I should break it down to make it quicker.

science_man_88 2015-07-31 15:25

What's new in 2.7.4 ? I have the 64 bit edition now but I don't see much different other than shortcuts to make it do things. and no I haven't made a parvecsort.

science_man_88 2015-07-31 17:11

and now I can't figure out the GPRC again I'm not sure how I messed it up again.

science_man_88 2015-08-07 19:47

alias
 
I know I can alias functions with one operand but I can't get anything other than add to work with more and I tend to get things like : (15:10) gp > [QUOTE]alias("x.pal","palindrome")
(15:11) gp > [2,3,4].pal
*** at top-level: [2,3,4].pal
*** ^---
*** not a function in function call
(15:11) gp > "west".pal
*** at top-level: "west".pal
*** ^---
*** not a function in function call[/QUOTE]

I can't get my RLE code to work though. is there something I'm completely missing ? edit:got my palindrome and RLE scripts working under a different aliases but still not sure about multiple input functions even through I've read about alias("add","_+_").

science_man_88 2015-12-05 21:40

lychrel numbers
 
I've done this before I think but I was watching a numberphile video on 196 and mentioned lychrel numbers so I made a code:

[CODE]revadd(x)= {
a=eval(Vec(Str(x)));
b=Vecrev(a);
until(b==a,
a=a+b;
forstep(x=#a,2,-1,
if(a[x]>9,
a[x-1]+=a[x]\10;a[x]=a[x]%10
)
);
if(a[1]>9,
a=concat(0,a);
a[1]+=a[2]\10;
a[2]=a[2]%10
);
b=Vecrev(a);
print(a","b)
);
}[/CODE]

and I tried to make a parallel version but I can't make it fast ( aka the normal version I just posted was 6 times the speed to get to 16584 digits could anyone help ? :

[CODE]revadd(x)={
a=eval(Vec(Str(x)));
while(a!=Vecrev(a),
print(#a);
a=parvector(#a,i,a[i]+a[#a-i+1]);
until(b==[] || b==[1],
b=parselect(q->q!=#a && a[q+1]>9,
Vec(parselect(r->r<20,a,1))
);
parfor(x=1,#b,
a[b[x]]+1,
c,
a[b[x]]=c;
a[b[x]+1]-=10
);
if(a[1]>9,
a=concat(0,a);
a[1]=a[2]\10;
a[2]-=10
)
)
);
}[/CODE]

CRGreathouse 2015-12-09 15:43

I don't think your binary has parallel support. What do you get with
[code]V=[2^256 + 1, 2^193 - 1];
parvector(#V,i,factor(V[i]))[/code]
(example from "Introduction to parallel GP")? What is
[code]default(nbthreads)[/code]
?

CRGreathouse 2015-12-09 15:48

[QUOTE=science_man_88;406976]What's new in 2.7.4 ? I have the 64 bit edition now but I don't see much different other than shortcuts to make it do things. and no I haven't made a parvecsort.[/QUOTE]

2.7.x are stable versions, there are by design no new features since 2.7.0. You can see the bug fixes and other changes here:
[url]http://pari.math.u-bordeaux.fr/pub/pari/unix/pari-2.7.5.changelog[/url]

science_man_88 2015-12-09 15:48

[QUOTE=CRGreathouse;418692]I don't think your binary has parallel support. What do you get with
[code]V=[2^256 + 1, 2^193 - 1];
parvector(#V,i,factor(V[i]))[/code]
(example from "Introduction to parallel GP")? What is
[code]default(nbthreads)[/code]
?[/QUOTE]

[CODE](11:44) gp > V=[2^256 + 1, 2^193 - 1];
(11:44) gp > parvector(#V,i,factor(V[i]))
%2 = [[1238926361552897, 1; 93461639715357977769163558199606896584051237541638188580280321, 1], [13821503, 1; 61654440233248340616559, 1; 14732265321145317331353282383, 1]]
(11:44) gp > ##
*** last result computed in 3,517 ms.
(11:44) gp > vector(#V,i,factor(V[i]))
%3 = [[1238926361552897, 1; 93461639715357977769163558199606896584051237541638188580280321, 1], [13821503, 1; 61654440233248340616559, 1; 14732265321145317331353282383, 1]]
(11:45) gp > ##
*** last result computed in 3,510 ms.
(11:45) gp > default(nbthreads)
%4 = 1
(11:45) gp > default(nbthreads,24)
(11:45) gp > parvector(#V,i,factor(V[i]))
%6 = [[1238926361552897, 1; 93461639715357977769163558199606896584051237541638188580280321, 1], [13821503, 1; 61654440233248340616559, 1; 14732265321145317331353282383, 1]]
(11:45) gp > ##
*** last result computed in 3,554 ms.
(11:45) gp > vector(#V,i,factor(V[i]))
%7 = [[1238926361552897, 1; 93461639715357977769163558199606896584051237541638188580280321, 1], [13821503, 1; 61654440233248340616559, 1; 14732265321145317331353282383, 1]]
(11:45) gp > ##
*** last result computed in 3,505 ms.
(11:45) gp > default(nbthreads,2)
(11:46) gp > vector(#V,i,factor(V[i]))
%9 = [[1238926361552897, 1; 93461639715357977769163558199606896584051237541638188580280321, 1], [13821503, 1; 61654440233248340616559, 1; 14732265321145317331353282383, 1]]
(11:46) gp > ##
*** last result computed in 3,505 ms.
(11:46) gp > parvector(#V,i,factor(V[i]))
%10 = [[1238926361552897, 1; 93461639715357977769163558199606896584051237541638188580280321, 1], [13821503, 1; 61654440233248340616559, 1; 14732265321145317331353282383, 1]]
(11:46) gp > ##
*** last result computed in 3,500 ms.[/CODE]

my computer has multiple cores so one thread per core could be fastest. I'm not smart enough to figure out the multithreaded engine I bet.

CRGreathouse 2015-12-10 15:43

If you had parallel support in your binary then the parvector call would have failed because V is not declared inline.

Also, I don't think you can meaningfully change default(nbthreads) after starting the program -- you can set the number but I don't think the number of threads used changes. (You could check the manual here.)

science_man_88 2015-12-10 15:51

[QUOTE=CRGreathouse;418808]If you had parallel support in your binary then the parvector call would have failed because V is not declared inline.

Also, I don't think you can meaningfully change default(nbthreads) after starting the program -- you can set the number but I don't think the number of threads used changes. (You could check the manual here.)[/QUOTE]

okay I guess I'll give up trying to make most codes.

science_man_88 2015-12-10 15:57

[QUOTE=CRGreathouse;418808] (You could check the manual here.)[/QUOTE]
[QUOTE="http://pari.math.u-bordeaux.fr/pub/pari/manuals/2.7.0/parallel.pdf"]
However the parallel GP interface does not depend on the multithread interface:[/QUOTE] ?

CRGreathouse 2015-12-10 19:10

Yes, you use the same functions whether you're using MPI or POSIX threads. MPI is best for networks of computers, POSIX threads are best for multicore computers (as a rule of thumb). You don't have support for either version as you can see from the test code I supplied, see p. 7 in Introduction to parallel GP which you linked above.

science_man_88 2015-12-10 19:51

[QUOTE=CRGreathouse;418831]Yes, you use the same functions whether you're using MPI or POSIX threads. MPI is best for networks of computers, POSIX threads are best for multicore computers (as a rule of thumb). You don't have support for either version as you can see from the test code I supplied, see p. 7 in Introduction to parallel GP which you linked above.[/QUOTE]

I don't know how to get anything working except the premade ones especially without a compiler or anything like that. my point with that quote was that you've always seemed premised that the threading interface.

CRGreathouse 2015-12-10 20:11

[QUOTE=science_man_88;418836]I don't know how to get anything working except the premade ones especially without a compiler or anything like that.[/QUOTE]

It's not easy on Windows, you're probably going to have to make do without parallel gp. You can always run multiple instances.

If you were on linux it would be easier.

[QUOTE=science_man_88;418836]my point with that quote was that you've always seemed premised that the threading interface.[/QUOTE]

I have no idea what you mean. :blush:

science_man_88 2015-12-11 23:46

[QUOTE=CRGreathouse;418838]It's not easy on Windows, you're probably going to have to make do without parallel gp. You can always run multiple instances.

If you were on linux it would be easier.



I have no idea what you mean. :blush:[/QUOTE]

okay parallel code aside is there any way to make the non parallel version even faster ( even with modifications to speed up the parallel code it takes about 3 times as long. before that roughly 6 times+ up to around 16000<length(x)<17000 at least in the x=196 sequence).

CRGreathouse 2015-12-12 06:07

OK, so you're asking about speeding up this code:
[code]revadd(x)={
a=eval(Vec(Str(x)));
while(a!=Vecrev(a),
print(#a);
a=parvector(#a,i,a[i]+a[#a-i+1]);
until(b==[] || b==[1],
b=parselect(q->q!=#a && a[q+1]>9,
Vec(parselect(r->r<20,a,1))
);
parfor(x=1,#b,
a[b[x]]+1,
c,
a[b[x]]=c;
a[b[x]+1]-=10
);
if(a[1]>9,
a=concat(0,a);
a[1]=a[2]\10;
a[2]-=10
)
)
);
}[/code]

I'd start with something like this:
[code]revadd(x)=
{
my(a=digits(x),b);
while(a!=Vecrev(a),
print(#a);
a=vector(#a,i,a[i]+a[#a-i+1]);
until(b==[] || b==[1],
b=select(q-> q!=#a && a[q+1]>9, select(r->r<20,a,1));
parfor(x=1,#b,a[b[x]]+1,c,a[b[x]]=c; a[b[x]+1]-=10);
if(a[1]>9,
a=concat(a[1]\10,a);
a[2]-=10
)
)
);
}[/code]

But I don't really understand what it's trying to do. Maybe it would be better to start over. The function takes a number as input, what do you want it to output and/or display?

science_man_88 2015-12-12 12:45

[QUOTE=CRGreathouse;419020]

But I don't really understand what it's trying to do. Maybe it would be better to start over. The function takes a number as input, what do you want it to output and/or display?[/QUOTE]

[YOUTUBE]bN8PE3eljdA[/YOUTUBE]

it's for doing the reverse and add procedure for lychrel numbers.

science_man_88 2015-12-12 13:54

the only other thing that comes up is that [CODE]a=vector(#a,i,a[i]+a[#a-i+1]);[/CODE] is just a way of adding the reverse but I guess I should start to use Vecrev. your's has the parallel version going down under 1 minute to get to 16000+ digits. never mind I can't get yours to copy into the console correctly it gets annoying. and when I add that line back in to the broken up code it seems to stall at 4 digits. never mind even more got lost in the copy but I got it working now. edit: now my thought comes to are there any improvements to the parallel version that I can't work into the not parallel version so that the ratio between the two come closer together still ( your version allowed it to drop from 6 times as slow to do in parallel down to about 2 times as slow.edit2 or 3: I'm running the Vecrev change now to see how it helps I think the parallel version gets caught with overhead and slows down quicker is the problem because the non parallel version doesn't seem as low around 15000 but the parallel starts slowing down around 7000 digits.

science_man_88 2015-12-12 15:26

doh I see something I must of changed back at one point since a[x]\10 can't be more than 1 at last check ( though maybe I did that because it messed up I can't remember) the adding a[x]\10 is adding 1 ( at least I think). if true that would make the non parallel code even faster. edit: so the real hard part is propagating the carries fast enough.

CRGreathouse 2015-12-14 04:54

[QUOTE=science_man_88;419050]it's for doing the reverse and add procedure for lychrel numbers.[/QUOTE]

Here's a program to reverse and add:
[code]A056964(n)=fromdigits(Vecrev(digits(n)))+n[/code]

If you want it to continue until you get a palindrome, try
[code]A033865(n)=my(k); while((k=fromdigits(Vecrev(digits(n)))) != n, n += k); n[/code]

Of course this will loop forever if you get a "Lychrel number".

science_man_88 2015-12-14 14:11

[QUOTE=CRGreathouse;419194]

Of course this will loop forever if you get a "Lychrel number".[/QUOTE]

without fromdigits defined neither is usable.

CRGreathouse 2015-12-14 16:02

[QUOTE=science_man_88;419220]without fromdigits defined neither is usable.[/QUOTE]

You have two choices: either update to a version that has that function, or define your own. The first is a lot faster than the second, if that matters to you.

First option:
[url]http://pari.math.u-bordeaux.fr/pub/pari/windows/snapshots/[/url] for Windows
[url]http://pari.math.u-bordeaux.fr/download.html[/url] otherwise

Second option:
[code]fromdigits(v,b=10)=subst(Pol(v),'x,b)[/code]
or
[code]s(v,b=10)=my(s);for(i=1,#v,s=b*s+v[i]);s[/code]

On my machine, testing with 10^4 digits, the built-in version takes < 1 ms, the first program takes about 4 ms, and the third about 6 ms. With 10^5 digits the timings are 5 ms, 330 ms, and 400 ms.

science_man_88 2015-12-14 16:17

[QUOTE=CRGreathouse;419228]You have two choices: either update to a version that has that function, or define your own. The first is a lot faster than the second, if that matters to you.

First option:
[url]http://pari.math.u-bordeaux.fr/pub/pari/windows/snapshots/[/url] for Windows
[url]http://pari.math.u-bordeaux.fr/download.html[/url] otherwise

Second option:
[code]fromdigits(v,b=10)=subst(Pol(v),'x,b)[/code]
or
[code]s(v,b=10)=my(s);for(i=1,#v,s=b*s+v[i]);s[/code]

On my machine, testing with 10^4 digits, the built-in version takes < 1 ms, the first program takes about 4 ms, and the third about 6 ms. With 10^5 digits the timings are 5 ms, 330 ms, and 400 ms.[/QUOTE]

on my machine without reinstalling an older version fo pari than I now have (2.7.5) I get that your first code gets to 60 digits in 14 seconds and unless I've made a programming error mine gets to over 3000 digits in that time.

CRGreathouse 2015-12-14 16:38

[QUOTE=science_man_88;419231]on my machine without reinstalling an older version fo pari than I now have (2.7.5) I get that your first code gets to 60 digits in 14 seconds and unless I've made a programming error mine gets to over 3000 digits in that time.[/QUOTE]

That should take about 5 microseconds = 0.005 ms, so I imagine you've made a mistake. Would you paste in what you typed that had these timings?

science_man_88 2015-12-14 16:46

[QUOTE=CRGreathouse;419233]That should take about 5 microseconds = 0.005 ms, so I imagine you've made a mistake. Would you paste in what you typed that had these timings?[/QUOTE]

it may be that I put in print(#n) to see where it got to that slowed it down because that's all I changed form the first fromdigits above and your second code that would use it. the only thing I changed from my non parallel code way up in this test was the print(a","b) went to print(#a)

LaurV 2015-12-14 17:06

[QUOTE=science_man_88;419237]the only thing I changed[/QUOTE]
Questions for [URL="https://en.wikipedia.org/wiki/Radio_Yerevan_jokes"]Radio Yerevan[/URL]: "Is is true that the communist party gave Ivan Ivanovich a car?", "Yes it is true, with the difference that it wasn't really a car, but a bicycle, and it wasn't given to him, but taken from him. The rest is true".

CRGreathouse 2015-12-14 18:06

[QUOTE=science_man_88;419237]it may be that I put in print(#n) to see where it got to that slowed it down because that's all I changed form the first fromdigits above and your second code that would use it. the only thing I changed from my non parallel code way up in this test was the print(a","b) went to print(#a)[/QUOTE]

Would you post the actual code? It's hard to imagine that my computer would be 2 million times faster than yours.

science_man_88 2015-12-14 18:15

[QUOTE=CRGreathouse;419244]Would you post the actual code? It's hard to imagine that my computer would be 2 million times faster than yours.[/QUOTE]

I don't have PARI open anymore I might as well not code I suck at it I might delete PARI soon since it doesn't make sense anymore.

CRGreathouse 2015-12-14 18:46

[QUOTE=science_man_88;419245]I don't have PARI open anymore I might as well not code I suck at it I might delete PARI soon since it doesn't make sense anymore.[/QUOTE]

Rather than try to figure out what scenario gave the unusual behavior, why don't you install the latest version of PARI/GP so your new programs can be faster.

science_man_88 2015-12-14 19:16

[QUOTE=CRGreathouse;419248]Rather than try to figure out what scenario gave the unusual behavior, why don't you install the latest version of PARI/GP so your new programs can be faster.[/QUOTE]

I ran the test in the 2.7.5 console. edit:in theory I could just transfer it to the android device I have and run that device into the ground.

CRGreathouse 2015-12-14 19:33

2.7.5 is fine, but it doesn't have the 2.8.0 features like fromdigits.

If you want to continue to use 2.7.5 then you should use one of the functions I provided. You can test their speed here:

[code]fromdigits1(v,b=10)=subst(Pol(v),'x,b);
fromdigits2(v,b=10)=my(s);for(i=1,#v,s=b*s+v[i]);s;
time(f, argument)=gettime(); for(i=1,100,f(argument)); gettime()/100.0;

v=digits(random(10^10^4));

time(fromdigits1, v)
time(fromdigits2, v)
time(fromdigits, v) \\ Only in 2.8.0+[/code]

This runs each version 100 times on the chosen vector v.

For copy/paste convenience, here are the original programs using fromdigits:
[code]A056964(n)=fromdigits(Vecrev(digits(n)))+n
A033865(n)=my(k); while((k=fromdigits(Vecrev(digits(n)))) != n, n += k); n[/code]

Rename the one you like "fromdigits" and this code will work for you.

science_man_88 2015-12-14 19:53

[CODE](15:48) gp > time(revadd,89)
%22 = 0.12000000000000000000000000000000000000
(15:48) gp > fromdigits(n)=my(k); while((k=fromdigits1(Vecrev(digits(n)))) != n, n += k); n
%23 = (n)->my(k);while((k=fromdigits1(Vecrev(digits(n))))!=n,n+=k);n
(15:49) gp > time(fromdigits,89)
%24 = 0.040000000000000000000000000000000000000[/CODE]

results from comparing the first revadd from

[URL="http://mersenneforum.org/showpost.php?p=418400&postcount=2552"]post 2552[/URL]

taking the print out and A033865(n) from your post.

CRGreathouse 2015-12-15 14:35

OK, well good so far -- my program is about 3 times faster. (Not sure why you called it fromdigits though.)

science_man_88 2015-12-15 14:40

[QUOTE=CRGreathouse;419322] (Not sure why you called it fromdigits though.)[/QUOTE]

that's what I thought you meant.

CRGreathouse 2015-12-15 15:54

I made two programs called fromdigits1/fromdigits2 which are replacements for fromdigits, and a program called A033865 which computes A033865. It's best to keep function names meaningful when possible.

science_man_88 2015-12-15 16:57

[QUOTE=CRGreathouse;419328]I made two programs called fromdigits1/fromdigits2 which are replacements for fromdigits, and a program called A033865 which computes A033865. It's best to keep function names meaningful when possible.[/QUOTE]

I found the right snapshot after a few tries and I realize why yours looked so bad I used #n instead of length(Str(n)). you are always too good at this speed things up stuff. and once I got fromdigits in the actual console I had to increase the 100 to 1000 to get a non zero result for timing.

CRGreathouse 2015-12-15 21:15

[QUOTE=science_man_88;419335]I found the right snapshot after a few tries and I realize why yours looked so bad I used #n instead of length(Str(n)). you are always too good at this speed things up stuff.[/QUOTE]

Ah, I see. You were counting words rather than digits for my code, so it looked a lot worse than it was.

Notice that the issue is the lack of Str, not the use of #. You can still write #Str(n) to get the number of digits of n. (Alternatively you could do #digits(n) which has the advantage of letting you change the base if needed.)

[QUOTE=science_man_88;419335]once I got fromdigits in the actual console I had to increase the 100 to 1000 to get a non zero result for timing.[/QUOTE]

Glad to hear it!

I think you'll like the new features in 2.8.0 -- self, fold, sinc, cotanh, maps, binary and hexadecimal constants, variadic functions, etc. It also introduces two large features (central algebras and L-functions) which are pretty technical.

science_man_88 2015-12-15 22:07

[QUOTE=CRGreathouse;419361]I think you'll like the new features in 2.8.0 -- self, fold, sinc, cotanh, maps, binary and hexadecimal constants, variadic functions, etc. It also introduces two large features (central algebras and L-functions) which are pretty technical.[/QUOTE]

maybe once I figure them all out, I'm not sure how to use self completely best I've got so far is a print statement that prints the print call. edit:

[CODE](18:22) gp > loop()=if(x<=y,print(x);x++;eval(self()))
%9 = ()->if(x<=y,print(x);x++;eval(self()))
(18:23) gp > x=2;y=10;loop()
2
3
4
5
6
7
8
9
10[/CODE]

in theory since things can form their own loops in theory any loop structure can now be made this way(except eval doesn't eval with parameters passed currently).

CRGreathouse 2015-12-16 15:04

[QUOTE=science_man_88;419365]maybe once I figure them all out, I'm not sure how to use self completely[/QUOTE]

See
[url]http://rosettacode.org/wiki/Fibonacci_number#Anonymous_recursion[/url]
for an example:
[code]fib(n)={
if(n<2,
n
,
my(s=self());
s(n-2)+s(n-1)
);
}[/code]

This could have been written
[code]fib(n)={
if(n<2,
n
,
my(s=fib);
s(n-2)+s(n-1)
);
}[/code]
but then you have to change the function whenever you rename it, and you can't use it anonymously like you can with self():
[code]n->{
if(n<2,
n
,
my(s=self());
s(n-2)+s(n-1)
)
}[/code]

danaj 2015-12-16 16:43

This reminds me a lot of the __SUB__ feature added to Perl 5.16.0. Since Perl has anonymous subroutines (functions are first class objects), this is especially useful.

Well duh (me), Charles already linked to the RosettaCode page that shows Pari/GP using self() and Perl 5.16+ using __SUB__.

CRGreathouse 2015-12-16 19:20

[QUOTE=danaj;419429]Well duh (me), Charles already linked to the RosettaCode page that shows Pari/GP using self() and Perl 5.16+ using __SUB__.[/QUOTE]

PARI/GP and Perl are alphabetical buddies. :smile:

science_man_88 2015-12-16 19:30

[QUOTE=danaj;419429]This reminds me a lot of the __SUB__ feature added to Perl 5.16.0. Since Perl has anonymous subroutines (functions are first class objects), this is especially useful.

Well duh (me), Charles already linked to the RosettaCode page that shows Pari/GP using self() and Perl 5.16+ using __SUB__.[/QUOTE]

at least you didn't try to reinvent the LL test in another way using fold and make it even slower, or try out [url]http://teknohog.godsong.org/math/dnals/[/url] with powers and parforvec or not even realize what call() is for.

CRGreathouse 2015-12-17 20:16

[QUOTE=science_man_88;419447]or not even realize what call() is for.[/QUOTE]

I haven't found a use for call() yet, even though it's been out for months. I just know what it does so when I do, I'll be able to use it.

But I did happen to use fold() just the other day:
[code]glue(a,b)=a*10^#Str(b)+b
fold(glue, big_vector_of_numbers_to_concatenate)[/code]

science_man_88 2015-12-17 20:41

[QUOTE=CRGreathouse;419543]I haven't found a use for call() yet, even though it's been out for months. I just know what it does so when I do, I'll be able to use it.

But I did happen to use fold() just the other day:
[code]glue(a,b)=a*10^#Str(b)+b
fold(glue, big_vector_of_numbers_to_concatenate)[/code][/QUOTE]

I've used call and fold a few ways but most functions that take more than 2 inputs usually take a vector as one. example:

[CODE]p=11;call(fold,[(a,b)->a^b-b,vector(p-1,i,if(i==1,Mod(4,1<<p-1),2))])[/CODE]

this is a fold version of the LL test called by call.

science_man_88 2015-12-17 22:49

variadic function attempt
 
[CODE]recurVec(a[..])={
my(b=a[1..#a\2-1],
c=a[ceil(#a/2)-1..#a-3],
d=a[#a-1],
f=a[#a]);
for(x=1,f,
e=vector(#a\2-1,
i,
if(d==1,
b[i]*c[i],
d==2,
if(i==1,
b[i]*c[i],
-1*b[i]*c[i]),
d==3,
(-1)^i*b[i]*c[i],
d==4,
(-1)^(i-1)*b[i]*c[i]
)
);
print(vecsum(e))
);
}[/CODE]

This isn't finished yet (ironing out details)! However this is an attempt at a variadic function to deal with difference types of recurrences like in the fibonacci n-step rosetta code ( though I see one there now that I found PARI/gp in the list) and lucas sequences all in one. Just under the first half of the input vector are the coefficients of the recurrence relation. Then all of but the last two remaining in the input vector become the start of the sequence ( they are made to be equal in length at the start) then the last two deal with what type of recurrence it is constantly adding/subtracting terms or alternating ( starting with subtraction or addition ( I know there's functions that do that already (edit: I think) that would be quicker)) edit:and how many new terms to calculate[\edit],and then it will concat the result of summing the new vectors contents as the next term in the recurrence ( yes I see an error in how I setup the vector I think that may be messing up the current output). Any thoughts ? edit: I plan on editing c as I go or changing e as I go so that's why e is declared in the loop. edit: I also made one that may be useful to the rosettacode stuff on ludic numbers.

[CODE]ludic()={
a=[1];
b=[2..100000];
until(#b==0,
a=concat(a,b[1]);
c=ceil(#b/b[1]);
b=setminus(b,vecextract(b,vector(c,i,b[1]*(i-1)+1)))
);
a
}[/CODE]

CRGreathouse 2015-12-18 16:01

[QUOTE=science_man_88;419557]However this is an attempt at a variadic function to deal with difference types of recurrences like in the fibonacci n-step rosetta code ( though I see one there now that I found PARI/gp in the list) and lucas sequences all in one. Just under the first half of the input vector are the coefficients of the recurrence relation. Then all of but the last two remaining in the input vector become the start of the sequence ( they are made to be equal in length at the start) then the last two deal with what type of recurrence it is constantly adding/subtracting terms or alternating ( starting with subtraction or addition ( I know there's functions that do that already (edit: I think) that would be quicker))[/QUOTE]

I think it would be cleared to put the fixed arguments up front and give them their own identifiers rather than keeping them in the variadic portion. Something like
[code]recurVal(d, f, v[..])[/code]

I also think that you should add help for your functions, to explain what the function is supposed to do and how it's called. Something like
[code]addhelp(recurVec, "recurVec(d, f, {x}*): Prints the first f terms of a linear recurrence with [describe more here]. d is a flag, where 1 means ..., 2 means ....");[/code]

In fact I would even suggest writing the addhelp like *before* you write the function so you have a clear idea of what you want before you start coding.

[QUOTE=science_man_88;419557][CODE]recurVec(a[..])={
my(b=a[1..#a\2-1],
c=a[ceil(#a/2)-1..#a-3],
d=a[#a-1],
f=a[#a]);
for(x=1,f,
e=vector(#a\2-1,
i,
if(d==1,
b[i]*c[i],
d==2,
if(i==1,
b[i]*c[i],
-1*b[i]*c[i]),
d==3,
(-1)^i*b[i]*c[i],
d==4,
(-1)^(i-1)*b[i]*c[i]
)
);
print(vecsum(e))
);
}[/CODE][/QUOTE]

I would avoid keeping large vectors in memory if you don't need to. If all you're doing with e is adding up its terms, use sum() so you don't have to spend extra time storing its values in memory and re-reading them.

You can also improve your code by [url=http://www.compileroptimizations.com/category/hoisting.htm]hoisting[/url] the comparison with d. Let's say d = 4 and f = 10^6; your code would check the value of d 4 million times. But if you put the comparison outside, you'd only need 4 comparisons. Like this:

[code]if (d==1,
e=vector(#a\2-1,i,b[i]*c[i])
,
d==2,
e=vector(#a\2-1,i, if(i>1,-1,1)*b[i]*c[i])
,
...[/code]

This should also make it easier for you to do further optimizations. For example, if you decide that all you need with e is to take the sum of its terms, then the case d==1 is just the dot product, which you can find like so:

[code]if (d==1,
esum=b*c~
,
...[/code]

which is faster but also cleaner.

[QUOTE=science_man_88;419557]I plan on editing c as I go or changing e as I go so that's why e is declared in the loop.[/QUOTE]

I can't offer any guidance on that until I know what you want to do.

[QUOTE=science_man_88;419557][CODE]ludic()={
a=[1];
b=[2..100000];
until(#b==0,
a=concat(a,b[1]);
c=ceil(#b/b[1]);
b=setminus(b,vecextract(b,vector(c,i,b[1]*(i-1)+1)))
);
a
}[/CODE][/QUOTE]

You should probably make the 100000 into a parameter so that the user can decide how large to make the returned value.

science_man_88 2015-12-18 17:37

[QUOTE=CRGreathouse;419594]I think it would be cleared to put the fixed arguments up front and give them their own identifiers rather than keeping them in the variadic portion. Something like
[code]recurVal(d, f, v[..])[/code]

I also think that you should add help for your functions, to explain what the function is supposed to do and how it's called. Something like
[code]addhelp(recurVec, "recurVec(d, f, {x}*): Prints the first f terms of a linear recurrence with [describe more here]. d is a flag, where 1 means ..., 2 means ....");[/code]

In fact I would even suggest writing the addhelp like *before* you write the function so you have a clear idea of what you want before you start coding.



I would avoid keeping large vectors in memory if you don't need to. If all you're doing with e is adding up its terms, use sum() so you don't have to spend extra time storing its values in memory and re-reading them.

You can also improve your code by [url=http://www.compileroptimizations.com/category/hoisting.htm]hoisting[/url] the comparison with d. Let's say d = 4 and f = 10^6; your code would check the value of d 4 million times. But if you put the comparison outside, you'd only need 4 comparisons. Like this:

[code]if (d==1,
e=vector(#a\2-1,i,b[i]*c[i])
,
d==2,
e=vector(#a\2-1,i, if(i>1,-1,1)*b[i]*c[i])
,
...[/code]

This should also make it easier for you to do further optimizations. For example, if you decide that all you need with e is to take the sum of its terms, then the case d==1 is just the dot product, which you can find like so:

[code]if (d==1,
esum=b*c~
,
...[/code]

which is faster but also cleaner.
[/QUOTE]

I thought the vector might depend on x at one point so I wasn't sure I'm adding more than linear recurrences I found a way to do the logistic maps (in theory). here's my current code with a few of your optimizations:

took out errored code...

except maybe with a modification of b depending on d.

doh now I see my error I'll go fix that:

[CODE]recurVec(d,f,a[..])={
my(b=a[1..#a/2],
c=a[#a/2+1..#a],
d=d,
f=f
);
if(d==2,
b*=-1;
b[1]*=-1,
d==3,
b*=-1;
b=apply(r->if(r%2==0,r*(-1)),b),
d==4,
b*=-1;
b=apply(r->if(r%2==1,r*(-1)),b),
d==5,
if(#b==1 && c[#c]>0 && c[#c]<1,
b=b[1]*(1-c[#c])
)
);
b=Vecrev(b);
for(x=1,f,
c=concat(c,c[#c-#b+1..#c]*b~)
);
c
}[/CODE]

and a form of addhelp that can have better formatting:

[QUOTE]recurVec(d,f,a[..]),given the coefficient(s) (modified by d,d=1 does nothing, d=2 multiplies all coefficients but the first by -1, d=3 multiplies all even indexed coefficients by -1,d=4 does the
same as d=3 but for odd indexed coefficients,(1 is considered an odd in this process),d=5 if a is appropriate length (2) and has relevant values( 0<a[2]<1 ) attempts the logistic map) and some/a
term (in that order within a) add f more terms to the recurrence, and return it in vector form[/QUOTE]

edit: I've now made d=5 throw an error when the appropriate form isn't given otherwise it seemed to default to something else or does it anyways.

science_man_88 2015-12-18 21:24

I just realized that there's something else I can trim with one condition if I move everything down one for the flag d because d==2,3,4 all have the multiply by -1. and the d==3,4 cases would also benefit from it in that we can then use another condition for apply because they'll become d==2,3 which fit with the mod 2 values r must be in the correct order so we can use r%2==d%2 and do both in one set of commands.

CRGreathouse 2016-01-15 13:22

I started a new thread at
[url]http://www.mersenneforum.org/showthread.php?p=422544[/url]

science_man_88 2016-08-21 12:28

want ot see how many functions you know ?
 
1 Attachment(s)
[CODE]funcIknow()=my(a=apply(r->print("do you know about "r"?");b=input();if(Str(b)=="y",1,0),%1));hammingweight(a)[/CODE]

yes I know that's a bad phrasing of the question as you can know about a function without knowing how to use it or using it correctly. attached is a list I made of 2.8.0 development PARI functions list.

CRGreathouse 2016-08-21 15:03

Maybe also
"algdegree", "alggroup", "alglathnf", "bitprecision", "charconj", "chardiv", "chareval", "charmul", "charorder", "ellsea", "lfuncost", "lfungenus2", "lfunthetacost", "msfromhecke", "msgetlevel", "msgetsign", "msgetweight", "msomseval", "mspadicL", "mspadicinit", "mspadicmoments", "mspadicseries", "mstooms", "rnfidealfactor", "rnfidealprimedec", "serprec", "zncharinduce", "zncharisodd", "znconreychar", "znconreyconductor", "znconreyexp", "znconreylog"

Not sure about
"a", "lfunmfpeters", "lfunsymsq", "lfunsymsqspec", "lgdegree"

CRGreathouse 2016-08-21 15:11

And to make it really fun, use (character-level) trigrams to make fake functions and penalize the player for claiming to know nonexistent functions. :cool:

science_man_88 2016-08-21 16:50

[QUOTE=CRGreathouse;440372]Maybe also
"algdegree", "alggroup", "alglathnf", "bitprecision", "charconj", "chardiv", "chareval", "charmul", "charorder", "ellsea", "lfuncost", "lfungenus2", "lfunthetacost", "msfromhecke", "msgetlevel", "msgetsign", "msgetweight", "msomseval", "mspadicL", "mspadicinit", "mspadicmoments", "mspadicseries", "mstooms", "rnfidealfactor", "rnfidealprimedec", "serprec", "zncharinduce", "zncharisodd", "znconreychar", "znconreyconductor", "znconreyexp", "znconreylog"

Not sure about
"a", "lfunmfpeters", "lfunsymsq", "lfunsymsqspec", "lgdegree"[/QUOTE]

[QUOTE]lfunmfpeters(ldata): ldata corresponding to a normalized eigenform but NOT to the output of lfunsymsq, returns the Petersson square of the form, computed using the symmetric square.

lfunsymsq(ldata,{known=[]}): creates the ldata corresponding to the symmetric square of the modular form ldata, including the search for the conductor and bad Euler factors. known, if present, is
the vector [conductor,[list of Euler factors]], where each Euler factor is of the form [p, a_p] corresponding to the factor 1/(1 - a_pp^(-s)). The result can then be used with the usual lfunxxx
functions. Warning: in the present implementation, only missing Euler factors of degree at most 1 are supported (this is sufficient in most cases, and always if N is squarefree).

lfunsymsqspec(ldata): ldata corresponding either to a normalized eigenform or to the output of lfunsymsq, returns [val, om2], where val is the vector of special values of the symmetric square of
the modular form on the right of the critical strip, and om2 is a period, not canonically normalized. For the moment, only for modular forms of even weight.
[/QUOTE]

not sure how the a and that end part separated on the list but I copied it from PARI and then just used a little find and replace to form that list. all the rest are valid functions.

CRGreathouse 2016-08-21 19:30

So you're running a recent but not bleeding edge version of gp. The master copy on the git repository had lfunmfpeters / lfunsymsq / lfunsymsqspec removed June 30:

[url=http://pari.math.u-bordeaux.fr/cgi-bin/gitweb.cgi?p=pari.git;a=commitdiff;h=8100d1840423249d4a5ceb577632772c629280ce;hp=3a0b53e508ba43b7f0a5f3c632b4970e277f9d8d]Remove lfunmfpeters / lfunsymsq / lfunsymsqspec[/url]
Thu, 30 Jun 2016 14:08:56 +0000 (16:08 +0200)
To be re-introduced
- with proper integration in mf package
- with proper algorithms (minimal twist at the eigenform level)

science_man_88 2016-08-21 19:53

[QUOTE=CRGreathouse;440393]So you're running a recent but not bleeding edge version of gp. The master copy on the git repository had lfunmfpeters / lfunsymsq / lfunsymsqspec removed June 30.[/QUOTE]

I downloaded one recently just I guess for some reason I thought it was auto-installed, installed it now but my antivirus won't let me open it right now.

CRGreathouse 2016-08-23 21:00

[QUOTE=science_man_88;440362][CODE]funcIknow()=my(a=apply(r->print("do you know about "r"?");b=input();if(Str(b)=="y",1,0),%1));hammingweight(a)[/CODE]

yes I know that's a bad phrasing of the question as you can know about a function without knowing how to use it or using it correctly. attached is a list I made of 2.8.0 development PARI functions list.[/QUOTE]

Just for fun I made a version, what do you think?
[url]http://math.crg4.com/gp-function-quiz.html[/url]

I threw in randomly-generated, somewhat plausible function names, though I don't entirely know what to do with them.

science_man_88 2016-08-23 21:29

[QUOTE=CRGreathouse;440534]Just for fun I made a version, what do you think?
[url]http://math.crg4.com/gp-function-quiz.html[/url]

I threw in randomly-generated, somewhat plausible function names, though I don't entirely know what to do with them.[/QUOTE]

my results:

You knew 2 out of 10 functions. A 95% confidence interval around these scores suggests that you know 4.6% to 52.0% of the GP functions.
Medium test You knew 3 out of 25 functions. A 95% confidence interval around these scores suggests that you know 3.3% to 30.8% of the GP functions.
Long test You knew 28 out of 100 functions. A 95% confidence interval around these scores suggests that you know 20.1% to 37.5% of the GP functions.
Full test (long) You knew 172 out of 859 functions. A 95% confidence interval around these scores suggests that you know 17.5% to 22.8% of the GP functions.

my criteria for knowing them wasn't always rigorous I mostly picked the one's that I've played around with most ( though some are one's I've read more up about)

CRGreathouse 2016-08-23 21:58

Very good -- you didn't trigger a single false positive! :smile:

I meant to take out the confidence interval on the full test -- obviously once you've exhausted all the possibilities you don't want/need an interval.

CRGreathouse 2016-08-23 22:00

Having somewhat less patience I got
[INDENT]You knew 18 out of 25 functions. A 95% confidence interval around these scores suggests that you know 52.2% to 85.9% of the GP functions.[/INDENT]
on the medium test. Incidentally, that seems way too high, but that's the luck of the draw.

science_man_88 2016-08-23 22:13

have you found a use for call yet like you said you might once it got a definition/example ? my weirdest thought so far is comparing multiple ways of doing things maybe call parapply or another grouping of functions to compare how fast they are.

CRGreathouse 2016-08-24 05:26

[QUOTE=science_man_88;440538]have you found a use for call yet like you said you might once it got a definition/example ? my weirdest thought so far is comparing multiple ways of doing things maybe call parapply or another grouping of functions to compare how fast they are.[/QUOTE]

In the process of writing the tool above I used the new ES6 [url=https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Operators/Spread_operator]spread[/url] operator, which can be used to simulate call. GP:
[code]call(f, v)[/code]
is essentially the same as JavaScript (ES6)
[code]f(...v)[/code]

It's a beautiful syntax, especially coming from a language not known for beauty. (It's a nice idea cluttered with lots of junk and some partial attempts to clean it up.)

CRGreathouse 2016-08-24 06:24

I updated the page, adding corrections for finite size on the confidence intervals (most important when doing the long test) and dropped the CI from the full test (because it's meaningless). I also added percentages and a link to the PARI/GP website.
[url]http://math.crg4.com/gp-function-quiz.html[/url]

science_man_88 2017-07-29 12:55

2.9.3 is here
 
if I did the copy and paste correctly it has 866 functions. what's changed and what hasn't :

I also started a[URL="https://math.stackexchange.com/questions/2375710/pari-gp-math-commands-and-how-many-duplications-are-there-and-what-are-there"] stack exchange question[/URL] but I fear it will be closed soon.

CRGreathouse 2017-07-30 03:07

If you'd like to update my page, all you need is to replace the list of functions on line 162 with your own list of functions. It has 859 functions, but there are a lot more than 7 new functions -- you'll just need to move to 2.10 to get them. :smile:

science_man_88 2017-07-30 10:31

[QUOTE=CRGreathouse;464500]If you'd like to update my page, all you need is to replace the list of functions on line 162 with your own list of functions. It has 859 functions, but there are a lot more than 7 new functions -- you'll just need to move to 2.10 to get them. :smile:[/QUOTE]

not sure how to find 2.10 also I suck at compiling etc. if that's what it takes. FWIW I also don't see the list you are talking about. if you mean the code tags in the inspect tab no wonder you are rated as unsecure

CRGreathouse 2017-08-01 04:56

[QUOTE=science_man_88;464518]not sure how to find 2.10 also I suck at compiling etc. if that's what it takes.[/QUOTE]

If you're on Windows, you can get them here:
[url]http://pari.math.u-bordeaux.fr/pub/pari/windows/snapshots/[/url]

I recommend installing the latest stable 64-bit version, then extracting the latest gp64-gmp binary into the installed directory. That gives a best of both worlds (for the most part): the newest features on top of the extra features not found in the standalone binary.

If you're on linux, I have install instructions here:
[url]http://math.crg4.com/software.html#pari[/url]
though there are simpler ways as well, for example using a package manager.

If you're on Mac OS, god help you. Braver souls than I have tried and failed.

[QUOTE=science_man_88;464518]FWIW I also don't see the list you are talking about. if you mean the code tags in the inspect tab no wonder you are rated as unsecure[/QUOTE]

It's in the source of the page I linked to. Most browsers will display this on Ctrl+U. I don't understand your last comment; a web page needs to pass HTML to the user to be displayed, having it accessible is a necessary part of making the Internet work.

science_man_88 2017-08-01 11:42

[QUOTE=CRGreathouse;464635]If you're on Windows, you can get them here:
[url]http://pari.math.u-bordeaux.fr/pub/pari/windows/snapshots/[/url]

I recommend installing the latest stable 64-bit version, then extracting the latest gp64-gmp binary into the installed directory. That gives a best of both worlds (for the most part): the newest features on top of the extra features not found in the standalone binary.

If you're on linux, I have install instructions here:
[url]http://math.crg4.com/software.html#pari[/url]
though there are simpler ways as well, for example using a package manager.

If you're on Mac OS, god help you. Braver souls than I have tried and failed.



It's in the source of the page I linked to. Most browsers will display this on Ctrl+U. I don't understand your last comment; a web page needs to pass HTML to the user to be displayed, having it accessible is a necessary part of making the Internet work.[/QUOTE]
okay so 2.10.0 has 973 functions , I guess my point is that I wasn't sure how I would get it onto the site I'm not a developer of the site and shouldn't have easy access to code. It's like on pinterest, they pop something up to stop you seeing more images if you aren't logged in but you can just delete it if you want because it's in the coding directly.

CRGreathouse 2017-08-01 15:16

[QUOTE=science_man_88;464654]okay so 2.10.0 has 973 functions[/QUOTE]

Wow, that jumped a lot more than I would have guessed! :smile:

[QUOTE=science_man_88;464654]I guess my point is that I wasn't sure how I would get it onto the site I'm not a developer of the site and shouldn't have easy access to code.[/QUOTE]

Ah, I understand where you're coming from now. You couldn't change it on my site (for the security reasons you were getting at), but you could change the copy on your computer.

danaj 2017-08-01 16:07

[QUOTE=CRGreathouse;464635]If you're on Mac OS, god help you. Braver souls than I have tried and failed.[/quote]Using 10.11.6, building from git seems to work fine, though there were a few hurdles with readline support -- had to get homebrew and install readline.

CRGreathouse 2017-08-01 17:34

[QUOTE=danaj;464690]Using 10.11.6, building from git seems to work fine, though there were a few hurdles with readline support -- had to get homebrew and install readline.[/QUOTE]

I heard that it got better in more recent versions, but I didn't know if it was true. That's exciting!

If you're interested, give me the steps (here or in a PM) and I'll include them in my build instructions. I need to update my Windows instructions anyway. Is it as simple as "brew install readline" and then "./Configure && make all && sudo make install"? I don't have a Mac handy to test.

danaj 2017-08-01 19:08

[QUOTE=CRGreathouse;464698]I heard that it got better in more recent versions, but I didn't know if it was true. That's exciting!

If you're interested, give me the steps (here or in a PM) and I'll include them in my build instructions. I need to update my Windows instructions anyway. Is it as simple as "brew install readline" and then "./Configure && make all && sudo make install"? I don't have a Mac handy to test.[/QUOTE]

I'd need a fresh system to really tell for instruction purposes. I installed it last October and I don't recall it being hard other than readline support which I just ignored until February (because getting it to work is indeed an exercise in frustration until you discover the key).

You definitely have to do "./Configure --with-readline=/usr/local" for Pari once "brew install readline" is done because without that it will find Apple's libedit (crippled readline) and make a non-readline version. For "make all" you'd also need TeX, but "make gp" works fine for me.

science_man_88 2017-08-01 20:16

1 Attachment(s)
[QUOTE=CRGreathouse;464682]Wow, that jumped a lot more than I would have guessed! :smile:



Ah, I understand where you're coming from now. You couldn't change it on my site (for the security reasons you were getting at), but you could change the copy on your computer.[/QUOTE]

I also don't care to edit most peoples stuff. I make it less useful. here's the list:


All times are UTC. The time now is 23:20.

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.