最速の素数判定アルゴリズム　　

function　LucasTest(n:int64):boolean;
var dt,odt:array of Word;
var cy:int64;
var i,hDt,hPg:integer;
var mask:LongWord;

　procedure　mul2; //２乗する処理
　var w1,w2:array of Word;
　var wk:LongWord;
　var i,j:integer;
　begin
　SetLength(w1,1+hDt );//メモリを確保して
　　SetLength(w2,1+1+hDt );
　　 cy:=0;
　　for i:=0　to　 hDt do　 begin
　　　for j:=0 to i do　cy:=cy+dt[j]*dt[i-j];
　　　 w1[i]:=cy;//下位桁だけを代入
　　　 cy:=cy shr 16;
　　end;
　　for i:=1　to　 hDt do　 begin
　　　for j:=i to hDt do　cy:=cy+dt[j]*dt[hDt+i-j];
　　　 w2[i-1]:=cy;//
　　　 cy:=cy shr 16;
　　end;
　　w2[hDt　]:=cy; cy:=cy shr 16;
　　w2[hDt+1]:=cy; cy:=0;
　　 wk:=(w2[0] shl (16- hPg)) and $ffff;
　　 cy:=cy+ wk +w1[0];
　　 dt[0]:=cy;//下位桁だけを代入
　　 cy:=cy shr 16;
　　for i:=1　to　 hDt do　begin //上位を下位に加算する
　　 wk:=(((w2[i] shl 16)+w2[i-1])shr hPg) and $ffff;
　　 cy:=cy+　wk　+w1[i];
　　 dt[i]:=cy;//下位桁だけを代入
　　 cy:=cy shr 16;
　　end;
　　 wk:=((w2[hDt+1] shl 16)+w2[hDt])shr hPg ;
　　cy:=wk+(( (cy shl 16) +dt[hDt]) shr hPg);
　　while cy >0 do　begin
　　　dt[hDt]:=dt[hDt] and mask;
　　　for i:=0　to　 hDt do　begin //さらに溢れた上位を加算する
　　　　　　　　cy:=cy+dt[i];
　　　　 dt[i]:=cy;//下位桁だけを代入
　　　　 cy:=cy shr 16;
　　　end;
　　cy:=( (cy shl 16)+ dt[hDt] ) shr hPg;
　　end;
　end;

----- つづく ----

var fsame,fzero:boolean;
var Rep:int64;
begin
hDt:=((n-1) div 16);//最大桁位置
hPg:= ((n-1) mod 16)+1 ;//ビット位置
mask:=(1 shl hPg)-1;
SetLength(dt,1+hDt );//メモリを確保して
for i:=0 to High(dt) do dt[i]:=0; dt[0]:=4; //初期値を４に
rep:=n;
while rep>0 do begin
　odt:=copy(dt,0,hdt+1);
　mul2; //２乗する
　fsame:=　　　　　dt[0]=odt[0];
　fzero:=true;
　for i:=1 to High(dt) do　begin
　　 fsame:=fsame and (dt[i]=odt[i]);//全て同じかチェック
　　 fzero:=fzero and (dt[i]=0); //全て同じかチェック
　　 if (not fsame) and (not fzero) then break;
　end;
　if fsame then begin Result:=false;exit;end;//もはや見込みなし
　cy:=-2;
　if fzero then
　begin
　　 if dt[0] =2　then begin Result:=True;exit;end;//オメデトウ
　　 if dt[0] <2　then begin //2以下なら負数になるから
　　　for i:=1 to　High(dt) do dt[i]:=$ffff;//全部が１は0と同じ
　　　　dt[hDt]:=mask;
　　　　dt[0]:=$ffff+dt[0]-2;cy:=0;　//で引算する
　　　end;
　　end;
　　for i:=0　to　 hDt do　begin //2を引く処理
　　　 cy:=cy+dt[i];
　　　 dt[i]:=cy;//下位桁だけを代入
　　　 cy:=(cy-dt[i]) div $10000;
　　　 if cy=0 then break;
　　 end;
　Dec(rep);
　end;
Result:=false;
end;