bc_libs/primes.bc

#!/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
}