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.

152 lines
4.4 KiB

#!/usr/local/bin/bc -l primes.bc primes_db_minipack.bc
### Primes_DB_Code.BC - prime counting / retrieval functions reliant on a DB
## Technical details:
## . Modulo 2310 [=primorial(11)] there are 480 possible primes
## . Each array entry in the pd_[] packed prime data array should represent
## . the primality of the 480 possibilities within the integers from
## . 2310*index to 2310*(index+1)-1
## . The last 24 bits of each array entry should be the cumulative count of
## . . actual primes up to the end of that entry.
## . Primes 2 to 11 are handled seperately as they are not relatively prime
## . . to 2310
## . This file and its minipack and bigpack siblings are designed to replace
## . . the older, unpacked, primes_db.bc. The prime[] array from the latter
## . . can be recreated using the fillprimearray() function herein.
max_array_ = 4^8-1
# Create metadata arrays which allow the below to operate correctly
define makemods2310_() {
auto oib,c,i;
oib=ibase;ibase=A
if(!pd_2310_[239]){c=0;for(i=1;i<=1155;i++)if(int_gcd(i,2310)==1)pd_2310_[c++]=i}
if(!pd_0112_[max_array_]){for(c=1;c<=32768;c+=c)for(i=c;i<c+c;i++)pd_0112_[i]=1+pd_0112_[i-c]}
ibase=oib;return 0
}
# Return the nth prime up to the limit of the database
# . TO DO: slow follow-on using nextprime()
define prime(n) {
auto oib,os,ma,i,l,m,h,k,hk,c,cm,b;
if(n<0)return -prime(-n)
if(n<1)return 1;if(n<2)return 2;if(n<3)return 3
if(n<4)return 5;if(n<5)return 7;if(n<6)return B
if(n<=max_array_&&k=prime[n]){return k} # deliberate assignment!
oib=ibase;ibase=A
os=scale;scale=0
n-=6;n/=1
h=pd_max_;cm=16777216#=16^6
if(n>=(c=pd_[h]%cm)){
print "prime() knows only the first ",c+5," primes\n";
scale=os;return -1
}
l=0;m=(l+h)/2;.=makemods2310_()
while(l<=m){
k=pd_[m];hk=k/cm;c=k-hk*cm;k=hk
if(c<=n){
l=m+1
} else {
b=0;if(m>0)b=pd_[m-1]%cm
if(b>n){
h=m
} else {
# main routine
ma = 1+max_array_
i=-1;while(k&&b<=n){
if(n-b<16){
hk=k/2;if(k-hk-hk).=b++;.=i++;k=hk
} else {
hk=k/ma;b+=pd_0112_[k-ma*hk];i+=16;k=hk # addition to i should be log2(max_array_+1)
}
}
if(i<240){
m= m *2310+pd_2310_[ i]
} else {
m=(m+1)*2310-pd_2310_[479-i]
} # end if i
scale=os;ibase=oib;return m
} # end if b
} # end if c
m=(l+h)/2
} # end while
scale=os;ibase=oib;return -1
}
# Return the number of primes up to (and including) p
# . as long as p is less than or equal to the maximum db entry
define primepi(p){
auto oib,os,ma,q,f,i,b,l,m,h,k;
if(p<0)return -primepi(-p);
p=prevprime_ifnotprime(p);
if(p<2)return 0;if(p<3)return 1;if(p<5)return 2
if(p<7)return 3;if(p<B)return 4;if(p<D)return 5
oib=ibase;ibase=A
os=scale;scale=0
if(q=prime[max_array_]&&p<=q){ # assignment of q is deliberate
l=1;h=max_array_;m=(max_array_+1)/2
while(p!=(q=prime[m])){
if(p>q){l=m+1}else{h=m-1}
m=(l+h)/2
}
ibase=oib;scale=os;return m
}
q=(1+pd_max_)*2310-1
if(p>q){
#q=prevprime_ifnotprime(q)
print "primepi() only knows the count of primes up to ";q
ibase=oib;scale=os;return -1
}
.=makemods2310_()
i=p/2310; q=p-i*2310
f=0;if(q>1155){f=1;q=2310-q}
#find q in pd_2310_
l=0;h=239;m=120
while(q!=(k=pd_2310_[m])){
if(q<k){
h=m-1;
} else {
l=m+1;
}
m=(l+h)/2
}
if(f)m=479-m
b=0;if(i)b=pd_[i-1]%16777216
k=pd_[i]/16777216;ma=1+max_array_
while(m){
if(m<16){
h=k/2;if(k-h-h).=b++;.=m--;k=h
} else {
h=k/ma;b+=pd_0112_[k-ma*h];m-=16;k=h # subtraction from m should be log2(max_array_+1)
}
}
ibase=oib;scale=os;return b+6;
}
# Unpack the first max_array_ primes from the global pd_ array into the
# . global prime array which is used by some functions in primes.bc
# . to speed up trial division code.
define fillprimearray() {
auto oib,os,ma,i,j,k,h,mi,s,ip
oib=ibase;ibase=A
os=scale;scale=0
.=makemods2310_()
prime[0]=1;prime[1]=2;prime[2]=3
prime[3]=5;prime[4]=7;prime[5]=B
ip=6;mi=0;ma=1+max_array_
for(i=0;i<356;i++){
k=pd_[i]/16777216;s=1
for(j=0;k&&j<240;j++){
h=k/2;if(k-h-h)prime[ip++]=mi+pd_2310_[j];k=h
if(ip==ma){i=ip;break}
}
mi+=2310
for(j=239;k&&j>=0;j--){
h=k/2;if(k-h-h)prime[ip++]=mi-pd_2310_[j];k=h
if(ip==ma){i=ip;break}
}
}
ibase=oib;scale=os;return 0
}