Add qsieve and qsafe, two programs used for generating the moduli file.

These were written by William Allen Simpson and submitted in PR 21983,
and are added with minor adjustments and nits from christos@ and myself.

Approved by christos@ and groo@.

9 files changed, 1105 insertions, 0 deletions

diff --git a/moduli/Makefile b/moduli/Makefile
new file mode 100644
index 00000000..4e8099bd
--- /dev/null
+++ b/moduli/Makefile
@@ -0,0 +1,5 @@
+# $NetBSD: Makefile,v 1.1 2006/01/19 23:23:58 elad Exp $
+
+SUBDIR=	qsieve qsafe
+
+.include <bsd.subdir.mk>
diff --git a/moduli/Makefile.inc b/moduli/Makefile.inc
new file mode 100644
index 00000000..a5346c37
--- /dev/null
+++ b/moduli/Makefile.inc
@@ -0,0 +1,7 @@
+# $NetBSD: Makefile.inc,v 1.1 2006/01/19 23:23:58 elad Exp $
+
+.include <bsd.own.mk>
+
+WARNS=4
+DPADD+= ${LIBCRYPTO}
+LDADD+= -lcrypto
diff --git a/moduli/qsafe/Makefile b/moduli/qsafe/Makefile
new file mode 100644
index 00000000..ab8484ae
--- /dev/null
+++ b/moduli/qsafe/Makefile
@@ -0,0 +1,10 @@
+# $NetBSD: Makefile,v 1.1 2006/01/19 23:23:58 elad Exp $
+
+NOMAN=yes
+PROG=   qsafe
+SRCS= qsafe.c qfile.c
+QSIEVE=${.CURDIR}/../qsieve
+.PATH: ${QSIEVE}
+CPPFLAGS+=-I${QSIEVE}
+
+.include <bsd.prog.mk>
diff --git a/moduli/qsafe/qsafe.c b/moduli/qsafe/qsafe.c
new file mode 100644
index 00000000..99c1e8c1
--- /dev/null
+++ b/moduli/qsafe/qsafe.c
@@ -0,0 +1,334 @@
+/* $NetBSD: qsafe.c,v 1.1 2006/01/19 23:23:58 elad Exp $ */
+
+/*-
+ * Copyright 1994 Phil Karn <karn@qualcomm.com>
+ * Copyright 1996-1998, 2003 William Allen206 Simpson <wsimpson@greendragon.com>
+ * Copyright 2000 Niels Provos <provos@citi.umich.edu>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/*
+ * Test probable "safe" primes,
+ *
+ *  suitable for use as Diffie-Hellman moduli;
+ *  that is, where q = (p-1)/2 is also prime.
+ *
+ * This is the second of two steps.
+ * This step is processor intensive.
+ *
+ * 1996 May     William Allen Simpson
+ *              extracted from earlier code by Phil Karn, April 1994.
+ *              read large prime candidates list (q),
+ *              and check prime probability of (p).
+ * 1998 May     William Allen Simpson
+ *              parameterized.
+ *              optionally limit to a single generator.
+ * 2000 Dec     Niels Provos
+ *              convert from GMP to openssl BN.
+ * 2003 Jun     William Allen Simpson
+ *              change outfile definition slightly to match openssh mistake.
+ *              move common file i/o to own file for better documentation.
+ *              redo debugprint again.
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <time.h>
+#include <openssl/bn.h>
+#include "qfile.h"
+
+/* define DEBUGPRINT     1 */
+#define TRIAL_MINIMUM           (4)
+
+static void     usage(void);
+
+/*
+ * perform a Miller-Rabin primality test
+ * on the list of candidates
+ * (checking both q and p)
+ * from standard input.
+ * The result is a list of so-call "safe" primes
+ * to standard output,
+ */
+int
+main(int argc, char *argv[])
+{
+	BIGNUM         *q, *p, *a;
+	BN_CTX         *ctx;
+	char           *cp;
+	char           *lp;
+	uint32_t        count_in = 0;
+	uint32_t        count_out = 0;
+	uint32_t        count_possible = 0;
+	uint32_t        generator_known;
+	uint32_t        generator_wanted = 0;
+	uint32_t        in_tests;
+	uint32_t        in_tries;
+	uint32_t        in_type;
+	uint32_t        in_size;
+	int             trials;
+	time_t          time_start;
+	time_t          time_stop;
+
+	setprogname(argv[0]);
+
+	if (argc < 2) {
+		usage();
+	}
+
+	if ((trials = strtoul(argv[1], NULL, 10)) < TRIAL_MINIMUM) {
+		trials = TRIAL_MINIMUM;
+	}
+
+	if (argc > 2) {
+		generator_wanted = strtoul(argv[2], NULL, 16);
+	}
+
+	time(&time_start);
+
+	p = BN_new();
+	q = BN_new();
+	ctx = BN_CTX_new();
+
+	(void)fprintf(stderr,
+		"%.24s Final %d Miller-Rabin trials (%x generator)\n",
+		ctime(&time_start), trials, generator_wanted);
+
+	lp = (char *) malloc((unsigned long) QLINESIZE + 1);
+
+	while (fgets(lp, QLINESIZE, stdin) != NULL) {
+		size_t          ll = strlen(lp);
+
+		count_in++;
+
+		if (ll < 14 || *lp == '!' || *lp == '#') {
+#ifdef  DEBUGPRINT
+			(void)fprintf(stderr, "%10lu: comment or short"
+				      " line\n", count_in);
+#endif
+			continue;
+		}
+
+		/* time */
+		cp = &lp[14];	/* (skip) */
+
+		/* type */
+		in_type = strtoul(cp, &cp, 10);
+
+		/* tests */
+		in_tests = strtoul(cp, &cp, 10);
+		if (in_tests & QTEST_COMPOSITE) {
+#ifdef  DEBUGPRINT
+			(void)fprintf(stderr, "%10lu: known composite\n",
+				count_in);
+#endif
+			continue;
+		}
+
+		/* tries */
+		in_tries = (uint32_t) strtoul(cp, &cp, 10);
+
+		/* size (most significant bit) */
+		in_size = (uint32_t) strtoul(cp, &cp, 10);
+
+		/* generator (hex) */
+		generator_known = (uint32_t) strtoul(cp, &cp, 16);
+
+		/* Skip white space */
+		cp += strspn(cp, " ");
+
+		/* modulus (hex) */
+		switch (in_type) {
+		case QTYPE_SOPHIE_GERMAINE:
+#ifdef  DEBUGPRINT
+			(void)fprintf(stderr, "%10lu: (%lu) "
+				      "Sophie-Germaine\n", count_in,
+				      in_type);
+#endif
+
+			a = q;
+			BN_hex2bn(&a, cp);
+			/* p = 2*q + 1 */
+			BN_lshift(p, q, 1);
+			BN_add_word(p, 1UL);
+			in_size += 1;
+			generator_known = 0;
+
+			break;
+
+		default:
+#ifdef  DEBUGPRINT
+			(void)fprintf(stderr, "%10lu: (%lu)\n",
+				      count_in,	in_type);
+#endif
+			a = p;
+			BN_hex2bn(&a, cp);
+			/* q = (p-1) / 2 */
+			BN_rshift(q, p, 1);
+
+			break;
+		}
+
+		/*
+		 * due to earlier inconsistencies in interpretation, check the
+		 * proposed bit size.
+		 */
+		if (BN_num_bits(p) != (in_size + 1)) {
+#ifdef  DEBUGPRINT
+			(void)fprintf(stderr, "%10lu: bit size %ul "
+				      "mismatch\n", count_in, in_size);
+#endif
+			continue;
+		}
+
+		if (in_size < QSIZE_MINIMUM) {
+#ifdef  DEBUGPRINT
+			(void)fprintf(stderr, "%10lu: bit size %ul "
+				      "too short\n", count_in, in_size);
+#endif
+			continue;
+		}
+
+		if (in_tests & QTEST_MILLER_RABIN)
+			in_tries += trials;
+		else
+			in_tries = trials;
+
+		/*
+		 * guess unknown generator
+		 */
+		if (generator_known == 0) {
+			if (BN_mod_word(p, 24UL) == 11)
+				generator_known = 2;
+			else if (BN_mod_word(p, 12UL) == 5)
+				generator_known = 3;
+			else {
+				BN_ULONG        r = BN_mod_word(p, 10UL);
+
+				if (r == 3 || r == 7) {
+					generator_known = 5;
+				}
+			}
+		}
+
+		/*
+		 * skip tests when desired generator doesn't match
+		 */
+		if (generator_wanted > 0 &&
+		    generator_wanted != generator_known) {
+#ifdef  DEBUGPRINT
+			(void)fprintf(stderr,
+				      "%10lu: generator %ld != %ld\n",
+				count_in, generator_known, generator_wanted);
+#endif
+			continue;
+		}
+
+		count_possible++;
+
+		/*
+		 * The (1/4)^N performance bound on Miller-Rabin is extremely
+		 * pessimistic, so don't spend a lot of time really verifying
+		 * that q is prime until after we know that p is also prime. A
+		 * single pass will weed out the vast majority of composite
+		 * q's.
+		 */
+		if (BN_is_prime(q, 1, NULL, ctx, NULL) <= 0) {
+#ifdef  DEBUGPRINT
+			(void)fprintf(stderr, "%10lu: q failed first "
+				      "possible prime test\n", count_in);
+#endif
+			continue;
+		}
+
+		/*
+		 * q is possibly prime, so go ahead and really make sure that
+		 * p is prime. If it is, then we can go back and do the same
+		 * for q. If p is composite, chances are that will show up on
+		 * the first Rabin-Miller iteration so it doesn't hurt to
+		 * specify a high iteration count.
+		 */
+		if (!BN_is_prime(p, trials, NULL, ctx, NULL)) {
+#ifdef  DEBUGPRINT
+			(void)fprintf(stderr, "%10lu: p is not prime\n",
+				      count_in);
+#endif
+			continue;
+		}
+
+#ifdef  DEBUGPRINT
+		(void)fprintf(stderr, "%10lu: p is almost certainly "
+			      "prime\n", count_in);
+#endif
+
+		/* recheck q more rigorously */
+		if (!BN_is_prime(q, trials - 1, NULL, ctx, NULL)) {
+#ifdef  DEBUGPRINT
+			(void)fprintf(stderr, "%10lu: q is not prime\n",
+				      count_in);
+#endif
+			continue;
+		}
+#ifdef  DEBUGPRINT
+		fprintf(stderr, "%10lu: q is almost certainly prime\n",
+			count_in);
+#endif
+		if (0 > qfileout(stdout,
+				 QTYPE_SAFE,
+				 (in_tests | QTEST_MILLER_RABIN),
+				 in_tries,
+				 in_size,
+				 generator_known,
+				 p)) {
+			break;
+		}
+		count_out++;
+#ifdef  DEBUGPRINT
+		fflush(stderr);
+		fflush(stdout);
+#endif
+	}
+
+	time(&time_stop);
+	free(lp);
+	BN_free(p);
+	BN_free(q);
+	BN_CTX_free(ctx);
+	fflush(stdout);		/* fclose(stdout); */
+	/* fclose(stdin); */
+	(void)fprintf(stderr,
+	   "%.24s Found %u safe primes of %u candidates in %lu seconds\n",
+		ctime(&time_stop), count_out, count_possible,
+		(long) (time_stop - time_start));
+
+	return (0);
+}
+
+static void
+usage(void)
+{
+	(void)fprintf(stderr, "Usage: %s <trials> [generator]\n",
+		      getprogname());
+	exit(1);
+}
diff --git a/moduli/qsieve/Makefile b/moduli/qsieve/Makefile
new file mode 100644
index 00000000..e6efe9da
--- /dev/null
+++ b/moduli/qsieve/Makefile
@@ -0,0 +1,8 @@
+# $NetBSD: Makefile,v 1.1 2006/01/19 23:23:58 elad Exp $
+
+PROG=   qsieve
+SRCS= qsieve.c qfile.c
+MAN= qsieve.6
+MLINKS = qsieve.6 qsafe.6
+
+.include <bsd.prog.mk>
diff --git a/moduli/qsieve/qfile.c b/moduli/qsieve/qfile.c
new file mode 100644
index 00000000..b6038794
--- /dev/null
+++ b/moduli/qsieve/qfile.c
@@ -0,0 +1,80 @@
+/* $NetBSD: qfile.c,v 1.1 2006/01/19 23:23:58 elad Exp $ */
+
+/*-
+ * Copyright 1994 Phil Karn <karn@qualcomm.com>
+ * Copyright 1996-1998, 2003 William Allen Simpson <wsimpson@greendragon.com>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/*
+ * File I/O for qsafe and qsieve.
+ * 
+ * 1996 May     William Allen Simpson extracted from earlier code by Phil Karn,
+ * April 1994. save large primes list for later processing. 1998 May
+ * William Allen Simpson parameterized. 2003 Jun     William Allen Simpson
+ * move common file i/o to own file for better documentation.
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <openssl/bn.h>
+#include "qfile.h"
+
+/*
+ * print moduli out in consistent form,
+ * normalizing error returns to printf-like expectations.
+ */
+int
+qfileout(FILE * ofile, uint32_t otype, uint32_t otests, uint32_t otries,
+	 uint32_t osize, uint32_t ogenerator, BIGNUM * omodulus)
+{
+	struct tm      *gtm;
+	time_t          time_now;
+
+	time(&time_now);
+	gtm = gmtime(&time_now);
+
+	if (0 > fprintf(ofile,
+			"%04d%02d%02d%02d%02d%02d %u %u %u %u %x ",
+			gtm->tm_year + 1900,
+			gtm->tm_mon + 1,
+			gtm->tm_mday,
+			gtm->tm_hour,
+			gtm->tm_min,
+			gtm->tm_sec,
+			otype,
+			otests,
+			otries,
+			osize,
+			ogenerator)) {
+
+		return (-1);
+	}
+
+	if (1 > BN_print_fp(ofile, omodulus)) {
+		return (-1);
+	}
+
+	return (fprintf(ofile, "\n"));
+}
diff --git a/moduli/qsieve/qfile.h b/moduli/qsieve/qfile.h
new file mode 100644
index 00000000..d24f6856
--- /dev/null
+++ b/moduli/qsieve/qfile.h
@@ -0,0 +1,71 @@
+/* $NetBSD: qfile.h,v 1.1 2006/01/19 23:23:58 elad Exp $ */
+
+/*-
+ * Copyright 1994 Phil Karn <karn@qualcomm.com>
+ * Copyright 1996-1998, 2003 William Allen Simpson <wsimpson@greendragon.com>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef _QFILE_H_
+#define _QFILE_H_
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <openssl/bn.h>
+
+/* need line long enough for largest moduli plus headers */
+#define QLINESIZE               (100+8192)
+
+/*
+ * Type: decimal. Specifies the internal structure of the prime modulus.
+ */
+#define QTYPE_UNKNOWN           (0)
+#define QTYPE_UNSTRUCTURED      (1)
+#define QTYPE_SAFE              (2)
+#define QTYPE_SCHNOOR           (3)
+#define QTYPE_SOPHIE_GERMAINE   (4)
+#define QTYPE_STRONG            (5)
+
+/*
+ * Tests: decimal (bit field). Specifies the methods used in checking for
+ * primality. Usually, more than one test is used.
+ */
+#define QTEST_UNTESTED          (0x00)
+#define QTEST_COMPOSITE         (0x01)
+#define QTEST_SIEVE             (0x02)
+#define QTEST_MILLER_RABIN      (0x04)
+#define QTEST_JACOBI            (0x08)
+#define QTEST_ELLIPTIC          (0x10)
+
+/*
+ * Size: decimal. Specifies the number of the most significant bit (0 to M). *
+ * WARNING: internally, usually 1 to N.
+ */
+#define QSIZE_MINIMUM           (511)
+
+
+int 
+qfileout(FILE *, uint32_t, uint32_t, uint32_t, uint32_t, uint32_t,
+	 BIGNUM *);
+
+#endif /* !_QFILE_H_ */
diff --git a/moduli/qsieve/qsieve.6 b/moduli/qsieve/qsieve.6
new file mode 100644
index 00000000..4279236d
--- /dev/null
+++ b/moduli/qsieve/qsieve.6
@@ -0,0 +1,117 @@
+.\"     $NetBSD: qsieve.6,v 1.1 2006/01/19 23:23:58 elad Exp $
+.\"
+.\" Copyright 1997, 2003 William Allen Simpson <wsimpson@greendragon.com>
+.\" All rights reserved.
+.\"
+.\" Redistribution and use in source and binary forms, with or without
+.\" modification, are permitted provided that the following conditions
+.\" are met:
+.\" 1. Redistributions of source code must retain the above copyright
+.\"    notice, this list of conditions and the following disclaimer.
+.\" 2. Redistributions in binary form must reproduce the above copyright
+.\"    notice, this list of conditions and the following disclaimer in the
+.\"    documentation and/or other materials provided with the distribution.
+.\" 3. All advertising materials mentioning features or use of this software
+.\"    must display the following acknowledgement:
+.\"      This product includes software designed by William Allen Simpson.
+.\" 4. The name of the author may not be used to endorse or promote products
+.\"    derived from this software without specific prior written permission.
+.\"
+.\" THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+.\" IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+.\" OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+.\" IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+.\" INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+.\" NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+.\" DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+.\" THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+.\" (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+.\" THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+.\"
+.\" Manual page, using -mandoc macros
+.\"
+.Dd July 28, 1997
+.Dt QSIEVE 1
+.Os
+.Sh NAME
+.Nm qsieve ,
+.Nm qsafe 
+.Nd generate system moduli file
+.Sh SYNOPSIS
+.Nm
+.Op Ar megabytes Ar bits Op Ar initial
+.br
+.Nm qsafe
+.Op Ar trials Op Ar generator
+.Sh DESCRIPTION
+The
+.Nm
+utility will list candidates for Sophie-Germaine primes
+(where q = (p-1)/2)
+to standard output.
+The list is checked against small known primes
+(less than 2**30).
+This step is both processor and memory intensive.
+.Pp
+The 
+.Ar megabytes
+value
+sets a limit for the internal sieve buffer.
+This should be small enough to remain entirely in memory.
+Swap thrashing can increase the run time
+from hours to days or weeks!
+When the 
+.Ar megabytes
+value is zero (0),
+.Nm
+will select a default suitable for the
+.Ar bits .
+.Pp
+The
+.Ar bits
+value
+sets the length of the generated possible primes
+(typically 768, 1024, 1536, 2048, 3072, or 4096,
+although others can be used for variety).
+.Pp
+The optional
+.Ar initial
+value (hex)
+specifies the beginning of the search.
+Otherwise,
+.Nm
+generates a randomly selected number.
+.Pp
+The
+.Nm qsafe
+utility will perform a Miller-Rabin primality test
+on the list of candidates
+(checking both q and p)
+from standard input.
+The result is a list of so-call "safe" primes
+to standard output,
+suitable for use as Diffie-Hellman moduli.
+This step is merely processor intensive.
+.Pp
+The
+.Ar trials
+value
+sets the number of Miller-Rabin interations
+(typically 16 to 128). 
+.Pp
+The optional
+.Ar generator
+value (hex)
+limits testing to candidates with a specific generator
+(usually 2). 
+Otherwise,
+.Nm qsafe
+will test each candidate
+and suggest a generator.
+.Sh SEE ALSO
+.Xr moduli 5
+.Sh HISTORY
+These programs were originally developed for
+the Photuris project,
+and later
+the OpenSSH project.
diff --git a/moduli/qsieve/qsieve.c b/moduli/qsieve/qsieve.c
new file mode 100644
index 00000000..34e5a339
--- /dev/null
+++ b/moduli/qsieve/qsieve.c
@@ -0,0 +1,473 @@
+/* $NetBSD: qsieve.c,v 1.1 2006/01/19 23:23:58 elad Exp $ */
+
+/*-
+ * Copyright 1994 Phil Karn <karn@qualcomm.com>
+ * Copyright 1996-1998, 2003 William Allen Simpson <wsimpson@greendragon.com>
+ * Copyright 2000 Niels Provos <provos@citi.umich.edu>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/*
+ * Sieve candidates for "safe" primes,
+ *  suitable for use as Diffie-Hellman moduli;
+ *  that is, where q = (p-1)/2 is also prime.
+ *
+ * This is the first of two steps.
+ * This step is memory intensive.
+ *
+ * 1996 May     William Allen Simpson
+ *              extracted from earlier code by Phil Karn, April 1994.
+ *              save large primes list for later processing.
+ * 1998 May     William Allen Simpson
+ *              parameterized.
+ * 2000 Dec     Niels Provos
+ *              convert from GMP to openssl BN.
+ * 2003 Jun     William Allen Simpson
+ *              change outfile definition slightly to match openssh mistake.
+ *              move common file i/o to own file for better documentation.
+ *              redo memory again.
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <openssl/bn.h>
+#include <string.h>
+#include <err.h>
+#include "qfile.h"
+
+/* define DEBUG_LARGE 1 */
+/* define DEBUG_SMALL 1 */
+
+/*
+ * Using virtual memory can cause thrashing.  This should be the largest
+ * number that is supported without a large amount of disk activity --
+ * that would increase the run time from hours to days or weeks!
+ */
+#define LARGE_MINIMUM   (8UL)	/* megabytes */
+
+/*
+ * Do not increase this number beyond the unsigned integer bit size.
+ * Due to a multiple of 4, it must be LESS than 128 (yielding 2**30 bits).
+ */
+#define LARGE_MAXIMUM   (127UL)	/* megabytes */
+
+/*
+ * Constant: assuming 8 bit bytes and 32 bit words
+ */
+#define SHIFT_BIT       (3)
+#define SHIFT_BYTE      (2)
+#define SHIFT_WORD      (SHIFT_BIT+SHIFT_BYTE)
+#define SHIFT_MEGABYTE  (20)
+#define SHIFT_MEGAWORD  (SHIFT_MEGABYTE-SHIFT_BYTE)
+
+/*
+ * Constant: when used with 32-bit integers, the largest sieve prime
+ * has to be less than 2**32.
+ */
+#define SMALL_MAXIMUM   (0xffffffffUL)
+
+/*
+ * Constant: can sieve all primes less than 2**32, as 65537**2 > 2**32-1.
+ */
+#define TINY_NUMBER     (1UL<<16)
+
+/*
+ * Ensure enough bit space for testing 2*q.
+ */
+#define TEST_MAXIMUM    (1UL<<16)
+#define TEST_MINIMUM    (QSIZE_MINIMUM + 1)
+/* real TEST_MINIMUM    (1UL << (SHIFT_WORD - TEST_POWER)) */
+#define TEST_POWER      (3)	/* 2**n, n < SHIFT_WORD */
+
+/*
+ * bit operations on 32-bit words
+ */
+#define BIT_CLEAR(a,n)  ((a)[(n)>>SHIFT_WORD] &= ~(1U << ((n) & 31)))
+#define BIT_SET(a,n)    ((a)[(n)>>SHIFT_WORD] |= (1U << ((n) & 31)))
+#define BIT_TEST(a,n)   ((a)[(n)>>SHIFT_WORD] & (1U << ((n) & 31)))
+
+/*
+ * sieve relative to the initial value
+ */
+uint32_t       *LargeSieve;
+uint32_t        largewords;
+uint32_t        largetries;
+uint32_t        largenumbers;
+uint32_t        largememory;	/* megabytes */
+uint32_t        largebits;
+BIGNUM         *largebase;
+
+/*
+ * sieve 2**30 in 2**16 parts
+ */
+uint32_t       *SmallSieve;
+uint32_t        smallbits;
+uint32_t        smallbase;
+
+/*
+ * sieve 2**16
+ */
+uint32_t       *TinySieve;
+uint32_t        tinybits;
+
+static void     usage(void);
+void            sieve_large(uint32_t);
+
+/*
+ * Sieve p's and q's with small factors
+ */
+void
+sieve_large(uint32_t s)
+{
+	BN_ULONG        r;
+	BN_ULONG        u;
+
+#ifdef  DEBUG_SMALL
+	(void)fprintf(stderr, "%lu\n", s);
+#endif
+	largetries++;
+	/* r = largebase mod s */
+	r = BN_mod_word(largebase, (BN_ULONG) s);
+	if (r == 0) {
+		/* s divides into largebase exactly */
+		u = 0;
+	} else {
+		/* largebase+u is first entry divisible by s */
+		u = s - r;
+	}
+
+	if (u < largebits * 2) {
+		/*
+		 * The sieve omits p's and q's divisible by 2, so ensure that
+		 * largebase+u is odd. Then, step through the sieve in
+		 * increments of 2*s
+		 */
+		if (u & 0x1) {
+			/* Make largebase+u odd, and u even */
+			u += s;
+		}
+
+		/* Mark all multiples of 2*s */
+		for (u /= 2; u < largebits; u += s) {
+			BIT_SET(LargeSieve, (uint32_t)u);
+		}
+	}
+
+	/* r = p mod s */
+	r = (2 * r + 1) % s;
+
+	if (r == 0) {
+		/* s divides p exactly */
+		u = 0;
+	} else {
+		/* p+u is first entry divisible by s */
+		u = s - r;
+	}
+
+	if (u < largebits * 4) {
+		/*
+		 * The sieve omits p's divisible by 4, so ensure that
+		 * largebase+u is not. Then, step through the sieve in
+		 * increments of 4*s
+		 */
+		while (u & 0x3) {
+			if (SMALL_MAXIMUM - u < s) {
+				return;
+			}
+
+			u += s;
+		}
+
+		/* Mark all multiples of 4*s */
+		for (u /= 4; u < largebits; u += s) {
+			BIT_SET(LargeSieve, (uint32_t)u);
+		}
+	}
+}
+
+/*
+ * list candidates for Sophie-Germaine primes
+ * (where q = (p-1)/2)
+ * to standard output.
+ * The list is checked against small known primes
+ * (less than 2**30).
+ */
+int
+main(int argc, char *argv[])
+{
+	BIGNUM         *q;
+	uint32_t        j;
+	int             power;
+	uint32_t        r;
+	uint32_t        s;
+	uint32_t        smallwords = TINY_NUMBER >> 6;
+	uint32_t        t;
+	time_t          time_start;
+	time_t          time_stop;
+	uint32_t        tinywords = TINY_NUMBER >> 6;
+	unsigned int    i;
+
+	setprogname(argv[0]);
+
+	if (argc < 3) {
+		usage();
+	}
+
+	/*
+         * Set power to the length in bits of the prime to be generated.
+         * This is changed to 1 less than the desired safe prime moduli p.
+         */
+	power = (int) strtoul(argv[2], NULL, 10);
+	if (power > TEST_MAXIMUM) {
+		errx(1, "Too many bits: %d > %lu.", power,
+		     (unsigned long)TEST_MAXIMUM);
+	} else if (power < TEST_MINIMUM) {
+		errx(1, "Too few bits: %d < %lu.", power,
+		     (unsigned long)TEST_MINIMUM);
+	}
+
+	power--;		/* decrement before squaring */
+
+	/*
+         * The density of ordinary primes is on the order of 1/bits, so the
+         * density of safe primes should be about (1/bits)**2. Set test range
+         * to something well above bits**2 to be reasonably sure (but not
+         * guaranteed) of catching at least one safe prime.
+	 */
+	largewords = (uint32_t)((unsigned long)
+			(power * power) >> (SHIFT_WORD - TEST_POWER));
+
+	/*
+         * Need idea of how much memory is available. We don't have to use all
+         * of it.
+	 */
+	largememory = (uint32_t)strtoul(argv[1], NULL, 10);
+	if (largememory > LARGE_MAXIMUM) {
+		warnx("Limited memory: %u MB; limit %lu MB.", largememory,
+		      LARGE_MAXIMUM);
+		largememory = LARGE_MAXIMUM;
+	}
+
+	if (largewords <= (largememory << SHIFT_MEGAWORD)) {
+		warnx("Increased memory: %u MB; need %u bytes.",
+		      largememory, (largewords << SHIFT_BYTE));
+		largewords = (largememory << SHIFT_MEGAWORD);
+	} else if (largememory > 0) {
+		warnx("Decreased memory: %u MB; want %u bytes.",
+		      largememory, (largewords << SHIFT_BYTE));
+		largewords = (largememory << SHIFT_MEGAWORD);
+	}
+
+	if ((TinySieve = (uint32_t *) calloc((size_t) tinywords, sizeof(uint32_t))) == NULL) {
+		errx(1, "Insufficient memory for tiny sieve: need %u byts.",
+		     tinywords << SHIFT_BYTE);
+	}
+	tinybits = tinywords << SHIFT_WORD;
+
+	if ((SmallSieve = (uint32_t *) calloc((size_t) smallwords, sizeof(uint32_t))) == NULL) {
+		errx(1, "Insufficient memory for small sieve: need %u bytes.",
+		     smallwords << SHIFT_BYTE);
+	}
+	smallbits = smallwords << SHIFT_WORD;
+
+	/*
+	 * dynamically determine available memory
+	 */
+	while ((LargeSieve = (uint32_t *)calloc((size_t)largewords,
+						sizeof(uint32_t))) == NULL) {
+		/* 1/4 MB chunks */
+		largewords -= (1L << (SHIFT_MEGAWORD - 2));
+	}
+	largebits = largewords << SHIFT_WORD;
+	largenumbers = largebits * 2;	/* even numbers excluded */
+
+	/* validation check: count the number of primes tried */
+	largetries = 0;
+
+	q = BN_new();
+	largebase = BN_new();
+
+	/*
+         * Generate random starting point for subprime search, or use
+         * specified parameter.
+	 */
+	if (argc < 4) {
+		BN_rand(largebase, power, 1, 1);
+	} else {
+		BIGNUM         *a;
+
+		a = largebase;
+		BN_hex2bn(&a, argv[2]);
+	}
+
+	/* ensure odd */
+	if (!BN_is_odd(largebase)) {
+		BN_set_bit(largebase, 0);
+	}
+
+	time(&time_start);
+	(void)fprintf(stderr,
+		"%.24s Sieve next %u plus %d-bit start point:\n# ",
+		ctime(&time_start), largenumbers, power);
+	BN_print_fp(stderr, largebase);
+	(void)fprintf(stderr, "\n");
+
+	/*
+         * TinySieve
+         */
+	for (i = 0; i < tinybits; i++) {
+		if (BIT_TEST(TinySieve, i)) {
+			/* 2*i+3 is composite */
+			continue;
+		}
+
+		/* The next tiny prime */
+		t = 2 * i + 3;
+
+		/* Mark all multiples of t */
+		for (j = i + t; j < tinybits; j += t) {
+			BIT_SET(TinySieve, j);
+		}
+
+		sieve_large(t);
+	}
+
+	/*
+         * Start the small block search at the next possible prime. To avoid
+         * fencepost errors, the last pass is skipped.
+         */
+	for (smallbase = TINY_NUMBER + 3;
+	     smallbase < (SMALL_MAXIMUM - TINY_NUMBER);
+	     smallbase += TINY_NUMBER) {
+		for (i = 0; i < tinybits; i++) {
+			if (BIT_TEST(TinySieve, i)) {
+				/* 2*i+3 is composite */
+				continue;
+			}
+
+			/* The next tiny prime */
+			t = 2 * i + 3;
+			r = smallbase % t;
+
+			if (r == 0) {
+				/* t divides into smallbase exactly */
+				s = 0;
+			} else {
+				/* smallbase+s is first entry divisible by t */
+				s = t - r;
+			}
+
+			/*
+			 * The sieve omits even numbers, so ensure that
+			 * smallbase+s is odd. Then, step through the sieve in
+			 * increments of 2*t
+			 */
+			if (s & 1) {
+				/* Make smallbase+s odd, and s even */
+				s += t;
+			}
+
+			/* Mark all multiples of 2*t */
+			for (s /= 2; s < smallbits; s += t) {
+				BIT_SET(SmallSieve, s);
+			}
+		}
+
+		/*
+                 * SmallSieve
+                 */
+		for (i = 0; i < smallbits; i++) {
+			if (BIT_TEST(SmallSieve, i)) {
+				/* 2*i+smallbase is composite */
+				continue;
+			}
+
+			/* The next small prime */
+			sieve_large((2 * i) + smallbase);
+		}
+
+		memset(SmallSieve, 0, (size_t)(smallwords << SHIFT_BYTE));
+	}
+
+	time(&time_stop);
+	(void)fprintf(stderr,
+		"%.24s Sieved with %u small primes in %lu seconds\n",
+		ctime(&time_stop), largetries,
+		(long) (time_stop - time_start));
+
+	for (j = r = 0; j < largebits; j++) {
+		if (BIT_TEST(LargeSieve, j)) {
+			/* Definitely composite, skip */
+			continue;
+		}
+
+#ifdef  DEBUG_LARGE
+		(void)fprintf(stderr, "test q = largebase+%lu\n", 2 * j);
+#endif
+
+		BN_set_word(q, (unsigned long)(2 * j));
+		BN_add(q, q, largebase);
+
+		if (0 > qfileout(stdout,
+				 (uint32_t) QTYPE_SOPHIE_GERMAINE,
+				 (uint32_t) QTEST_SIEVE,
+				 largetries,
+				 (uint32_t) (power - 1), /* MSB */
+				 (uint32_t) (0), /* generator unknown */
+				 q)) {
+			break;
+		}
+
+		r++;		/* count q */
+	}
+
+	time(&time_stop);
+
+	free(LargeSieve);
+	free(SmallSieve);
+	free(TinySieve);
+
+	fflush(stdout);
+	/* fclose(stdout); */
+
+	(void) fprintf(stderr, "%.24s Found %u candidates\n",
+	    ctime(&time_stop), r);
+
+	return (0);
+}
+
+static void
+usage(void)
+{
+	(void)fprintf(stderr, "Usage: %s <megabytes> <bits> [initial]\n"
+		"Possible values for <megabytes>: 0, %lu to %lu\n"
+		"Possible values for <bits>: %lu to %lu\n",
+		getprogname(),
+		LARGE_MINIMUM,
+		LARGE_MAXIMUM,
+		(unsigned long) TEST_MINIMUM,
+		(unsigned long) TEST_MAXIMUM);
+
+	exit(1);
+}