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