mirror of
https://github.com/lcn2/calc.git
synced 2025-08-16 01:03:29 +03:00
Release calc version 2.11.1
This commit is contained in:
Landon Curt Noll
319
cal/mfactor.cal
Normal file
319
cal/mfactor.cal
Normal file
@@ -0,0 +1,319 @@
|
||||
/*
|
||||
* mfactor - return the lowest factor of 2^n-1, for n > 0
|
||||
*
|
||||
* Copyright (C) 1999 Landon Curt Noll
|
||||
*
|
||||
* Calc is open software; you can redistribute it and/or modify it under
|
||||
* the terms of the version 2.1 of the GNU Lesser General Public License
|
||||
* as published by the Free Software Foundation.
|
||||
*
|
||||
* Calc is distributed in the hope that it will be useful, but WITHOUT
|
||||
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
||||
* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
|
||||
* Public License for more details.
|
||||
*
|
||||
* A copy of version 2.1 of the GNU Lesser General Public License is
|
||||
* distributed with calc under the filename COPYING-LGPL. You should have
|
||||
* received a copy with calc; if not, write to Free Software Foundation, Inc.
|
||||
* 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
|
||||
*
|
||||
* @(#) $Revision: 29.1 $
|
||||
* @(#) $Id: mfactor.cal,v 29.1 1999/12/14 09:15:31 chongo Exp $
|
||||
* @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/mfactor.cal,v $
|
||||
*
|
||||
* Under source code control: 1996/07/06 06:09:40
|
||||
* File existed as early as: 1996
|
||||
*
|
||||
* chongo <was here> /\oo/\ http://reality.sgi.com/chongo/
|
||||
* Share and enjoy! :-) http://reality.sgi.com/chongo/tech/comp/calc/
|
||||
*/
|
||||
|
||||
/*
|
||||
* hset method
|
||||
*
|
||||
* We will assume that mfactor is called with p_elim == 17.
|
||||
*
|
||||
* n = (the Mersenne exponent we are testing)
|
||||
* Q = 4*2*3*5*7*11*13*17 (4 * pfact(of some reasonable integer))
|
||||
*
|
||||
* We first determine all values of h mod Q such that:
|
||||
*
|
||||
* gcd(h*n+1, Q) == 1 and h*n+1 == +/-1 mod 8
|
||||
*
|
||||
* There will be 2*1*2*4*6*10*12*16 such values of h.
|
||||
*
|
||||
* For efficiency, we keep the difference between consecutive h values
|
||||
* in the hset[] difference array with hset[0] being the first h value.
|
||||
* Last, we multiply the hset[] values by n so that we only need
|
||||
* to add sequential values of hset[] to get factor candidates.
|
||||
*
|
||||
* We need only test factors of the form:
|
||||
*
|
||||
* (Q*g*n + hx) + 1
|
||||
*
|
||||
* where:
|
||||
*
|
||||
* g is an integer >= 0
|
||||
* hx is computed from hset[] difference value described above
|
||||
*
|
||||
* Note that (Q*g*n + hx) is always even and that hx is a multiple
|
||||
* of n. Thus the typical factor form:
|
||||
*
|
||||
* 2*k*n + 1
|
||||
*
|
||||
* implies that:
|
||||
*
|
||||
* k = (Q*g + hx/n)/2
|
||||
*
|
||||
* This allows us to quickly eliminate factor values that are divisible
|
||||
* by 2, 3, 5, 7, 11, 13 or 17. (well <= p value found below)
|
||||
*
|
||||
* The following loop shows how test_factor is advanced to higher test
|
||||
* values using hset[]. Here, hcount is the number of elements in hset[].
|
||||
* It can be shown that hset[0] == 0. We add hset[hcount] to the hset[]
|
||||
* array for looping control convenience.
|
||||
*
|
||||
* (* increase test_factor thru other possible test values *)
|
||||
* test_factor = 0;
|
||||
* hindx = 0;
|
||||
* do {
|
||||
* while (++hindx <= hcount) {
|
||||
* test_factor += hset[hindx];
|
||||
* }
|
||||
* hindx = 0;
|
||||
* } while (test_factor < some_limit);
|
||||
*
|
||||
* The test, mfactor(67, 1, 10000) took on an 200 Mhz r4k (user CPU seconds):
|
||||
*
|
||||
* 210.83 (prior to use of hset[])
|
||||
* 78.35 (hset[] for p_elim = 7)
|
||||
* 73.87 (hset[] for p_elim = 11)
|
||||
* 73.92 (hset[] for p_elim = 13)
|
||||
* 234.16 (hset[] for p_elim = 17)
|
||||
* p_elim == 19 requires over 190 Megs of memory
|
||||
*
|
||||
* Over a long period of time, the call to load_hset() becomes insignificant.
|
||||
* If we look at the user CPU seconds from the first 10000 cycle to the
|
||||
* end of the test we find:
|
||||
*
|
||||
* 205.00 (prior to use of hset[])
|
||||
* 75.89 (hset[] for p_elim = 7)
|
||||
* 73.74 (hset[] for p_elim = 11)
|
||||
* 70.61 (hset[] for p_elim = 13)
|
||||
* 57.78 (hset[] for p_elim = 17)
|
||||
* p_elim == 19 rejected because of memory size
|
||||
*
|
||||
* The p_elim == 17 overhead takes ~3 minutes on an 200 Mhz r4k CPU and
|
||||
* requires about ~13 Megs of memory. The p_elim == 13 overhead
|
||||
* takes about 3 seconds and requires ~1.5 Megs of memory.
|
||||
*
|
||||
* The value p_elim == 17 is best for long factorizations. It is the
|
||||
* fastest even thought the initial startup overhead is larger than
|
||||
* for p_elim == 13.
|
||||
*
|
||||
* NOTE: The values above are prior to optimizations where hset[] was
|
||||
* multiplied by n plus other optimizations. Thus, the CPU
|
||||
* times you may get will not likely match the above values.
|
||||
*/
|
||||
|
||||
|
||||
/*
|
||||
* mfactor - find a factor of a Mersenne Number
|
||||
*
|
||||
* Mersenne numbers are numbers of the form:
|
||||
*
|
||||
* 2^n-1 for integer n > 0
|
||||
*
|
||||
* We know that factors of a Mersenne number are of the form:
|
||||
*
|
||||
* 2*k*n+1 and +/- 1 mod 8
|
||||
*
|
||||
* We make use of the hset[] difference array to eliminate factor
|
||||
* candidates that would otherwise be divisible by 2, 3, 5, 7 ... p_elim.
|
||||
*
|
||||
* given:
|
||||
* n attempt to factor M(n) = 2^n-1
|
||||
* start_k the value k in 2*k*n+1 to start the search (def: 1)
|
||||
* rept_loop loop cycle reporting (def: 10000)
|
||||
* p_elim largest prime to eliminate from test factors (def: 17)
|
||||
*
|
||||
* returns:
|
||||
* factor of (2^n)-1
|
||||
*
|
||||
* NOTE: The p_elim argument is optional and defaults to 17. A p_elim value
|
||||
* of 17 is faster than 13 for even medium length runs. However 13
|
||||
* uses less memory and has a shorter startup time.
|
||||
*/
|
||||
define mfactor(n, start_k, rept_loop, p_elim)
|
||||
{
|
||||
local Q; /* 4*pfact(p_elim), hset[] cycle size */
|
||||
local hcount; /* elements in the hset[] difference array */
|
||||
local loop; /* report loop count */
|
||||
local q; /* test factor of 2^n-1 */
|
||||
local g; /* g as in test candidate form: (Q*g*hset[i])*n + 1 */
|
||||
local hindx; /* hset[] index */
|
||||
local i;
|
||||
local tmp;
|
||||
local tmp2;
|
||||
|
||||
/*
|
||||
* firewall
|
||||
*/
|
||||
if (!isint(n) || n <= 0) {
|
||||
quit "n must be an integer > 0";
|
||||
}
|
||||
if (!isint(start_k)) {
|
||||
start_k = 1;
|
||||
} else if (!isint(start_k) || start_k <= 0) {
|
||||
quit "start_k must be an integer > 0";
|
||||
}
|
||||
if (!isint(rept_loop)) {
|
||||
rept_loop = 10000;
|
||||
}
|
||||
if (rept_loop < 1) {
|
||||
quit "rept_loop must be an integer > 0";
|
||||
}
|
||||
if (!isint(p_elim)) {
|
||||
p_elim = 17;
|
||||
}
|
||||
if (p_elim < 3) {
|
||||
quit "p_elim must be an integer > 2 (try 13 or 17)";
|
||||
}
|
||||
|
||||
/*
|
||||
* declare our global values
|
||||
*/
|
||||
Q = 4*pfact(p_elim);
|
||||
hcount = 2;
|
||||
/* allocate the h difference array */
|
||||
for (i=2; i <= p_elim; i = nextcand(i)) {
|
||||
hcount *= (i-1);
|
||||
}
|
||||
local mat hset[hcount+1];
|
||||
|
||||
/*
|
||||
* load the hset[] difference array
|
||||
*/
|
||||
{
|
||||
local x; /* h*n+1 mod 8 */
|
||||
local h; /* potential h value */
|
||||
local last_h; /* previous valid h value */
|
||||
|
||||
last_h = 0;
|
||||
for (i=0,h=0; h < Q; ++h) {
|
||||
if (gcd(h*n+1,Q) == 1) {
|
||||
x = (h*n+1) % 8;
|
||||
if (x == 1 || x == 7) {
|
||||
hset[i++] = (h-last_h) * n;
|
||||
last_h = h;
|
||||
}
|
||||
}
|
||||
}
|
||||
hset[hcount] = Q*n - last_h*n;
|
||||
}
|
||||
|
||||
/*
|
||||
* setup
|
||||
*
|
||||
* determine the next g and hset[] index (hindx) values such that:
|
||||
*
|
||||
* 2*start_k <= (Q*g + hset[hindx])
|
||||
*
|
||||
* and (Q*g + hset[hindx]) is a minimum and where:
|
||||
*
|
||||
* Q = (4 * pfact(of some reasonable integer))
|
||||
* g = (some integer) (hset[] cycle number)
|
||||
*
|
||||
* We also compute 'q', the next test candidate.
|
||||
*/
|
||||
g = (2*start_k) // Q;
|
||||
tmp = 2*start_k - Q*g;
|
||||
for (tmp2=0, hindx=0;
|
||||
hindx < hcount && (tmp2 += hset[hindx]/n) < tmp;
|
||||
++hindx) {
|
||||
}
|
||||
if (hindx == hcount) {
|
||||
/* we are beyond the end of a hset[] cycle, start at the next */
|
||||
++g;
|
||||
hindx = 0;
|
||||
tmp2 = hset[0]/n;
|
||||
}
|
||||
q = (Q*g + tmp2)*n + 1;
|
||||
|
||||
/*
|
||||
* look for a factor
|
||||
*
|
||||
* We ignore factors that themselves are divisible by a prime <=
|
||||
* some small prime p.
|
||||
*
|
||||
* This process is guaranteed to find the smallest factor
|
||||
* of 2^n-1. A smallest factor of 2^n-1 must be prime, otherwise
|
||||
* the divisors of that factor would also be factors of 2^n-1.
|
||||
* Thus we know that if a test factor itself is not prime, it
|
||||
* cannot be the smallest factor of 2^n-1.
|
||||
*
|
||||
* Eliminating all non-prime test factors would take too long.
|
||||
* However we can eliminate 80.81% of the test factors
|
||||
* by not using test factors that are divisible by a prime <= 17.
|
||||
*/
|
||||
if (pmod(2,n,q) == 1) {
|
||||
return q;
|
||||
} else {
|
||||
/* report this loop */
|
||||
printf("at 2*%d*%d+1, cpu: %f\n",
|
||||
(q-1)/(2*n), n, runtime());
|
||||
fflush(files(1));
|
||||
loop = 0;
|
||||
}
|
||||
do {
|
||||
/*
|
||||
* determine if we need to report
|
||||
*
|
||||
* NOTE: (q-1)/(2*n) is the k value from 2*k*n + 1.
|
||||
*/
|
||||
if (rept_loop <= ++loop) {
|
||||
/* report this loop */
|
||||
printf("at 2*%d*%d+1, cpu: %f\n",
|
||||
(q-1)/(2*n), n, runtime());
|
||||
fflush(files(1));
|
||||
loop = 0;
|
||||
}
|
||||
|
||||
/*
|
||||
* skip if divisable by a prime <= 449
|
||||
*
|
||||
* The value 281 was determined by timing loops
|
||||
* which found that 281 was at or near the
|
||||
* minimum time to factor 2^(2^127-1)-1.
|
||||
*
|
||||
* The addition of the do { ... } while (factor(q, 449)>1);
|
||||
* loop reduced the factoring loop time (36504 k values with
|
||||
* the hset[] initialization time removed) from 25.69 sec to
|
||||
* 15.62 sec of CPU time on a 200Mhz r4k.
|
||||
*/
|
||||
do {
|
||||
/*
|
||||
* determine the next factor candidate
|
||||
*/
|
||||
q += hset[++hindx];
|
||||
if (hindx >= hcount) {
|
||||
hindx = 0;
|
||||
/*
|
||||
* if we cared about g,
|
||||
* then we wound ++g here too
|
||||
*/
|
||||
}
|
||||
} while (factor(q, 449) > 1);
|
||||
} while (pmod(2,n,q) != 1);
|
||||
|
||||
/*
|
||||
* return the factor found
|
||||
*
|
||||
* q is a factor of (2^n)-1
|
||||
*/
|
||||
return q;
|
||||
}
|
||||
|
||||
if (config("resource_debug") & 3) {
|
||||
print "mfactor(n [, start_k=1 [, rept_loop=10000 [, p_elim=17]]])"
|
||||
}
|
||||
Reference in New Issue
Block a user