1
0
Fork 0
Libraries for bc and dc.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

731 lines
21 KiB

#!/usr/local/bin/bc -l
### Primes.BC - Primes and factorisation (rudimentary)
## All factor finding is done by trial division meaning that many
## functions will eat CPU for long periods when encountering
## certain numbers. Primality testing uses better techniques and
## is much faster if no factors are required.
## e.g. 2^503-1 is identified as non-prime by the primality testers
## but no factors will be found in any sensible amount of time
## through trial division.
##
## Steps have been taken to make the trial division as fast as
## possible, meaning much code re-use.
max_array_ = 4^8-1
# Greatest common divisor of x and y - stolen from funcs.bc
define int_gcd(x,y) {
auto r,os;
os=scale;scale=0
x/=1;y/=1
while(y>0){r=x%y;x=y;y=r}
scale=os
return(x)
}
### Primality testing ###
# workhorse function for int_modpow and others
define int_modpow_(x,y,m) {
auto r, y2;
if(y==0)return(1)
if(y==1)return(x%m)
y2=y/2
r=int_modpow_(x,y2,m); if(r>=m)r%=m
r*=r ; if(r>=m)r%=m
if(y%2){r*=x ; if(r>=m)r%=m}
return( r )
}
# Raise x to the y-th power, modulo m
define int_modpow(x,y,m) {
auto os;
os=scale;scale=0
x/=1;y/=1;m/=1
if(x< 0){print "int_modpow error: base is negative\n"; x=-x}
if(y< 0){print "int_modpow error: exponent is negative\n";y=-y}
if(m< 0){print "int_modpow error: modulus is negative\n"; m=-m}
if(m==0){print "int_modpow error: modulus is zero\n"; return 0}
x=int_modpow_(x,y,m)
scale=os
return( x )
}
## Pseudoprime tests
# Global variable to limit the number of Rabin-Miller iterations to try
# 0 => Run until sure number is prime
rabin_miller_maxtests_=0
# Uses the Rabin-Miller test for primality
# uses a shortcut for numbers < 300 decimal digits
define is_rabin_miller_pseudoprime(p) {
auto os,a,inc,top,next_a,q,r,s,d,x,c4;
os=scale;scale=0
if(p!=p/1){scale=os;return 0}
if(p<=(q=F+2)){scale=os;return(p==2||p==3||p==5||p==7||p==B||p==D||p==q)}
s=0;d=q=p-1;x=d/2;while(0==d-x-x){.=s++;d=x;x=d/2}
if(p<A^(1+3*A*A)){
# Takes a few liberties for the sake of speed compared to a
# . full RM; Assumes small primes and perrin tests have been run
inc=(p-4)/(C+length(p));if(inc<1)inc=1
top=q
} else {
# This bizarre construct sets inc to an approximation of
# . log to base 4 of p, meaning the loop below should
# . run enough times to guarantee that p is prime
inc=((B*B+F+F)*B*(length(p)+1))/A^3;if(inc<1)inc=1
if(!inc%2).=inc++ # inc needs to be odd to ensure good mix of candidates
top=inc*(inc+1)
}
if(rabin_miller_maxtests_){
rabin_miller_maxtests_/=1
if(rabin_miller_maxtests_<=0){
rabin_miller_maxtests_=0
print "Warning: rabin_miller_maxtests_ set to invalid value. Now = 0\n"
print " This calculation may take longer to run than expected\n"
}else{
inc=top/rabin_miller_maxtests_+1
}
}
for(a=2;a<top;a+=inc){
next_a=0
x=int_modpow_(a,d,p)
if(x!=1&&x!=q){
for(r=1;r<s;r++){
x*=x;x%=p
if(x==1){scale=os;return 0}#composite
if(x==q){next_a=1;break}
}#end for
if(!next_a){scale=os;return 0}
}#end if
}#end for
scale=os;return 1
}
# Determine whether p is a Perrin pseudoprime
# returns 0 if definitely composite
# returns 1 if possibly prime (but not definitely)
define is_perrin_pseudoprime(p) {
auto os,i,h,rp,m[],r[],t[];#,m[];
os=scale;scale=0
if(p!=p/1){scale=os;return 0}
if(p==2){scale=os;return 1}
if(p<2||p%2==0){scale=os;return 0}
#set rp to reverse of p
# would love to use int_modpow to calculate powers but it doesn't work
# on matrices!
# could use an array to store bits of p, but arrays are limited to
# 65536 elements; integers have no such limits
rp=0;i=p;while(i){h=i/2;rp+=rp+i-h-h;i=h}
# Quick Mersenne test; if rp == p then p could be 2^i - 1 for some i.
# if i is not prime then neither is p, and it's faster to check i first.
# this test is unnecessary, but since half the work is already done
# we might as well check
if(rp==p){
for(i=0;rp;i++)rp/=2
if(2^i-1==p&&!is_perrin_pseudoprime(i)){scale=os;return 0}
rp=p # restore rp; if we're here then the test failed
}
# m[]={0,1,1;1,0,0;0,1,0}
#m[0]=0; m[1]=1; m[2]=1;
#m[3]=1; m[4]=0; m[5]=0; # Never actually used
#m[6]=0; m[7]=1; m[8]=0;
# r[]={1,0,0;0,1,0;0,0,1}=Unit
r[0]=1; r[1]=0; r[2]=0;
r[3]=0; r[4]=1; r[5]=0;
r[6]=0; r[7]=0; r[8]=1;
while(rp) {
# square r
t[0]=r[1]*r[3]; t[1]=r[2]*r[6]; t[2]=r[0]+r[4]
t[3]=r[2]*r[7]; t[4]=r[0]+r[8]
t[5]=r[1]*r[5]; t[6]=r[5]*r[6]; t[7]=r[5]*r[7]; t[8]=r[4]+r[8]
t[9]=r[2]*r[3]; t[A]=r[3]*r[7]; t[B]=r[1]*r[6]
r[0]*=r[0];r[0]+=t[0]+t[1]; r[1]*=t[2];r[1]+=t[3]; r[2]*=t[4];r[2]+=t[5]
r[3]*=t[2];r[3]+=t[6]; r[4]*=r[4];r[4]+=t[0]+t[7]; r[5]*=t[8];r[5]+=t[9]
r[6]*=t[4];r[6]+=t[A]; r[7]*=t[8];r[7]+=t[B]; r[8]*=r[8];r[8]+=t[1]+t[7]
h=rp/2
if(rp-h-h){# odd
# multiply r by m
# this is a hack that assumes m is {0,1,1;1,0,0;0,1,0}
# without actually using m.
# if m is changed, this will need to be rewritten
t[0]=r[0]+r[2]; t[1]=r[3]+r[5]; t[2]=r[6]+r[8]
r[2]=r[0]; r[0]=r[1]; r[1]=t[0]
r[5]=r[3]; r[3]=r[4]; r[4]=t[1]
r[8]=r[6]; r[6]=r[7]; r[7]=t[2]
}
for(i=0;i<9;i++)r[i]%=p
rp=h
}
r[0]=(2*r[6]+3*r[8])%p
scale=os;
return (r[0]==0)
}
# Determine whether x is possibly prime through division by small numbers
define is_small_division_pseudoprime(x) {
auto os,j[],ji,sx,p,n;#oldscale,jump,jump-index,sqrtx,prime,nth
if(x<2)return 0;
os=scale;scale=0
if(x!=x/1){scale=os;return 0}
j[0]=4;j[1]=2;j[2]=4;j[3]=2;j[4]=4;j[5]=6;j[6]=2;j[7]=6
for(p=2;p<7;p+=p-1){if(x==p){scale=os;return 1};if(x%p==0){scale=os;return 0}}
sx=sqrt(x);if(sx>(n=A^5))sx=n;# 100000(decimal) upper limit
if(prime[A^4]){ #primes-db is present
for(n=4;(p=prime[n])<=sx;n++)if(x%p==0){scale=os;return 0}
} else {
ji=7;p=7;n=2*A-1
while(p<=sx){
if(x%p==0){scale=os;return 0}
if(ji++==8)ji=0;p+=j[ji];
if(p>n)while(p%7==0||p%B==0||p%D==0||p%n==0){if(ji++==8)ji=0;p+=j[ji]}
}
}
scale=os;return(1)
}
## Primality / Power-freedom tests
define is_prime(x) {
# It is estimated that all numbers will not be misidentified
# using the tests below, but it may take time
if(!is_small_division_pseudoprime(x))return 0
if(x<A^A)return 1 # pairs with the A^5 in is_s.d.pp()
if(x<A^(1+3*A*A)||rabin_miller_maxtests_){
# 300 digits or less; faster doing RM, then PP
# and shortcut RM may miss something (hard to prove)
if(!is_rabin_miller_pseudoprime(x))return 0
if(!is_perrin_pseudoprime(x))return 0
} else {
# Tests reversed because full RM is slower than PP
if(!is_perrin_pseudoprime(x))return 0
if(!is_rabin_miller_pseudoprime(x))return 0
}
return 1
}
# Determine whether x is prime through trial division only
define is_prime_td(x) {
auto os,j[],ji,sx,p,n;#oldscale,jump,jump-index,sqrtx,prime,nth
if(x<2)return 0;
os=scale;scale=0
if(x!=x/1){scale=os;return 0}
j[0]=4;j[1]=2;j[2]=4;j[3]=2;j[4]=4;j[5]=6;j[6]=2;j[7]=6
for(p=2;p<7;p+=p-1){if(x==p){scale=os;return 1};if(x%p==0){scale=os;return 0}}
if(!is_perrin_pseudoprime(x)){scale=os;return 0}#cheat a bit
sx=sqrt(x);p=7;ji=7
if(prime[max_array_])for(n=4;n<=max_array_&&(p=prime[n])<=sx;n++)if(x%p==0){scale=os;return 0}
if(p>7)ji=4#assume p is now prime[max_array_]
n=2*A-1
while(p<=sx){
if(x%p==0){scale=os;return 0}
if(ji++==8)ji=0;p+=j[ji];
if(p>n)while(p%7==0||p%B==0||p%D==0||p%n==0){if(ji++==8)ji=0;p+=j[ji]}
}
scale=os;return(1)
}
### Storage and output of prime factorisations ###
# Output the given array interpreted as prime factors and powers thereof
# . this function plus fac_store() make for a "delayed" equivalent
# . to the fac_print() function
define printfactorpow(fp[]) {
auto i,c;
for(i=0;fp[i];i+=2){
print fp[i]
if((c=fp[i+1])>1)print "^",c
if(fp[i+2])print " * "
}
print"\n"
return (fp[1]==1&&fp[2]==0) # fp[] is prime?
}
# Workhorse function for the below
# . for retaining a copy of the last calculated factorisation
# . in factorpow global array to save time if further functions
# . are to be called on same number
define factorpow_set_(fp[]) {
auto i;
for(i=0;fp[i];i++)factorpow[i]=fp[i]
return factorpow[i]=factorpow[i+1]=0;
}
# Workhorse function for the below
# . appends newly found factor and power thereof to the provided array
# . outputs that information if the print_ flag is set
define fac_store_(*fp[],m,p,c,print_) {
auto z;
if(!m%2).=m++ # m should be position of last element and thus odd
# even elements are prime factor, odd elements are how many.
# 9 -> {3,2} -> 3^2 , 60 -> {2,2,3,1,5,1} -> 2^2*3^1*5^1
# negative c means we know this is the end and we can write two zeroes
z=0;if(c<0){z=1;c=-c}
fp[++m]=p;fp[++m]=c
if(print_){
print p;if(c>1)print "^",c
if(!z){print " * "}else{print "\n"};
}
if(z){fp[++m]=0;fp[++m]=0}
return m
}
# Workhorse function for the below
# . performs action that otherwise occurs three times
# . relies on inherited scope (POSIX style)
# . returns 0 if parent should also return
define fac_sp_innerloop_() {
for(c=0;x%p==0;c++)x/=p
if(c){
if(x==1)c=-c
m=fac_store_(fp[],m,p,c,print_);
if(x==1)return factorpow_set_(fp[]); # = 0
if(is_prime(x)){
m=fac_store_(fp[],m,x,-1,print_);
return factorpow_set_(fp[]); # = 0
}
}
return 1;
}
# Workhorse function for the below
# . factorises x through trial division, using the above functions
# . for output, storage, retention, etc.
define fac_sp_(*fp[],x,print_) {
auto os,j[],ji,sx,p,c,n,m,f;#oldscale,jump,jump-index,sqrtx,prime,count,nth,mth
os=scale;scale=0;x/=1
# Check to see if last calculation was the same as this one - save work
f=1;for(m=0;p=factorpow[m]&&f<=x;m+=2)f*=(fp[m]=p)^(fp[m+1]=factorpow[m+1]);
if(f==x){
if(print_).=printfactorpow(fp[]);
scale=os;return fp[m]=fp[m+1]=0;
}
# Main algorithm
m=-1
if(x<0){m=fac_store_(fp[],m,-1,1,print_);x=-x}
if(x<=1||is_prime(x)){m=fac_store_(fp[],m,x,-1,print_);scale=os;return (x>1)}
j[0]=4;j[1]=2;j[2]=4;j[3]=2;j[4]=4;j[5]=6;j[6]=2;j[7]=6
for(p=2;p<7;p+=p-1)if(!fac_sp_innerloop_()){scale=os;return 0} #1
sx=sqrt(x);p=7;ji=7
if(prime[max_array_])for(n=4;n<=max_array_&&(p=prime[n])<=sx;n++){
if(!fac_sp_innerloop_()){scale=os;return 0} #2
}
if(p>7)ji=4#assume p is now prime[max_array_]
n=2*A-1;sx=sqrt(x)
while(p<=sx){
if(!fac_sp_innerloop_()){scale=os;return 0} #3
if(c)sx=sqrt(x)
if(ji++==8)ji=0;p+=j[ji];
if(p>n)while(p%7==0||p%B==0||p%D==0||p%n==0){if(ji++==8)ji=0;p+=j[ji]}
}
if(x>1)m=fac_store_(fp[],m,x,-1,print_)
scale=os;return factorpow_set_(fp[]);
}
# Output the prime factors of x with powers thereof
# . displays factors as they are found which allows more chance
# . of some factors being output before becoming bogged down
# . finding larger factors by trial division
define fac_print( x) { auto a[];x=fac_sp_( a[],x,1);return ( a[1]==1&& a[2]==0); }
# Store the prime factors of x into the given array
define fac_store(*fp[],x) { x=fac_sp_(fp[],x,0);return (fp[1]==1&&fp[2]==0); }
# Can be slow in the case of a large spf
define smallest_prime_factor(x) {
auto os,j[],ji,sx,p,n;#oldscale,jump,jump-index,sqrtx,prime,nth
os=scale;scale=0;x/=1
if(x<0){scale=os;return -1}
if(x<=1||is_prime(x)){scale=os;return x}
j[0]=4;j[1]=2;j[2]=4;j[3]=2;j[4]=4;j[5]=6;j[6]=2;j[7]=6
for(p=2;p<7;p+=p-1)if(x%p==0){scale=os;return p}
sx=sqrt(x);p=7;ji=7
if(prime[max_array_])for(n=4;n<=max_array_&&(p=prime[n])<=sx;n++)if(x%p==0){scale=os;return p}
if(p>7)ji=4#assume p is now prime[max_array_]
n=2*A-1;sx=sqrt(x)
while(p<=sx){
if(x%p==0){scale=os;return p}
if(ji++==8)ji=0;p+=j[ji];
if(p>n)while(p%7==0||p%B==0||p%D==0||p%n==0){if(ji++==8)ji=0;p+=j[ji]}
}
scale=os;return(-sx) #if we ever get here, something is wrong!
}
# Non trivial = slow
define largest_prime_factor(x) {
auto i,fp[];
.=fac_store(fp[],x);
for(i=0;fp[i];i+=2){}
return fp[i-2];
}
# Largest prime power within x
# e.g. 992 = 2^5*31 and 2^5>31 so 992 -> 2^5 = 32
define largest_prime_power(x) {
auto i,fp[],p,max;
.=fac_store(fp[],x);
for(i=0;(p=fp[i]^fp[i+1])!=1;i+=2)if(max<p)max=p
return max;
}
# Return powerfree kernel of x (largest powerfree number which divides x)
define rad(x) {
auto i,r,f,fp[];
.=fac_store(fp[],x)
r=1;for(i=0;f=fp[i];i+=2)r*=f
return r;
}
# Return square part of x
define squarepart(x) {
auto os,i,r,f,p,fp[];
.=fac_store(fp[],x)
os=scale;scale=0
r=1;for(i=0;f=fp[i];i+=2){p=fp[i+1];p-=p%2;r*=f^p}
scale=os;return r
}
# Return squarefree core of x
define core(x) {
auto os;
os=scale;scale=0
x/=squarepart(x)
scale=os;return x
}
# Count the number of (non-unique) prime factors of x
# e.g. 60 = 2^2*3^1*5^1 -> 2 + 1 + 1 = 4
define count_factors(x) {
auto i,c,fp[];
if(x<0)return count_factors(-x)+1;
if(x==0||x==1)return 0;
.=fac_store(fp[],x)
for(i=0;fp[i];i+=2)c+=fp[i+1]
return c;
}
# Count the number of unique prime factors of x
# e.g. 84 = [2]^2*[3]^1*[7]^1 -> #{2,3,7} = 3
define count_unique_factors(x) {
auto i,d,fp[];
if(x<0)return count_unique_factors(-x)+1;
if(x==0||x==1)return 0;
.=fac_store(fp[],x);
for(i=0;fp[i];i+=2).=d++
return d;
}
# Determine whether x is y-th power-free
# e.g. has_freedom(51, 2) will return whether 51 is square-free
# + sign of result indicates (-1)^[number of prime factors]
# making has_freedom(x,2) equivalent to the mobius function
# Special case: has_freedom(x, 1) returns whether x is prime
# Pseudo-boolean, since always returns 0 for false, but not always 1 for true
define has_freedom(x,y) {
auto os,i,p,m,fp[];
os=scale;scale=0;if(x<0)x=-x
if(x!=x/1){scale=os;return 0}
if(x==1){scale=os;return 1}
if(y==1){scale=os;return is_prime(x)}
if(x>A^A)if(is_prime(x)){scale=os;return -1}
if(x<0||y<1){scale=os;return 0}
.=fac_store(fp[],x);
m=1
for(i=1;p=fp[i];i+=2){
if(p>=y){scale=os;return 0}
m*=p*(-1)^p
}
return m
}
# Returns 0 if non-squarefree,
# 1 for 1 and all products of an even number of unique primes
# -1 otherwise
define mobius(x) {
return has_freedom(x,2);
}
# Determine the so-called arithmetic derivative of x
define arithmetic_derivative(x) {
auto os,i,f,n,d,fp[];
if(x<0)return -arithmetic_derivative(-x)
os=scale;scale=0;x/=1
if(x==0||x==1){scale=os;return 0}
.=fac_store(fp[],x);n=0;d=1
for(i=0;f=fp[i];i+=2){n=n*f+d*fp[i+1];d*=f}
n=(n*x)/d
scale=os;return n
}
### Storage and output of divisors of a number ###
# Count the number of divisors of x (prime or otherwise)
define count_divisors(x) {
auto i,c,p,fp[];
i=scale;x/=1;scale=i
if(x<0)return 2*count_divisors(-x);
if(x==0||x==1)return x
.=fac_store(fp[],x);
c=1;for(i=1;p=fp[i];i+=2)c*=1+p
return c
}
# Workhorse function for the below
define divisors_sp_(*divs[],x,print_) {
# opts: 1 -> print; 0 -> store
auto os,s,sx,c,ch,p,m,i,j,k,f,fp[]
os=scale;scale=0;x/=1
if(x==0||x==1){
if(!print_){divs[0]=x;divs[1]=0}
scale=os;return x;
}
.=fac_store(fp[],x)
c=1;for(i=1;(p=++fp[i++])>1;i++)c*=p
.=c--
s=1;if(x<0){s=-1;x=-x}
if(print_){
print s,", "
} else {
divs[0]=s
sx=0
ch=(c+1)/2
if(c>max_array_){
sx=sqrt(x)
print "too many divisors to store. storing as many as possible\n"
}
}
for(i=1;i<c;i++){
j=i;k=0;f=1
while(j){if(m=j%(p=fp[k+1]))f*=fp[k]^m;j/=p;k+=2} # generate a divisor
if(print_){print s*f,", "}else{
if(sx){
# Only applies if all divisors won't fit in the array
# Divisors <= sqrt(x) can be used to reconstruct missing divisors
# These can be overlooked due to the way the generator works
k=f;if(k<0)k=-k
if(k>sx){
# Store divisors > sqrt(x) in any space that would otherwise have
# been available at the high end of the array or else skip them
if(ch<max_array_)divs[ch++]=s*f
continue;
}
}
if(i<=max_array_){divs[i]=s*f}else{break}
}
}
if(print_){x}else{if(c>max_array_-1)c=max_array_-1;divs[c]=s*x;divs[c+1]=0}
scale=os;return s*x
}
# Displays all divisors of x in a logical but unsorted order
define divisors_print( x) { auto d[]; return divisors_sp_(d[],x,1); }
# Store calculated divisors in given array in same logical but unsorted order
# . see array.bc for sorting and rudimentary printing of arrays
define divisors_store(*d[],x) { return divisors_sp_(d[],x,0); }
# Returns the sum of the divisors of x where each divisor is raised
# to the power n; e.g. sigma(2,6) -> 1^2 + 2^2 + 3^2 + 6^2 = 50
# . only supports integer n at present
define sigma(n,x) {
auto os,i,u,d,f,fp[];
os=scale;scale=0;x/=1;n/=1
if(n==0){scale=os;return count_divisors(x)}
if(n<0){scale=os;n=-n;return sigma(n,x)/x^n}
.=fac_store(fp[],x);u=d=1
if(x<0)x=0
if(x==0||x==1)return x;
for(i=0;fp[i];i+=2){u*=(f=fp[i]^n)^(fp[i+1]+1)-1;d*=f-1}
u/=d;scale=os;return u;
}
# Old slow version of sigma
# supports integer and half-integer n
#define sigma_slow(n,x) {
# auto os,c,p,m,h,i,j,k,f,sum,fp[];
# if(n==0)return count_divisors(x);
# n+=n
# os=scale;scale=0;x/=1;n/=1
# if(x<0)x=0
# if(x==0||x==1){scale=os;return x;}
#
# p=h=n;n/=2;h-=n+n
# if(p<0){scale=os;n=-n;sum=sigma(-p/2,x)/x^n;if(h)sum/=sqrt(x);return sum}
# .=fac_store(fp[],x)
# c=1;for(i=1;(p=++fp[i++])>1;i++)c*=p
# .=c--;p=x^n;if(h){scale=os;p*=sqrt(x);scale=0};sum=1+p
# for(i=1;i<c;i++){
# j=i;k=0;f=1
# while(j){if(m=j%(p=fp[k+1]))f*=fp[k]^m;j/=p;k+=2}
# p=f^n;if(h){scale=os;p*=sqrt(f);scale=0}
# sum+=p
# }
# scale=os;return sum
#}
# Returns the sum of the divisors of x, inclusive of x
# e.g. primes -> prime + 1, 2^x -> 2^(x+1)-1, 6 -> 12
define sum_of_divisors(x) {
return sigma(1,x);
}
# Computes the Euler totient function for x
define totient(x) {
auto i,t,f,fp[];
.=fac_store(fp[],x);t=1
if(x==0||x==1)return x;
for(i=0;fp[i];i+=2)t*=((f=fp[i])-1)*f^(fp[i+1]-1)
return t
}
# Number of iterations of the totient function to reach 1
define totient_itercount(x) {auto t;while(x>1){x=totient(x);.=t++};return t}
# Sum of iterating the totient function down to 1
define totient_itersum(x) {auto t;while(x>1){x=totient(x);t+=x};return t}
# Returns if x is a perfect totient number
define is_perfect_totient(x) { return totient_itersum(x)==x }
# Computes Ramanujan's c_q function for n (given q)
define ramanujan_c(q,n) {
auto i,c,d,t,f,p,fp[];
if(q<0||n<0)return 0;
if(q==1)return 1;
if(n==1)return mobius(q);
if(n==q)return totient(q);
n=q/int_gcd(q,n)
if(n==1)return totient(q);
.=fac_store(fp[],n);t=1;c=d=0;
for(i=0;fp[i];i+=2){
t*=((f=fp[i])-1)*f^((p=fp[i+1])-1)
c+=p;.=d++
}
if(c!=d){c=0}else{c=(-1)^d}
return c*totient(q)/t # mobius(n')*totient(q)/totient(n')
}
# Determines whether a number is a practical number
define is_practical(x) {
auto os,i,ni,s,p,f,nf,fp[];
if(x<0)return 0;
if(x==1)return 1;
os=scale;scale=0;i=x%2;scale=os
if(i)return 0
.=fac_store(fp[],x)
if(!fp[2])return 1;# x is power of 2
scale=0
#fp[0] has to be 2 here, so replace with 2
s=2^(fp[1]+1)-1
f=fp[i=2]
while(1){
if(f>1+s){scale=os;return 0}
if((nf=fp[ni=i+2])==0){scale=os;return 1}
s*=(f^(fp[i+1]+1)-1)/(f-1)
f=nf;i=ni
}
return -1;# should never get here
}
### Exploring prime neighbourhood ###
# Finds and returns the nearest prime > x
define nextprime(x) {
auto os,ox;
if(x<0)return -prevprime(-x)
if(x<2)return 2
if(x<3)return 3
os=scale;scale=0
ox=x
if(x/1<x).=x++ #ceiling
x/=1 #truncate
x+=1-x%2 #make odd
if(x==ox)x+=2 #same as we started with
#while(!is_prime(x))x+=2 # could use jumps here, but is not much faster
while(1){
while(!is_small_division_pseudoprime(x))x+=2
if(x<A^A)break; # pairs with the A^5 in is_s.d.pp()
if(is_rabin_miller_pseudoprime(x)){
if(rabin_miller_maxtests_){
if(is_perrin_pseudoprime(x))break;
} else {
break;
}
}
x+=2
}
scale=os;return x
}
# nearest prime >= x
define nextprime_ifnotprime(x) {
if(is_prime(x))return x;
return nextprime(x)
}
# Finds and returns the nearest prime < x
define prevprime(x) {
auto os,ox;
if(x<0)return -nextprime(-x)
if(x<=2)return -2
if(x<=3)return 2
if(x<=5)return 3
os=scale;scale=0
ox=x;x/=1
if(x%2){if(x==ox)x-=2}else{.=x--}
#while(!is_prime(x)&&x>0)x-=2
while(x>0){
while(!is_small_division_pseudoprime(x))x-=2
if(x<A^A)break; # pairs with the A^5 in is_s.d.pp()
if(is_rabin_miller_pseudoprime(x)){
if(rabin_miller_maxtests_){
if(is_perrin_pseudoprime(x))break;
} else {
break;
}
}
x-=2
}
if(x<2)return x=-2
scale=os;return x
}
# nearest prime <= x
define prevprime_ifnotprime(x) {
if(is_prime(x))return x;
return prevprime(x)
}
# Find the nearest prime to x (inclusive)
define nearestprime(x) {
auto np,pp;
if(is_prime(x))return x;
np=nextprime(x)
pp=prevprime(x)
if(np-x>x-pp)return pp
return np
}
### Relatives of the Primorial / Factorial
# Primorial
define primorial(n) {
auto i,pm,p;
pm=1;p=2
if(prime[max_array_])for(i=1;i<=max_array_&&p=prime[i]<=n;i++)pm*=p
for(p=p;p<=n;p=nextprime(p))pm*=p
return pm
}
# Subprimorial
define subprimorial(n) {
auto i,pm,p;
pm=1;p=2
if(prime[max_array_])for(i=2;i<=max_array_&&(p=prime[i])<=n;i++)pm*=p-1
for(p=p;p<=n;p=nextprime(p))pm*=p-1
return pm
}