c_cpp/etc/calc/cal/randmprime.cal

   1 /*
   2  * randmprime - generate a random prime of the form h*2^n-1
   3  *
   4  * Copyright (C) 1999  Landon Curt Noll
   5  *
   6  * Calc is open software; you can redistribute it and/or modify it under
   7  * the terms of the version 2.1 of the GNU Lesser General Public License
   8  * as published by the Free Software Foundation.
   9  *
  10  * Calc is distributed in the hope that it will be useful, but WITHOUT
  11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  13  * Public License for more details.
  14  *
  15  * A copy of version 2.1 of the GNU Lesser General Public License is
  16  * distributed with calc under the filename COPYING-LGPL.  You should have
  17  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  18  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  19  *
  20  * @(#) $Revision: 30.1 $
  21  * @(#) $Id: randmprime.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
  22  * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/randmprime.cal,v $
  23  *
  24  * Under source code control:   1994/03/14 23:11:21
  25  * File existed as early as:    1994
  26  *
  27  * chongo <was here> /\oo/\     http://www.isthe.com/chongo/
  28  * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
  29  */
  30
  31
  32 /* obtain our required libs */
  33 read -once "lucas.cal"
  34
  35 /*
  36  * randmprime - find a random prime of the form h*2^n-1 of a given size
  37  *
  38  * given:
  39  *      bits    minimum bits in prime to return
  40  *      seed    random seed for srandom()
  41  *      [dbg]   if given, enable debugging
  42  *
  43  * returns:
  44  *      a prime of the form h*2^n-1
  45  */
  46 define
  47 randmprime(bits, seed, dbg)
  48 {
  49     local n;            /* n as in h*2^n-1 */
  50     local h;            /* h as in h*2^n-1 */
  51     local plush;        /* value added to h since the beginning */
  52     local init;         /* initial cpu time */
  53     local start;        /* cpu time before last test */
  54     local stop;         /* cpu time afte last test */
  55     local tmp;          /* just a tmp place holder value */
  56     local ret;          /* h*2^n-1 that is prime */
  57
  58     /* firewall */
  59     if (param(0) < 2 || param(0) > 3) {
  60         quit "bad usage: rndprime(dig, seed [,dbg])";
  61     }
  62     if (!isint(bits) || bits < 0 || !isint(seed) || seed < 0) {
  63         quit "args must be non-negative integers";
  64     }
  65     if (bits < 1) {
  66         bits = 1;
  67     }
  68     if (param(0) == 2 || dbg < 0) {
  69         dbg = 0;
  70     }
  71
  72     /* seed generator */
  73     tmp = srandom(seed, 13);
  74
  75     /* determine initial h and n values */
  76     n = random(bits>>1, highbit(bits)+bits>>1+1);
  77     h = randombit(n);
  78     h += iseven(h);
  79     while (highbit(h) >= n) {
  80        ++n;
  81     }
  82     if (dbg >= 1) {
  83         print "DEBUG3: initial h =", h;
  84         print "DEBUG3: initial n =", n;
  85     }
  86
  87     /*
  88      * loop until we find a prime
  89      */
  90     if (dbg >= 1) {
  91         start = usertime();
  92         init = usertime();
  93         plush = 0;
  94         print "DEBUG1: testing (h+" : plush : ")*2^" : n : "-1";
  95     }
  96     while (lucas(h,n) == 0) {
  97
  98         /* bump h, and n if needed */
  99         if (dbg >= 2) {
 100             stop = usertime();
 101             print "DEBUG2: last test:", stop-start, "   total time:", stop-init;
 102         }
 103         if (dbg >= 1) {
 104             print "DEBUG1: composite: (h+" : plush : ")*2^" : n : "-1";
 105             plush += 2;
 106         }
 107         h += 2;
 108         while (highbit(h) >= n) {
 109            ++n;
 110         }
 111         if (dbg >= 1) {
 112             print "DEBUG1: testing (h+" : plush : ")*2^" : n : "-1";
 113             start = stop;
 114         }
 115     }
 116
 117     /* found a prime */
 118     if (dbg >= 2) {
 119         stop = usertime();
 120         print "DEBUG2: last test:", stop-start, "   total time:", stop-init;
 121         print "DEBUG3: " : h : "*2^" : n : "-1 is prime";
 122     }
 123     if (dbg >= 1) {
 124         print "DEBUG1: prime: (h+" : plush : ")*2^" : n : "-1";
 125     }
 126     ret = h*2^n-1;
 127     if (dbg >= 3) {
 128         print "DEBUG3: highbit(h):", highbit(h);
 129         print "DEBUG3: digits(h):", digits(h);
 130         print "DEBUG3: highbit(n):", highbit(n);
 131         print "DEBUG3: digits(2^n):", int(n*ln(10)/ln(2)+1);
 132         print "DEBUG3: highbit(h*2^n-1):", highbit(ret);
 133         print "DEBUG3: digits(h*2^n)-1:", digits(ret);
 134     }
 135     return ret;
 136 }