\\ https://groups.yahoo.com/neo/groups/seqfun/conversations/messages/42 \\ Pari utilities and functions not in the origional distribution \\ by Cino Hilliard Rev. June, 2007 \\ \\ Highlight and copy the script below or copy the entire text for \\ handy reference. Then paste into a blank util.gp file in a \\ working directory. Include the line read "util" to the .gprc file. \\ If you do not want to put this in the .gprc file then just read \\ it in a session \r util. Now when you open a gp session you will \\ have the many useful functions that are listed below. \\ String handling functions \\ asc(chr) Get the ASCII code of character chr \\ Countchr(str,chr) Get the num of occurrences of chr in str. \\ Countmatch(str,m) Get the number of unique occurrences of match \\ string, m in string str. Countmatch2(str,m) Get the number of \\ all occurrences of match string m, in string str. \\ instr(str,m) Get the position of occurrence of match string, \\ m in string str. \\ left(str,n) Get the left n characters from string str. \\ mid(str,s,n) Get a substring of length n from string, str \\ starting at position s in str. \\ rev(str) Get the reverse of the input string str. \\ revvec(v) Get the reverse of the input vector. \\ right(str,n) Get the right n characters from string str. \\ Bit operations \\ bitclear(a,b) Set Result is a with bit b cleared or 0 if it \\ was a 1 \\ bitset(a,b) Set Result is a with bit b set or 1 \\ Number theory functions \\ base(r1,r2,n) Return n base r1 to base r2. bases 2 - 255 are \\ allowed. \\ bdiv(n) The biggest prime divisor of n. \\ composite(n) The n-th composite number \\ diffdigits(n) Return 1 if all digits of n are different else 0 \\ droot(n) Return the digital root of the number n. \\ exponent(n) Return the exponent of n if n is a power > 1 \\ ifactor(n) The vector of integer factors of n with multiplicity. \\ ifactord(n,m=0) The vector of the distinct integer factors of n \\ (without multiplicity). \\ isfourth(n) Return 1 if n is a fourt power. \\ iscube(n) Return 1 if n is a perfect cube else return 0 \\ isprime2(n) Return 1 if f:\sieve\isprime.exe calling \\ primes2-1trill.bin is prime \\ istwin(n) Return -1 if n is a lower twin,1 if upper \\ otherwise return 0 \\ isinteger(n) Return 1 if n is an integer \\ ispal(n) Return 1 if n is a palindromic number \\ issamedr(n) Return 1 if the factors of n with multiplicity \\ have the same digital root else 0. \\ kprimes(n,k) List the indices and primes that differ by k \\ less than or including n \\ Li(x) Logarithmic integral \\ logx(n,b) Logarithm of n to the base b \\ nextsquare(n) The next square number greater than or equal to n. \\ nexttwin(n) Next twin prime pair including or after n. \\ pfgw(n) Return the pfgw output for n \\ pi(n) Prime count function - 2.2.9 has primepi() \\ picics(n) The number of composite-index-composits (cics \\ pronounced kiks) less than or equal to n \\ picomposites(n) The number of composite numbers less than or \\equal to n \\ pip(n) The n-th prime-index-prime (pip) \\ pipips(n) The number of prime-index-primes (pips) less \\ than or equal to n \\ pitwin(n) The number of twin primes less than n. \\ pitwins(n) The number of twin prime pairs <= n. \\ power(n) Return the largest factor for which n is a \\ power \\ primorial(n) The product of primes <= n using the pari primelimit. \\ primorial2(n) The product of primes <= n using isprime \\ primorial3(n) The product of primes <= n using ispseudoprime \\ prime2(n) The nth prime using h:\sieve\prime.exe calling \\ primes2-1trill.bin \\ primeapigcc(n) The nth prime using f:\sieve\primeapigcc.exe \\ calling primes2-1trill.bin \\ primes2(m,n) List the primes from m to n using prime2(n) \\ R(x) Riemann's approx of pi(x) \\ sdiv(n) The smallest prime divisor of n. \\ sumprimes(m,n) Return the sum of the primes(m) to prime(n) \\ sumprimes2(m,n) Return the sum of the primes2(m) to prime2(n) \\ twins(n) List the index and the first n twin primes. \\ twinl(n) The nth lower twin prime \\ twinu(n) The nth upper twin prime \\ \\ Some comments and the actual code begins here \\ These scripts are in util.gp which should be copied into your \\ working directory. Add read "util" to the .gprc file. \\ However the pari ascii set is different from the dos set. Eg., \\ code 227 = symbol for pi in dos. Here we get a. Enclose str,chr, \\ in "" to avoid errors. You may have to use kill to clear the \\ variable if you don't use the "". asc(chr) = \\Get the ascii code of chr { setsearch(ascii,chr) } acodes() = \\ Store the Pari ascii codes globally. \\ User may have to kill ascii in scripts in the session. { local(x); ascii=vector(255); for(x=1,255,ascii[x]=Strchr(x)); } acodes() \\ Define the Pari ascii table globally in standard \\Pari ascii order. base(r1,r2,num) = \\Base conversion 2 to 255. \\Use for compression and encryption { local(RDX,asci,Q,j,ln,p,i,x,d,c); local(pwr,dec,dec2,Qq,k,lnq,r1q,r2q); num=Str(num); if(num == "0",return("0")); d=""; RDX=""; r1q=r1; r2q=r2; dec=0; ln = length(num); for(j=1,ln, c=0; d = mid(num,j,1); until(d==code[c],c++); asci = c-1; lnq = ln-j; dec += asci*r1q^lnq; \\Convert the num in base r1 to decimal); j=1; while(mid(num,j,1)=="0", RDX=concat(RDX,"0"); j++; ); j=0; dec2=dec; while(dec2 > 0, dec2 = floor(dec2/r2); j++; ); forstep(k=j-1,0,-1, pwr = r2q^k; \\ This is just going from base 10 to another base. Q = floor(dec / pwr); \\ convert the decimal number to r2. dec = dec - pwr*Q; RDX = concat(RDX,code[Q+1]); \\ in base r2. Call code[] \\ RDX = concat(RDX,Str(Q)); \\ base. Do it like Maple except \\ RDX = concat(RDX,","); \\seperate by comma leave off brackets ); return(RDX) } bitclear(a,b) = \\ Result is a with bit b in a cleared or 0 if 1. { bitand(a,bitneg(1<1,i=1); x=1; while(x<=lns-lnm+1, if(mid(str,x,lnm)== match,c++;x+=lnm,x++); ); return(c) } countmatch2(str,match) = \\ Count the number of all occurrences of \\ string match in string str. { local(lnm,lns,x,c,i); str=Str(str); \\ This allows leaving quotes off input match=Str(match); c=0; i=0; lns=length(str); lnm=length(match); if(lnm>1,i=1); x=1; while(x<=lns-lnm+1, if(mid(str,x,lnm)==match,c++); x++; ); return(c) } diffdigits(n) = \\ Return 1 if all digits of n are different otherwise return 0 { local(a,i,j,ln,f=1); ln = length(Str(n)); a=eval(Vec(Str(n))); for(i=1,ln, for(j=1,ln, if(i!=j&&a[i]==a[j],f=0); ); ); f; } digits(n) = \\ The vector of the digits of n { return(eval(Vec(Str(n)))) } droot(n) = \\ the digital root of a number. { local(x); x= n%9; if(x>0,return(x),return(9)) } exponent(n) = \\ Return the exponent if n is a power { local(x,ln,j,e=0); ln=omega(n); x=factor(n); e=x[1,2]; for(j=2,ln, if(x[j,2] < e,e=x[j,2]) ); return(e) } find(str,match) = \\ Return the position of the first occurrence of \\string match in string str { local(lnm,lns,tstr,vstr,x,j); vstr=Vec(Str(str)); match=Str(match); lns=length(str); lnm=length(match); for(x=1,lns-lnm+1, tstr=""; for(j=x,x+lnm-1, tstr=concat(tstr,vstr[j]); ); if(match==tstr,return(x)) ); return(0); } getcodes() = \\ Load Ascii character map for base conversion { local(j,ct); \\ vec ->ascii char code = vector(258); \\ code to be called by base(r1,r2,num) ct=0; for(j=48,57, \\ 01-10->0123456789 ct++; code[ct] = Strchr(j); ); for(j=97,122, \\ 11-36->abcdefghijklmnopqrstuvwxyz good to base 36. ct++; \\ Lower case first to be compatible with Lcc code[ct] = Strchr(j); ); for(j=65,90, \\ 37-62->ABCDEFGHIJKLMNOPQRSTUVWXYZ good to base 62. ct++; code[ct] = Strchr(j); ); for(j=32,47, \\ 63-76-> !"#\$%&'()*+,-./ ct++; code[ct] = Strchr(j); ); for(j=58,64, \\ 77-83->:;<=>?@ ct++; code[ct] = Strchr(j); ); for(j=91,96, \\ 84-90->:;<=>?@ ct++; code[ct] = Strchr(j); ); for(j=123,255, \\ 91-224->special ascii char ct++; code[ct] = Strchr(j); ); for(j=1,31, \\ 225-255->special ascii char ct++; code[ct] = Strchr(j); ); } getcodes() \\ This the ascii tree suited for numbers compression ifactor(n) = \\ The vector of the integer factors of n with mult. { local(f,j,k,flist); flist=[]; f=Vec(factor(n)); for(j=1,length(f[1]), for(k = 1,f[2][j],flist = concat(flist,f[1][j]) ); ); return(flist) } ifactord(n) = \\ The vector of the distinct factors of n w/o mult. { local(f,j,k,flist); flist=[]; f=Vec(factor(n)); for(j=1,length(f[1]), flist = concat(flist,f[1][j]) ); return(flist) } instr(str,match) = \\ Get the position of occurrence of match string \\in string str. { local(lns,lnm,x); str=Str(str); \\This allows leaving quotes off input match=Str(match); lnm=length(match); lns=length(str); for(x=1,lns-lnm+1, if(mid(str,x,lnm)==match,return(x)) ) } iscube(n) = { local(r); r = n^(1/3); if(floor(r+.5)^3== n,1,0) } isfourth(n) = \\ Return 1 if n is a fourt power { local(r); r = sqrt(sqrt(n)); if(floor(r+.5)^4== n,1,0) } isinteger(n) = { if(n==floor(n),return(1),return(0)) } issamedr(n) = \\ Return 1 if the factors of n with multiplicity have \\the same digital root else 0. { local(f,a,ln,x); f=0; a=ifactor(n); ln=length(a); for(x=1,ln-1, if(droot(a[x])<>droot(a[x+1]), f=1;break)); if(f==0&ln>1,return(1),return(0)) } ispal(n) = { local(len,s1,s2); n=Str(n); len=length(n)\2; s1 = left(n,len); s2 = right(n,len); if(s1 == rev(s2),return(1),return(0)) } ispal2(n) = { local(flag,len,v,i); n=Str(n); len=length(n); v=Vec(n); flag=1; for(i=1,len\2, if(v[i] <> v[len-i+1],flag=0;break) ); return(flag) } isprime2(n) = \\ 1 if f:\sieve\prime.exe finds n prime { local(x,s); s=concat("f:/sieve/isprime ",Str(n)); s=concat(s," > temp.txt"); \\Must save to a file for correct output system(s); return(read("temp.txt")) } istwin(n) = \\ Return -1 if n is a lower twin,1 if upper,2 if both, otherwise return 0 { local(p1,p2); if(n==5,return(2)); if(isprime(n), p1=n-2; p2=n+2; if(isprime(p1),return(1)); if(isprime(p2),return(-1)); return(0) ) } kprimes(n,k) = \\ List the indices and primes that differ by k less \\than or including n { local(c,x); c=0;forprime(x=3,n,if(isprime(x+k),c++;print(c" "x" "x+k))) } left(str,n) = \\ Get the left n characters from string str. { local(v,tmp,x,ln); v =""; tmp = Vec(str); ln=length(tmp); if(n > ln,n=ln); for(x=1,n, v=concat(v,tmp[x]); ); return(v) } Li(x) = \\ Logarithmic integral { -eint1(log(1/x)) } logx(n,b) = \\ logx(n,b) Logarithm of n to the base b { if(b==0,b=10); \\ Default to base 10 return(log(n)/log(b)) } mid(str,s,n) = \\Get a substring of length n from string str \\starting at position s in str. { local(v,ln,x,tmp); v =""; tmp = Vec(str); ln=length(tmp); for(x=s,s+n-1, v=concat(v,tmp[x]); ); return(v) } nextsquare(n) = \\ The next square number >= n { return(if(issquare(n),n,(floor(sqrt(n))+1)^2)) } nexttwin(n) = \\ Next twin prime pair including or after n. { local(x); forprime(x=n,n+n,if(isprime(x+2),return(x))) } pfgw(n) = { local(s); s=concat("c:/pfgw/pfgw -t -f -e -q",Str(n)); system(s) } pitwins(n) = \\ The number of twin prime pairs <= n. { local(c,x); c=0; forprime(x=3,n, if(isprime(x+2),c++) ); return(c) } pikprimes(n,k) = \\ The number of k difference prime pairs \\less than n. { local(c,x); c=0; forprime(x=3,n, if(isprime(x+k),c++) ); return(c) } pip(n) = \\ The n-th prime-index-prime pip { return(prime(prime(n))) } pipips(n) = \\ The number of prime-index-primes (pips) less than \\or equal to n { local(x,y); x=0; y=0; while(y<=n, x++; y=prime(prime(x)); ); return(x-1); } pipips3(n) = \\ The number of prime-index-prime-index-primes (pips) \\less than or equal to n { local(x,y); x=0; y=0; while(y<=n, x++; y=prime(prime(prime(x))); ); return(x-1); } pipipsR(n) = \\ The number of prime-index-primes (pips) less than \\ or equal to n. Uses the Riemann Integral approximation to pi(n) { local(x,y,yR,ppyR); x=0; y=0; yR = round(R(R(n))); ppyR =prime2(prime2(yR)); if(ppyR > n, while(prime2(prime2(yR)) > n, yR -=1); ); if(ppyR < n, while(prime2(prime2(yR)) < n, yR +=1); if(n <= 10^5,yR--); ); return(yR); } piR(n) = \\ The number of primes less than or equal to n. \\ Uses the Riemann Integral approximation to pi(n) as a boundary. { local(x,y,yR,ppyR); x=0; y=0; yR = round(R(n)); ppyR =prime2(yR); if(ppyR > n, while(prime2(yR) > n, yR -=1); ); if(ppyR < n, while(prime2(yR) > n, yR +=1); if(n <= 10^5,yR--); ); return(yR); } pipipsR3(n) = \\ The number of prime-index-primes (pips) less than \\or equal to n Uses the Riemann Integral \\approximation to pi(n) as a boundary { local(x,y,yR,ppyR); x=0; y=0; yR = R(R(n)); ppyR =prime2(prime2(prime2(yR))); if(ppyR > n, while(prime2(prime2(prime2(yR))) > n, yR -=1); ); if(ppyR < n, while(prime2(prime2(prime2(yR))) < n, yR +=1); if(n <= 10^5,yR--); ); return(yR); } picics(n) = \\ The number of composite-index-composites (cics) pronounced kiks <= n. { local(x,y); x=0; y=0; while(y<=n, x++; y=composite(composite(x)); ); return(x-1); } picomposites(n) = \\ The number of composite numbers less than or \\equal to n { return(n-pi(n)) } pi(n) = \\ pi(n) Prime count function. 2.2.9 has primepi(). \\use pi() to preserve existing code { return(primepi(n)) } pitwin(n) = \\ The number of twin primes less than n. { local(x,c); forprime(x=3,n,if(isprime(x+2),c++)); return(c) } power(n) = \\ Return the largest factor for which n is a power of { if(ispower(n), return(round(n^(1/exponent(n)))) ) } prime2(n) = \\ the nth prime using c:\sieve\prime.exe calling \\ Dell8200\\h8200\\sievedata\\prime2-1trill.bin" { local(x,s); s=concat("c:/sieve/prime ",Str(n)); s=concat(s," > temp.txt"); \\Must save to file for correct output system(s); return(read("temp.txt")) } primorial(n) = \\ The product pf primes <= n using the pari primelimit. { local(p1,x); if(n==0||n==1,return(2)); p1=1; forprime(x=2,n,p1*=x); return(p1) } primorial2(n) = \\ The product of primes <= n using isprime. { local(p1,x); p1=2; forstep(x=3,n,2,if(isprime(x),p1*=x)); return(p1) } primorial3(n) = \\ The product of primes <= n using ispseudoprime. { local(p1,x); p1=2; forstep(x=3,n,2,if(ispseudoprime(x),p1*=x)); return(p1) } primes2(m,n) = \\ List primes from m to n using prime2(n) { local(x,a); a=[]; for(x=m,n, a=concat(a,prime2(x)); ); return(a) } primeGR(n) = \\ A good approximation for the nth prime number using \\the Gram-Riemann approximation of Pi(x) { local(x,px,r1,r2,r,p10,b,e); b=10; p10=log(n)/log(10); if(Rg(b^p10*log(b^(p10+1)))< b^p10,m=p10+1,m=p10); r1 = 0; r2 = 7.18281828; for(x=1,400, r=(r1+r2)/2; px = Rg(b^p10*log(b^(m+r))); if(px <= b^p10,r1=r,r2=r); r=(r1+r2)/2; ); floor(b^p10*log(b^(m+r))+.5); } R(x) = \\ Riemann's approx of Pi(x) { local(j); (sum(j=1,300,moebius(j)*Li(x^(1/j))/j)) } Rg(x) = \\ Gram's Riemann's Approx of Pi(x) { local(n=1,L,s=1,r); L=r=log(x); while(s<10^100*r, s=s+r/zeta(n+1)/n; n=n+1; r=r*L/n); (s) } Rgb(x) = \\ Gram's Riemann's Approx of Pi(x) Broadhurst style { local(n=1,L,s=1,r); L=r=log(x); while(s<10^30*r, s=s+r/zeta(n+1)/n; n=n+1; r=r*L/n); s } rev(str) = \\ Get the reverse of the input string { local(tmp,s,j); tmp = Vec(Str(str)); s=""; forstep(j=length(tmp),1,-1, s=concat(s,tmp[j])); return(s) } revi(str) = \\ Get the integer of the reverse of the input string { local(tmp,s,j); tmp = Vec(Str(str)); s=""; forstep(j=length(tmp),1,-1, s=concat(s,tmp[j])); return(eval(s)) } revvec(v) = \\ Get the reverse of the input vector { local(tmp,j); tmp = v; ln=length(tmp); vector v(ln); forstep(j=ln,1,-1, v[ln-j+1] = tmp[j]); return(v) } right(str,n) = \\ Get the right n characters from string str. { local(v,ln,s,x); v =""; tmp = Vec(str); ln=length(tmp); if(n > ln,n=ln); s = ln-n+1; for(x=s,ln, v=concat(v,tmp[x]); ); return(v) } sumdigits(n) = \\ The sum of the digits of n { local(x,j,s=0); x=digits(n); for(j=1,length(x), s+=x[j]; ); return(s) } sumprimes(m,n) = \\ Return the sum of the mth to the nth prime { local(x); return(sum(x=m,n,prime(x))) } sumprimes2(m,n) = \\ Return the sum of the mth to the nth prime using \\ h:\sievedata\prime2-1trill.bin. { local(x); return(sum(x=m,n,prime2(x))) } sdiv(n) = \\ The smallest prime divisor of n. { local(x); x=ifactor(n); return(x[1]) } twinl(n) = \\ The nth lower twin prime. { local(c,x); c=0; x=1; while(c