/* BC program primepatterns 26-28th March 2014 */
pglobal[0]=2
pglobal[1]=3
pglobal[2]=5
pglobal[3]=7
pglobal[4]=11
pglobal[5]=13
pglobal[6]=17
pglobal[7]=19
pglobal[8]=23
pglobal[9]=29
pglobal[10]=31
pglobal[11]=37
pglobal[12]=41
pglobal[13]=43
pglobal[14]=47
pglobal[15]=53
pglobal[16]=59
pglobal[17]=61
pglobal[18]=67
pglobal[19]=71
nglobal=19 /* this is one less than the number of primes in the above list. */

/* See K. Prachar, Primzahlverteilung, Satz 4.2. p. 45 for the concept of admissible sequence of 
 * primes a[0]m+b[0],...,a[s-1]m+b[s-1]. It was conjectured by L.E. Dickson (1904) that these 
 * numbers will be simultaneously prime for infinitely many m.  See Ribenboim, The New Book of 
 * Prime Number Records, p.372.
 */
define testpattern(a[],b[],s){
auto i,k,count,l,p,m
  for(i=0;i<s;i++){
     if(a[i]<=0){
       print "a[",i,"]<=0\n"
       return(0)
     }
     if(gcd(a[i],b[i])>1){
       print "gcd(a[",i,"],b[",i,"])>1\n"
       return(0)
     }
     for(k=0;k<s && k!=i;k++){
        if(a[k]==a[i] && b[k]==b[i]){
          print "a[",k,"]=a[",i,"] and b[",k,"]=b[",i,"]\n"
          return(0)
        }
     }
  }
  product=1
  for(l=0;l<=nglobal && pglobal[l]<=s;l++){
     p=pglobal[l]
     count=0
     for(m=0;m<p;m++){
       count=count+poly(a[],b[],s,m,p)
     }
     if(count==p){
       print "omega(",p,")=",p,"\n"
       return(0)
     }
  }
  print "we have an admissible constellation "
  for(i=0;i<s;i++){
   t=printp(a[i],b[i])
   print ","
   }
   print "\n"
  return(1)
}

define poly(a[],b[],s,m,p){
auto prod,i
  prod=1
  for(i=0;i<s;i++){
    prod=prod*(a[i]*m+b[i])
  }
  prod=prod%p
  if(prod==0){
    return(1)
  }else{
    return(0)
  }
}

/* This is the main program here. It requires BC programs gcd, jacobi and a file containing the 
 * a[0],...,a[s-1] and b[0],...,b[s-1]. Here 2<=s<=20.
 * Here 1<=u<=v. The a[i] and b[] are tested to be an admissible sequence, in which case 
 * we replace u by the least m>=1 such that all of a[i]m+b[i] are > 1.
 * Then the numbers a[i]m+b[i], 0<=i<s, are tested for pseudoprimality by the Lucas-strong base 2  
 * pseudoprime test. 
 */

define primepatterns(a[],b[],s,u,v){
 auto t,i,m,j,temp,count,flag
 if(s>20){
    print "number of linear forms is greater than 20\n"
    return(-1)
 }
 if(testpattern(a[],b[],s)==0){
    print "inadmissible test pattern\n"
    return(-1)
 }
 if(u<1||v<1){
   print "u < 1 or v< 1\n"
   return(-1)
 }
 if(u>v){
   temp=v
   v=u
   u=temp
 }
 m=1
 while(1){
 for(i=0;i<s;i++){
   if(a[i]*m+b[i]<=1){
     m=m+1
     break
   }
 }
 if(i==s){
   u=max(u,m)
   break
 }
}
 count=0
 for(m=u;m<=v;m++){
   for(i=0;i<s;i++){
     temp=a[i]*m+b[i]
     if(lucas0(temp)==0){
       break
     }
   }
   if(i==s){
       count=count+1
  /*   print "we have found ",s," primes\n" */
     t=s-1
     for(j=0;j<t;j++){
       temp=a[j]*m+b[j]
       print temp,","
     }
     temp=a[j]*m+b[j]
     print temp,"\n"
   }
 }
 print "The number of prime ",s," constellations "
 for(i=0;i<s;i++){
   t=printp(a[i],b[i])
   print ","
 }
 print "\n"
 print "in the range \n"
 print u, "<=m<=",v,"\n"
 print "is "
 return(count)
}

/* This prints am+b in a nice way. */
define printp(a,b){
  if(a>1){
    if(b>0){
      print a,"m+",b
    }
    if(b<0){
      print a,"m",b
    }
    if(b==0){
      print a,"m"
    }
  }
  if(a==1){
    if(b>0){
      print "m+",b
    }
    if(b<0){
      print "m",b
    }
    if(b==0){
      print "m"
    }
  }
  if(a<-1){
    if(b>0){
      print a,"m+",b
    }
    if(b<0){
      print a,"m",b
    }
    if(b==0){
      print a,"m"
    }
  }
  if(a==-1){
    if(b>0){
      print "-m+",b
    }
    if(b<0){
      print "-m",b
    }
    if(b==0){
      print "-m"
    }
  }
}