Looking for primes
Math 302 Advanced Cryptography
March 8, 2019

T. Ratliff, Wheaton College
ratliff_thomas@wheatoncollege.edu

mrWitness[a, n] determines if a is a Miller-Rabin witness for n

mrWitness[a_,n_]:=Module[{temp=n-1,k=0,q,aq,i},

While[Mod[temp,2]==0,temp=temp/2;k++];
q=(n-1)/2^k;
aq=PowerMod[a,q,n];

If[Mod[aq,n]==1,Return[False]];

i=0;
While[i<k,If[Mod[aq+1,n]==0,Return[False]];
aq=Mod[aq^2,n];
i++];

Return[True]
]

mrWitness[2,561]

mrWitness[2,409]

mrTest[n, r] tests if n is a pseudoprime by checking r random values as MR witnesses for n
	Returns False if n detected as composite, and True none of r values are MR witnesses

mrTest[n_,r_]:=Module[{a,i},
If[n==2,Return[True]];
If[EvenQ[n],Return[False]];

For[i=1,i≤r,i++,
a=RandomInteger[{1,n-1}];

If[GCD[a,n]>1,Return[False]];

If[mrWitness[a,n]==True,Return[False]]
];

Return[True]];

mrTest[243242342432421,10]

kbitPrime[k] returns a k-bit number where 30 random values failed to be a MR witness

kbitPrime[k_]:=Module[{n,i=1},
n=2*RandomInteger[ {2^(k-2),2^(k-1)}]+1;

While[mrTest[n,30]==False, 
n=2*RandomInteger[ {2^(k-2),2^(k-1)}]+1;
i++];
Print["   Required ",i," random choices to find prime"];
Return[n]
]

kbitPrime[160]

kbitPrime[1024]

kbitPrime[4096]

kbitPrime[4096]

parametersDSA[qbits, pbits] returns a pair {q,p} where 
	q is a qbits-bit prime, 
	p is a pbits-bit prime, 
	and q | (p-1)
	
By default, checks 30 random values for MR witnesses

parametersDSA[qbits_,pbits_ ]:=Module[ {q,p,m,i= 0,j=0,foundp=False}, 
While[ foundp==False, 
Print["Looking for q. . ."];
q=kbitPrime[qbits];
j++;
For[i=1,(i≤ 4096)&& (foundp==False), i++,
m=RandomInteger[{2^(pbits-1),2^pbits}];
p=m-Mod[m,2*q]+1;

If[mrTest[p,30]==True,
foundp=True];
];
];
Print["Found {q,p} on choice #",j," of q and choice #",i," of p for this q"];

Return[{q,p}];
]

parametersDSA[160,1024]

parametersDSA[160,1024]

Find an element of order q in (Subscript[F, p]^( * ))

q=760846931543193671156140389968568795335166685423;
p=;

Let’s verify that {q,p} have the desired properties

{Length[IntegerDigits[q,2]], Length[IntegerDigits[p,2]]}
PrimeQ[{q,p}]
Mod[p-1,q]

b=PowerMod[2, (p-1)/q, p]

PowerMod[b,q,p]

Since (b^q)=1 mod p, we have ord(b) | q. 
	Since q is prime, ord(b)=1 or ord(b)=q.
	But b ≠ 1 so ord(b)=q