From a5fdb9121c31938559a9df7174d97965dbcbff0a Mon Sep 17 00:00:00 2001
From: Mark Wooding <mdw@distorted.org.uk>
Date: Mon, 25 Nov 2019 13:12:43 +0000
Subject: [PATCH] catacomb/__init__.py: Rewrite `kcdsaprime' to return the
 correct length.

The old version wouldn't return numbers of the correct length reliably.
This new one, written to follow Catacomb's own algorithm, will.
---
 catacomb/__init__.py | 22 ++++++++++++++--------
 t/t-pgen.py          |  2 +-
 2 files changed, 15 insertions(+), 9 deletions(-)

diff --git a/catacomb/__init__.py b/catacomb/__init__.py
index 08ec3d7..2b433cf 100644
--- a/catacomb/__init__.py
+++ b/catacomb/__init__.py
@@ -1186,13 +1186,19 @@ def findprimitive(mod, hh = [], exp = None, name = 'g', event = pgen_nullev):
 
 def kcdsaprime(pbits, qbits, rng = rand,
                event = pgen_nullev, name = 'p', nsteps = 0):
-  hbits = pbits - qbits
-  h = pgen(rng.mp(hbits, 1), name + ' [h]',
-           PrimeGenStepper(2), PrimeGenTester(),
-           event, nsteps, RabinMiller.iters(hbits))
-  q = pgen(rng.mp(qbits, 1), name, SimulStepper(2 * h, 1, 2),
-           SimulTester(2 * h, 1), event, nsteps, RabinMiller.iters(qbits))
-  p = 2 * q * h + 1
-  return p, q, h
+  hbits = pbits - qbits - 1
+  while True:
+    h = pgen(rng.mp(hbits, 1), name + ' [h]',
+             PrimeGenStepper(2), PrimeGenTester(),
+             event, nsteps, RabinMiller.iters(hbits))
+    while True:
+      q0 = rng.mp(qbits, 1)
+      p0 = 2*q0*h + 1
+      if p0.nbits == pbits: break
+    q = pgen(q0, name, SimulStepper(2*h, 1, 2),
+             SimulTester(2 * h, 1), event, nsteps, RabinMiller.iters(qbits))
+    p = 2*q*h + 1
+    if q.nbits == qbits and p.nbits == pbits: return p, q, h
+    elif nsteps: raise ValueError("prime generation failed")
 
 #----- That's all, folks ----------------------------------------------------
diff --git a/t/t-pgen.py b/t/t-pgen.py
index 651ce70..e200add 100644
--- a/t/t-pgen.py
+++ b/t/t-pgen.py
@@ -221,7 +221,7 @@ class TestPGen (U.TestCase):
     me.assertEqual(p, 2*q*h + 1)
     me.assertTrue(p.primep())
     me.assertEqual(p.nbits, 512)
-    me.assertEqual(ev.events, "[p [h]:F17/P6/D][p:F60/P26/D]")
+    me.assertEqual(ev.events, "[p [h]:F53/P6/D][p:F32/P26/D]")
 
 ###--------------------------------------------------------------------------
 class TestPrimeIter (U.TestCase):
-- 
2.11.0