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.

117 lines
2.6 KiB

#!/usr/local/bin/bc -l primes.bc
### Primes-Other.BC - Extra functions to go along with Primes.BC
# Both Funcs.BC and Primes.BC are required to use functions herein
# 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((remainder(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( 3*A*(e+rh)-s*(r*(r-7)-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=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 la,b,lx,k,oib;
if(x<=0)return 0
if(x<A){return x*aprimepi(A)/A}
oib=ibase;ibase=A;scale+=4
lx=l(x)/2.3026 #l(10)
la=2*l(lx)/(3*lx)
#b=1+2/(17*pow(cosh(lx-e(2)),1/32))
b=e(lx-e(2))
b=(b+1/b)/2 #cosh b
b=17*sqrt(sqrt(sqrt(sqrt(sqrt(b))))) #17.b^(1/32)
b=2/b+1
k=1+la-l(b)
ibase=oib;scale-=4
return(x*k/l(x))
}
# Use the above approximation to find a
# number close to the nth prime
define fastguessprime(n) {
auto os,l,x,ox,i;
os=scale;scale=10
s=1;if(n<0)n*=(s=-1)
n+=.5;l=l(n);x=n*l
ox=1
while(ox!=i){
ox=i;x+=l*(n-aprimepi(x))
scale=0;i=x/1;scale=10
}
scale=os;return ox
}
# Use the above to find a prime close to the nth prime
# (is almost always wrong, but is generally within 0.5%)
define guessprime(n) {
n=fastguessprime(n)
return nearestprime(n)
}
# Sum of prime factors of a number
# . e.g. 150 = 2*3*5^2 -> 2+3+5*2 = 15
define sum_of_factors(x) {
auto i,c,fp[];
if(x<0)return sum_of_factors(-x)-1;
if(x==0||x==1)return 0;
.=fac_store(fp[],x)
for(i=0;fp[i];i++)c+=fp[i]*fp[++i]
return c;
}
# As above but with no splitting of powers into multiplies
# . e.g. 150 = 2*3*5^2 -> 2+3+5^2 = 30
define sum_of_factor_terms(x) {
auto i,c,fp[];
if(x<0)return sum_of_factor_terms(-x)-1;
if(x==0||x==1)return 0;
.=fac_store(fp[],x)
for(i=0;fp[i];i++)c+=fp[i]^fp[++i]
return c;
}
# Raise the powers of the prime factors to the
# power of their primes and multiply
# . e.g. 150 = 2*3*5^2 -> 1^2*1^3*2^5 = 1*1*32 = 32
define factor_invert(x) {
auto i,c,fp[];
if(x<0)return factor_invert(-x)+1;
if(x==0||x==1)return 0;
.=fac_store(fp[],x)
c=1;for(i=0;fp[i];i+=2)c*=fp[i+1]^fp[i]
return c;
}