mirror of
https://github.com/lcn2/calc.git
synced 2025-08-16 01:03:29 +03:00
135 lines
3.7 KiB
Plaintext
135 lines
3.7 KiB
Plaintext
/*
|
|
* randmprime - generate a random prime of the form h*2^n-1
|
|
*
|
|
* Copyright (c) 1997 by Landon Curt Noll. All Rights Reserved.
|
|
*
|
|
* Permission to use, copy, modify, and distribute this software and
|
|
* its documentation for any purpose and without fee is hereby granted,
|
|
* provided that the above copyright, this permission notice and text
|
|
* this comment, and the disclaimer below appear in all of the following:
|
|
*
|
|
* supporting documentation
|
|
* source copies
|
|
* source works derived from this source
|
|
* binaries derived from this source or from derived source
|
|
*
|
|
* LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
|
|
* INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO
|
|
* EVENT SHALL LANDON CURT NOLL BE LIABLE FOR ANY SPECIAL, INDIRECT OR
|
|
* CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF
|
|
* USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
|
|
* OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
|
|
* PERFORMANCE OF THIS SOFTWARE.
|
|
*
|
|
* Landon Curt Noll
|
|
* http://reality.sgi.com/chongo
|
|
*
|
|
* chongo <was here> /\../\
|
|
*/
|
|
|
|
/* obtain our required libs */
|
|
read -once "lucas.cal"
|
|
|
|
/*
|
|
* randmprime - find a random prime of the form h*2^n-1 of a given size
|
|
*
|
|
* given:
|
|
* bits minimum bits in prime to return
|
|
* seed random seed for srandom()
|
|
* [dbg] if given, enable debugging
|
|
*
|
|
* returns:
|
|
* a prime of the form h*2^n-1
|
|
*/
|
|
define
|
|
randmprime(bits, seed, dbg)
|
|
{
|
|
local n; /* n as in h*2^n-1 */
|
|
local h; /* h as in h*2^n-1 */
|
|
local plush; /* value added to h since the beginning */
|
|
local init; /* initial cpu time */
|
|
local start; /* cpu time before last test */
|
|
local stop; /* cpu time afte last test */
|
|
local tmp; /* just a tmp place holder value */
|
|
local ret; /* h*2^n-1 that is prime */
|
|
|
|
/* firewall */
|
|
if (param(0) < 2 || param(0) > 3) {
|
|
quit "bad usage: rndprime(dig, seed [,dbg])";
|
|
}
|
|
if (!isint(bits) || bits < 0 || !isint(seed) || seed < 0) {
|
|
quit "args must be non-negative integers";
|
|
}
|
|
if (bits < 1) {
|
|
bits = 1;
|
|
}
|
|
if (param(0) == 2 || dbg < 0) {
|
|
dbg = 0;
|
|
}
|
|
|
|
/* seed generator */
|
|
tmp = srandom(seed, 13);
|
|
|
|
/* determine initial h and n values */
|
|
n = random(bits>>1, highbit(bits)+bits>>1+1);
|
|
h = randombit(n);
|
|
h += iseven(h);
|
|
while (highbit(h) >= n) {
|
|
++n;
|
|
}
|
|
if (dbg >= 1) {
|
|
print "DEBUG3: initial h =", h;
|
|
print "DEBUG3: initial n =", n;
|
|
}
|
|
|
|
/*
|
|
* loop until we find a prime
|
|
*/
|
|
if (dbg >= 1) {
|
|
start = runtime();
|
|
init = runtime();
|
|
plush = 0;
|
|
print "DEBUG1: testing (h+" : plush : ")*2^" : n : "-1";
|
|
}
|
|
while (lucas(h,n) == 0) {
|
|
|
|
/* bump h, and n if needed */
|
|
if (dbg >= 2) {
|
|
stop = runtime();
|
|
print "DEBUG2: last test:", stop-start, " total time:", stop-init;
|
|
}
|
|
if (dbg >= 1) {
|
|
print "DEBUG1: composite: (h+" : plush : ")*2^" : n : "-1";
|
|
plush += 2;
|
|
}
|
|
h += 2;
|
|
while (highbit(h) >= n) {
|
|
++n;
|
|
}
|
|
if (dbg >= 1) {
|
|
print "DEBUG1: testing (h+" : plush : ")*2^" : n : "-1";
|
|
start = stop;
|
|
}
|
|
}
|
|
|
|
/* found a prime */
|
|
if (dbg >= 2) {
|
|
stop = runtime();
|
|
print "DEBUG2: last test:", stop-start, " total time:", stop-init;
|
|
print "DEBUG3: " : h : "*2^" : n : "-1 is prime";
|
|
}
|
|
if (dbg >= 1) {
|
|
print "DEBUG1: prime: (h+" : plush : ")*2^" : n : "-1";
|
|
}
|
|
ret = h*2^n-1;
|
|
if (dbg >= 3) {
|
|
print "DEBUG3: highbit(h):", highbit(h);
|
|
print "DEBUG3: digits(h):", digits(h);
|
|
print "DEBUG3: highbit(n):", highbit(n);
|
|
print "DEBUG3: digits(2^n):", int(n*ln(10)/ln(2)+1);
|
|
print "DEBUG3: highbit(h*2^n-1):", highbit(ret);
|
|
print "DEBUG3: digits(h*2^n)-1:", digits(ret);
|
|
}
|
|
return ret;
|
|
}
|