germtest.c 8.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336
  1. /*
  2. * Copyright (c) 1995 Colin Plumb. All rights reserved.
  3. * For licensing and other legal details, see the file legal.c.
  4. *
  5. * germtest.c - Random Sophie Germain prime generator.
  6. *
  7. * This generates random Sophie Germain primes using the command line
  8. * as a seed value. It uses George Marsaglia's "mother of all random
  9. * number generators" to (using the command line as a seed) to pick the
  10. * starting search value and then searches sequentially for the next
  11. * Sophie Germain prime p (a prime such that 2*p+1 is also prime).
  12. *
  13. * This is a really good way to burn a lot of CPU cycles.
  14. */
  15. #if HAVE_CONFIG_H
  16. #include "bnconfig.h"
  17. #endif
  18. #include <stdio.h>
  19. #if !NO_STRING_H
  20. #include <string.h>
  21. #elif HAVE_STRINGS_H
  22. #include <strings.h>
  23. #endif
  24. #include <stdlib.h> /* For malloc() */
  25. #include "bn.h"
  26. #include "germain.h"
  27. #include "sieve.h"
  28. #include "cputime.h"
  29. #define BNDEBUG 1
  30. #include "bnprint.h"
  31. #define bnPut(prompt, bn) bnPrint(stdout, prompt, bn, "\n")
  32. /*
  33. * Generate random numbers according to George Marsaglia's
  34. * Mother Of All Random Number Generators. This has a
  35. * period of 0x17768215025F82EA0378038A03A203CA7FFF,
  36. * or decimal 2043908804452974490458343567652678881935359.
  37. */
  38. static unsigned mstate[8];
  39. static unsigned mcarry;
  40. static unsigned mindex;
  41. static unsigned
  42. mRandom_16(void)
  43. {
  44. unsigned long t;
  45. t = mcarry +
  46. mstate[ mindex ] * 1941ul +
  47. mstate[(mindex+1)&7] * 1860ul +
  48. mstate[(mindex+2)&7] * 1812ul +
  49. mstate[(mindex+3)&7] * 1776ul +
  50. mstate[(mindex+4)&7] * 1492ul +
  51. mstate[(mindex+5)&7] * 1215ul +
  52. mstate[(mindex+6)&7] * 1066ul +
  53. mstate[(mindex+7)&7] * 12013ul;
  54. mcarry = (unsigned)(t >> 16); /* 0 <= mcarry <= 0x5a87 */
  55. mindex = (mindex-1) & 7;
  56. return mstate[mindex] = (unsigned)(t & 0xffff);
  57. }
  58. /*
  59. * Initialize the RNG based on the given seed.
  60. * A zero-length seed will produce pretty lousy numbers,
  61. * but it will work.
  62. */
  63. static void
  64. mSeed(unsigned char const *seed, unsigned len)
  65. {
  66. unsigned i;
  67. for (i = 0; i < 8; i++)
  68. mstate[i] = 0;
  69. mcarry = 1;
  70. while (len--) {
  71. mcarry += *seed++;
  72. (void)mRandom_16();
  73. }
  74. }
  75. /*
  76. * Generate a bignum of a specified length, with the given
  77. * high and low 8 bits. "High" is merged into the high 8 bits of the
  78. * number. For example, set it to 0x80 to ensure that the number is
  79. * exactly "bits" bits long (i.e. 2^(bits-1) <= bn < 2^bits).
  80. * "Low" is merged into the low 8 bits. For example, set it to
  81. * 1 to ensure that you generate an odd number. "High" is merged
  82. * into the high bits; set it to 0x80 to ensure that the high bit
  83. * is set in the returned value.
  84. */
  85. static int
  86. genRandBn(struct BigNum *bn, unsigned bits, unsigned char high,
  87. unsigned char low, unsigned char const *seed, unsigned len)
  88. {
  89. unsigned char buf[64];
  90. unsigned bytes;
  91. unsigned l = 0; /* Current position */
  92. unsigned t, i;
  93. bnSetQ(bn, 0);
  94. if (bnPrealloc(bn, bits) < 0)
  95. return -1;
  96. mSeed(seed, len);
  97. bytes = (bits+7) / 8; /* Number of bytes to use */
  98. for (i = 0; i < sizeof(buf); i += 2) {
  99. t = mRandom_16();
  100. buf[i] = (unsigned char)(t >> 8);
  101. buf[i+1] = (unsigned char)t;
  102. }
  103. buf[sizeof(buf)-1] |= low;
  104. while (bytes > sizeof(buf)) {
  105. bytes -= sizeof(buf);
  106. /* Merge in low half of high bits, if necessary */
  107. if (bytes == 1 && (bits & 7))
  108. buf[0] |= high << (bits & 7);
  109. if (bnInsertBigBytes(bn, buf, l, sizeof(buf)) < 0)
  110. return -1;
  111. l += sizeof(buf);
  112. for (i = 0; i < sizeof(buf); i += 2) {
  113. t = mRandom_16();
  114. buf[i] = (unsigned char)t;
  115. buf[i+1] = (unsigned char)(t >> 8);
  116. }
  117. }
  118. /* Do the final "bytes"-long section, using the tail bytes in buf */
  119. /* Mask off excess high bits */
  120. buf[sizeof(buf)-bytes] &= 255 >> (-bits & 7);
  121. /* Merge in specified high bits */
  122. buf[sizeof(buf)-bytes] |= high >> (-bits & 7);
  123. if (bytes > 1 && (bits & 7))
  124. buf[sizeof(buf)-bytes+1] |= high << (bits & 7);
  125. /* Merge in the appropriate bytes of the buffer */
  126. if (bnInsertBigBytes(bn, buf+sizeof(buf)-bytes, l, bytes) < 0)
  127. return -1;
  128. return 0;
  129. }
  130. struct Progress {
  131. FILE *f;
  132. unsigned column;
  133. unsigned wrap;
  134. };
  135. /* Print a progress indicator, with line-wrap */
  136. static int
  137. genProgress(void *arg, int c)
  138. {
  139. struct Progress *p = arg;
  140. if (++p->column > p->wrap) {
  141. putc('\n', p->f);
  142. p->column = 1;
  143. }
  144. putc(c, p->f);
  145. fflush(p->f);
  146. return 0;
  147. }
  148. static int
  149. genSophieGermain(struct BigNum *bn, unsigned bits, unsigned order,
  150. unsigned char const *seed, unsigned len, FILE *f)
  151. {
  152. #if CLOCK_AVAIL
  153. timetype start, stop;
  154. unsigned long s;
  155. #endif
  156. int i;
  157. #if BNDEBUG
  158. unsigned char s1[1024], s2[1024];
  159. #endif
  160. char buf[40];
  161. unsigned p1, p2;
  162. struct BigNum step;
  163. struct Progress progress;
  164. if (f)
  165. fprintf(f, "Generating a %u-bit order-%u Sophie Germain prime with \"%.*s\"\n",
  166. bits, order, (int)len, (char *)seed);
  167. progress.f = f;
  168. progress.column = 0;
  169. progress.wrap = 78;
  170. /* Find p - choose a starting place */
  171. if (genRandBn(bn, bits, 0xC0, 3, seed, len) < 0)
  172. return -1;
  173. #if BNDEBUG /* DEBUG - check that sieve works properly */
  174. bnBegin(&step);
  175. bnSetQ(&step, 2);
  176. sieveBuild(s1, 1024, bn, 2, order);
  177. sieveBuildBig(s2, 1024, bn, &step, order);
  178. p1 = p2 = 0;
  179. if (s1[0] != s2[0])
  180. printf("Difference: s1[0] = %x s2[0] = %x\n", s1[0], s2[0]);
  181. do {
  182. p1 = sieveSearch(s1, 1024, p1);
  183. p2 = sieveSearch(s2, 1024, p2);
  184. if (p1 != p2)
  185. printf("Difference: p1 = %u p2 = %u\n", p1, p2);
  186. } while (p1 && p2);
  187. bnEnd(&step);
  188. #endif
  189. /* And search for a prime */
  190. #if CLOCK_AVAIL
  191. gettime(&start);
  192. #endif
  193. i = germainPrimeGen(bn, order, f ? genProgress : 0, (void *)&progress);
  194. if (i < 0)
  195. return -1;
  196. #if CLOCK_AVAIL
  197. gettime(&stop);
  198. #endif
  199. if (f) {
  200. putc('\n', f);
  201. fprintf(f, "%d modular exponentiations performed.\n", i);
  202. }
  203. #if CLOCK_AVAIL
  204. subtime(stop, start);
  205. s = sec(stop);
  206. printf("%u-bit time = %lu.%03u sec.", bits, s, msec(stop));
  207. if (s > 60) {
  208. putchar(' ');
  209. putchar('(');
  210. if (s > 3600)
  211. printf("%u:%02u", (unsigned)(s/3600),
  212. (unsigned)(s/60%60));
  213. else
  214. printf("%u", (unsigned)(s/60));
  215. printf(":%02u)", (unsigned)(s%60));
  216. }
  217. putchar('\n');
  218. #endif
  219. bnPut(" p = ", bn);
  220. for (p1 = 0; p1 < order; p1++) {
  221. if (bnLShift(bn, 1) <0)
  222. return -1;
  223. (void)bnAddQ(bn, 1);
  224. sprintf(buf, "%u*p+%u = ", 2u<<p1, (2u<<p1) - 1);
  225. bnPut(buf, bn);
  226. }
  227. return 0;
  228. }
  229. /* Copy the command line to the buffer. */
  230. static unsigned char *
  231. copy(int argc, char **argv, size_t *lenp)
  232. {
  233. size_t len;
  234. int i;
  235. unsigned char *buf, *p;
  236. len = argc > 2 ? (size_t)(argc-2) : 0;
  237. for (i = 1; i < argc; i++)
  238. len += strlen(argv[i]);
  239. *lenp = len;
  240. buf = malloc(len+!len); /* Can't malloc 0 bytes... */
  241. if (buf) {
  242. p = buf;
  243. for (i = 1; i < argc; i++) {
  244. if (i > 1)
  245. *p++ = ' ';
  246. len = strlen(argv[i]);
  247. memcpy(p, argv[i], len);
  248. p += len;
  249. }
  250. }
  251. return buf;
  252. }
  253. int
  254. main(int argc, char **argv)
  255. {
  256. unsigned len;
  257. struct BigNum bn;
  258. unsigned char *buf;
  259. if (argc < 2) {
  260. fprintf(stderr, "Usage: %s <seed>\n", argv[0]);
  261. fputs("\
  262. <seed> should be a a string of bytes to be hashed to seed the prime\n\
  263. generator. Note that unquoted whitespace between words will be counted\n\
  264. as a single space. To include multiple spaces, quote them.\n", stderr);
  265. return 1;
  266. }
  267. buf = copy(argc, argv, &len);
  268. if (!buf) {
  269. fprintf(stderr, "Out of memory!\n");
  270. return 1;
  271. }
  272. bnBegin(&bn);
  273. genSophieGermain(&bn, 0x100, 0, buf, len, stdout);
  274. genSophieGermain(&bn, 0x100, 1, buf, len, stdout);
  275. genSophieGermain(&bn, 0x100, 2, buf, len, stdout);
  276. genSophieGermain(&bn, 0x100, 3, buf, len, stdout);
  277. genSophieGermain(&bn, 0x200, 0, buf, len, stdout);
  278. genSophieGermain(&bn, 0x200, 1, buf, len, stdout);
  279. genSophieGermain(&bn, 0x200, 2, buf, len, stdout);
  280. genSophieGermain(&bn, 0x300, 0, buf, len, stdout);
  281. genSophieGermain(&bn, 0x300, 1, buf, len, stdout);
  282. genSophieGermain(&bn, 0x400, 0, buf, len, stdout);
  283. genSophieGermain(&bn, 0x400, 1, buf, len, stdout);
  284. genSophieGermain(&bn, 0x500, 0, buf, len, stdout);
  285. genSophieGermain(&bn, 0x500, 1, buf, len, stdout);
  286. genSophieGermain(&bn, 0x600, 0, buf, len, stdout);
  287. genSophieGermain(&bn, 0x600, 1, buf, len, stdout);
  288. #if 0
  289. /* These get *really* slow */
  290. genSophieGermain(&bn, 0x800, 0, buf, len, stdout);
  291. genSophieGermain(&bn, 0x800, 1, buf, len, stdout);
  292. genSophieGermain(&bn, 0xc00, 0, buf, len, stdout);
  293. genSophieGermain(&bn, 0xc00, 1, buf, len, stdout);
  294. /* Like, plan on a *week* or more for this one. */
  295. genSophieGermain(&bn, 0x1000, 0, buf, len, stdout);
  296. genSophieGermain(&bn, 0x1000, 1, buf, len, stdout);
  297. #endif
  298. bnEnd(&bn);
  299. free(buf);
  300. return 0;
  301. }