### Primes.BC - Primes and factorisation (rudimentary)

# funcs.bc is required to run this library

# Returns 2, 3, or number of form 6n[+-]1
define aq(x) {
 if(x<0)return(-aq(-x))
 if(x<3)return(x+1)
 x-=3;x+=int(x/2)
 return(x+x+5)
}

# Inverse of the above
define iaq(x) {
 if(x<0)return(-iaq(-x))
 if(x<4)return(x-1)
 return((rem(x+3,6)+x+x)/6+1)
}

# Returns 2, 3, 5 or number of form 30n[+-]{1,7,11,13}
define aq30(x) {
 auto os, e, r, rh, s
 os=scale;scale=0;x/=1
 if(x<0){x=-aq30(-x);scale=os;return(x)}
 if(x<4){x=x+1+x/3  ;scale=os;return(x)}
 x-=3     ; e=x/8
 r=x%8    ; rh=r/4
 s=1-2*rh ; r=s*r+7*rh
 scale=os;return( 30*(e+rh)-s*(r^2-7*r-1) )
}
 
# Inverse of the above
define iaq30(x) {
 auto os, e, r
 os=scale;scale=0;x/=1
 if(x<0){x=-iaq30(-x);scale=os;return(x)}
 if(x<7){x=x-1-x/5   ;scale=os;return(x)}
 e=x/30;r=x%30
 #r+=(31-r)/30#not needed with bc integer division
 r=r/6+(r-2)/7
 scale=os;return(8*e+r +3)
}

# Cyrek's Approximation to the Prime Counting Function pi(x)
define aprimepi(x) { 
 auto a,b,lx,k
 lx=l(x)/l(10)
 a=pow(lx,2/(3*lx))
 b=1+2/(17*pow(ch(lx-e(2)),1/32))
 k=1+l(a/b)
 return(x*k/l(x))
}

# Output the prime factors of x
define fac(x) {
 auto os,j[],ji,p,n,c,xx;#oldscale,jump,jump-index,prime,nth,count,xx
 os=scale;scale=0;x/=1
 if(x<0){print"-1*\\\n";x*=-1}
 if(x==0){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){
  xx=x/p
  c=0;while(x==p*xx){c+=1;x=xx;xx=x/p}
  if(c){print p;if(c>1){print "^",c};print "*\\\n"}
 }
 ji=7;p=7
 while(p^2<=x){
  xx=x/p
  c=0;while(x==p*xx){c+=1;x=xx;xx=x/p}
  if(c){print p;if(c>1){print "^",c};print "*\\\n"}
  p+=j[ji=(ji+1)%8]
 }
 if(x>1)print x,"*\\\n"
 scale=os;return(1)
}

# Determine whether x is prime
define is_prime(x) {
 auto os,j[],ji,p,n,xx;#oldscale,jump,jump-index,prime,nth,xx
 os=scale;scale=0;x/=1
 if(x<0){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){return 1};xx=x/p;if(x==p*xx){return 0}}
 ji=7;p=7
 while(p^2<=x){
  xx=x/p;if(x==p*xx){return 0}
  p+=j[ji=(ji+1)%8]
 }
 scale=os;return(1)
}

# Determine whether x is y-th power-free
#  e.g. has_freedom(51, 2) will return whether 51 is square-free
#  Special case: has_freedom(x, 1) returns whether x is prime
define has_freedom(x,y) {
 auto os,j[],ji,p,py,n,xx;#oldscale,jump,jump-index,prime,nth,xx
 os=scale;scale=0;x/=1
 if(y==1){return is_prime(x)}
 if(x<0||y<1){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){return 1};py=p^y;xx=x/py;if(x==py*xx){return 0}}
 ji=7;p=7;py=p^y
 while(py<=x){
  xx=x/py;if(x==py*xx){return 0}
  p+=j[ji=(ji+1)%8];py=p^y
 }
 scale=os;return(1)
}

# Fill the global prime[] array with g primes
define genprimes(g) {
 auto os,j[],ji,i,x,xx,p,p2,pflag
 #os-oldscale,j[]-jumps,ji-jumpindex,i-primeindex,x-primecandidate
 #xx-x/p,p-prime,pflag-is.x.prime?
 os=scale;scale=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
 prime[0]=1;prime[1]=2;prime[2]=3;prime[3]=5;prime[4]=7;primetop=4
 if(g<100)g=100;if(g>65535)g=65535
 for(i=primetop+1;i<=g;i++)prime[i]=0
 ji=0;x=11
 while(primetop<g){
  i=2;p=prime[i];pflag=1
  while(p^2<=x){
   xx=x/p;if(x==p*xx){pflag=0;break}
   i+=1;p=prime[i]
  }
  if(pflag){primetop+=1;prime[primetop]=x;}#print primetop,": ",x,"\n"}
  x+=j[ji=(ji+1)%8]
 }
 scale=os;return(prime[primetop])
}