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
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; |
|
}
|
|
|