########################
# SCRIPT: dh_gen.bdscr #
########################
# Generate domain parameters (p,q,g)
# with primes (p,q) where q is a prime divisor of p-1; p = jq+1 with j>=2;
# and g has order q mod p; i.e. g^q mod p = 1 if g!=1
L=1024
N=160
println("Generating p,q for L=", L," N=",N)
q=genprime(N)
println('q=',hex(q))
# Note that p = x - (x mod (2*q)) + 1 => p is congruent to 1 mod 2q => 2q | (p-1)
i = 0
repeat 
  i++; 
  println('Iteration for p ',i," "); 
  x=randbits(L); 
  p=x-(x mod (2*q))+1; 
until (bitlen(p)==L and isprime(p))
# Generate generator, g
puts "Generating g..."
e = (p-1)/q
printf("e=%#x\n", e)
i = 0
repeat 
  i++;
  print('Iteration for g ',i, " ");
  h=random(p-1);
  g=modexp(h,e,p);
until (g ne 1)
# Output results and checks...
printf("\nFinal p=%#x\n",p)
println("bitlen(p)=",bitlen(p))
println("``p is prime`` is ", bool(isprime(p)))
printf("q=%#x\n",q)
println("bitlen(q)=",bitlen(q))
println("``q is prime`` is ", bool(isprime(q)))
println("(p-1) mod q = ", (p-1) mod q, " (expecting 0)")
printf("g=%#x\n",g)
println("``g < p`` is ", bool((g < p))) 
println("g^q mod p = ", modexp(g,q,p), " (expecting 1)")
save  # Save (p,q,g) variables for later