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 Raw Permalink Blame History

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