#!/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(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(x7)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 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;isx){ # 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(chmax_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 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 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(xx-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 }