/* bc program nprimeap */ /* needs lucas (which needs jacobi) */ print "\n" print "Finds the first Lucas probable prime p >= m, p=ak+b. \n" print "Here 00){ a=b b=c c=a%b /* c=r[j]=r[j-2] mod(r[j-1]) */ } return(b) } define nprimeap(m,a,b){ auto flag,i,h,r,pp,t,y0 y0=9591 if(b<=0){ print "b<=0\n" return(0) } if(b>=a){ print "b>=a\n" return(0) } if (gcd(a,b)>1){ print "gcd(a,b)>1\n" return(0) } if(m<=p[y0]){ for(i=0;i<=y0;i++){ if(m<=p[i] && p[i]%a==b){ print p[i]," is the first prime >= ",m," and of the form ",a,"k+",b,"\n" return(1) } } } /* m > prime[y0] */ r=m%a q=m+b-r if(b= ",m," and of the form ",a,"k+",b,"\n" return(1) }else{ q=q+a } } }