1/*
2 * mfactor - return the lowest factor of 2^n-1, for n > 0
3 *
4 * Copyright (C) 1999,2021  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 * Under source code control:	1996/07/06 06:09:40
21 * File existed as early as:	1996
22 *
23 * chongo <was here> /\oo/\	http://www.isthe.com/chongo/
24 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
25 */
26
27/*
28 * hset method
29 *
30 * We will assume that mfactor is called with p_elim == 17.
31 *
32 *	n = (the Mersenne exponent we are testing)
33 *	Q = 4*2*3*5*7*11*13*17 (4 * pfact(of some reasonable integer))
34 *
35 * We first determine all values of h mod Q such that:
36 *
37 *	gcd(h*n+1, Q) == 1   and   h*n+1 == +/-1 mod 8
38 *
39 * There will be 2*1*2*4*6*10*12*16 such values of h.
40 *
41 * For efficiency, we keep the difference between consecutive h values
42 * in the hset[] difference array with hset[0] being the first h value.
43 * Last, we multiply the hset[] values by n so that we only need
44 * to add sequential values of hset[] to get factor candidates.
45 *
46 * We need only test factors of the form:
47 *
48 *	(Q*g*n + hx) + 1
49 *
50 * where:
51 *
52 *	g is an integer >= 0
53 *	hx is computed from hset[] difference value described above
54 *
55 * Note that (Q*g*n + hx) is always even and that hx is a multiple
56 * of n.  Thus the typical factor form:
57 *
58 *	2*k*n + 1
59 *
60 * implies that:
61 *
62 *	k = (Q*g + hx/n)/2
63 *
64 * This allows us to quickly eliminate factor values that are divisible
65 * by 2, 3, 5, 7, 11, 13 or 17.	 (well <= p value found below)
66 *
67 * The following loop shows how test_factor is advanced to higher test
68 * values using hset[].	 Here, hcount is the number of elements in hset[].
69 * It can be shown that hset[0] == 0.  We add hset[hcount] to the hset[]
70 * array for looping control convenience.
71 *
72 *	(* increase test_factor thru other possible test values *)
73 *	test_factor = 0;
74 *	hindx = 0;
75 *	do {
76 *		while (++hindx <= hcount) {
77 *			test_factor += hset[hindx];
78 *		}
79 *		hindx = 0;
80 *	} while (test_factor < some_limit);
81 *
82 * The test, mfactor(67, 1, 10000) took on an 200 MHz r4k (user CPU seconds):
83 *
84 *	210.83	(prior to use of hset[])
85 *	 78.35	(hset[] for p_elim = 7)
86 *	 73.87	(hset[] for p_elim = 11)
87 *	 73.92	(hset[] for p_elim = 13)
88 *	234.16	(hset[] for p_elim = 17)
89 *	p_elim == 19 requires over 190 Megs of memory
90 *
91 * Over a long period of time, the call to load_hset() becomes insignificant.
92 * If we look at the user CPU seconds from the first 10000 cycle to the
93 * end of the test we find:
94 *
95 *	205.00	(prior to use of hset[])
96 *	 75.89	(hset[] for p_elim = 7)
97 *	 73.74	(hset[] for p_elim = 11)
98 *	 70.61	(hset[] for p_elim = 13)
99 *	 57.78	(hset[] for p_elim = 17)
100 *	 p_elim == 19 rejected because of memory size
101 *
102 * The p_elim == 17 overhead takes ~3 minutes on an 200 MHz r4k CPU and
103 * requires about ~13 Megs of memory.  The p_elim == 13 overhead
104 * takes about 3 seconds and requires ~1.5 Megs of memory.
105 *
106 * The value p_elim == 17 is best for long factorizations.  It is the
107 * fastest even thought the initial startup overhead is larger than
108 * for p_elim == 13.
109 *
110 * NOTE: The values above are prior to optimizations where hset[] was
111 *	 multiplied by n plus other optimizations.  Thus, the CPU
112 *	 times you may get will not likely match the above values.
113 */
114
115
116/*
117 * mfactor - find a factor of a Mersenne Number
118 *
119 * Mersenne numbers are numbers of the form:
120 *
121 *	2^n-1	for integer n > 0
122 *
123 * We know that factors of a Mersenne number are of the form:
124 *
125 *	2*k*n+1	  and	+/- 1 mod 8
126 *
127 * We make use of the hset[] difference array to eliminate factor
128 * candidates that would otherwise be divisible by 2, 3, 5, 7 ... p_elim.
129 *
130 * given:
131 *	n		attempt to factor M(n) = 2^n-1
132 *	start_k		the value k in 2*k*n+1 to start the search (def: 1)
133 *	rept_loop	loop cycle reporting (def: 10000)
134 *	p_elim		largest prime to eliminate from test factors (def: 17)
135 *
136 * returns:
137 *	factor of (2^n)-1
138 *
139 * NOTE: The p_elim argument is optional and defaults to 17.  A p_elim value
140 *	 of 17 is faster than 13 for even medium length runs.  However 13
141 *	 uses less memory and has a shorter startup time.
142 */
143define mfactor(n, start_k, rept_loop, p_elim)
144{
145	local Q;	/* 4*pfact(p_elim), hset[] cycle size */
146	local hcount;	/* elements in the hset[] difference array */
147	local loop;	/* report loop count */
148	local q;	/* test factor of 2^n-1 */
149	local g;	/* g as in test candidate form: (Q*g*hset[i])*n + 1 */
150	local hindx;	/* hset[] index */
151	local i;
152	local tmp;
153	local tmp2;
154
155	/*
156	 * firewall
157	 */
158	if (!isint(n) || n <= 0) {
159		quit "n must be an integer > 0";
160	}
161	if (!isint(start_k)) {
162		start_k = 1;
163	} else if (!isint(start_k) || start_k <= 0) {
164		quit "start_k must be an integer > 0";
165	}
166	if (!isint(rept_loop)) {
167		rept_loop = 10000;
168	}
169	if (rept_loop < 1) {
170		quit "rept_loop must be an integer > 0";
171	}
172	if (!isint(p_elim)) {
173		p_elim = 17;
174	}
175	if (p_elim < 3) {
176		quit "p_elim must be an integer > 2 (try 13 or 17)";
177	}
178
179	/*
180	 * declare our global values
181	 */
182	Q = 4*pfact(p_elim);
183	hcount = 2;
184	/* allocate the h difference array */
185	for (i=2; i <= p_elim; i = nextcand(i)) {
186		hcount *= (i-1);
187	}
188	local mat hset[hcount+1];
189
190	/*
191	 * load the hset[] difference array
192	 */
193	{
194		local x;	/* h*n+1 mod 8 */
195		local h;	/* potential h value */
196		local last_h;	/* previous valid h value */
197
198		last_h = 0;
199		for (i=0,h=0; h < Q; ++h) {
200			if (gcd(h*n+1,Q) == 1) {
201				x = (h*n+1) % 8;
202				if (x == 1 || x == 7) {
203					hset[i++] = (h-last_h) * n;
204					last_h = h;
205				}
206			}
207		}
208		hset[hcount] = Q*n - last_h*n;
209	}
210
211	/*
212	 * setup
213	 *
214	 * determine the next g and hset[] index (hindx) values such that:
215	 *
216	 *	2*start_k <= (Q*g + hset[hindx])
217	 *
218	 * and (Q*g + hset[hindx]) is a minimum and where:
219	 *
220	 *	Q = (4 * pfact(of some reasonable integer))
221	 *	g = (some integer) (hset[] cycle number)
222	 *
223	 * We also compute 'q', the next test candidate.
224	 */
225	g = (2*start_k) // Q;
226	tmp = 2*start_k - Q*g;
227	for (tmp2=0, hindx=0;
228	     hindx < hcount && (tmp2 += hset[hindx]/n) < tmp;
229	     ++hindx) {
230	}
231	if (hindx == hcount) {
232		/* we are beyond the end of a hset[] cycle, start at the next */
233		++g;
234		hindx = 0;
235		tmp2 = hset[0]/n;
236	}
237	q = (Q*g + tmp2)*n + 1;
238
239	/*
240	 * look for a factor
241	 *
242	 * We ignore factors that themselves are divisible by a prime <=
243	 * some small prime p.
244	 *
245	 * This process is guaranteed to find the smallest factor
246	 * of 2^n-1.  A smallest factor of 2^n-1 must be prime, otherwise
247	 * the divisors of that factor would also be factors of 2^n-1.
248	 * Thus we know that if a test factor itself is not prime, it
249	 * cannot be the smallest factor of 2^n-1.
250	 *
251	 * Eliminating all non-prime test factors would take too long.
252	 * However we can eliminate 80.81% of the test factors
253	 * by not using test factors that are divisible by a prime <= 17.
254	 */
255	if (pmod(2,n,q) == 1) {
256		return q;
257	} else {
258		/* report this loop */
259		printf("at 2*%d*%d+1, CPU: %f\n",
260			(q-1)/(2*n), n, usertime());
261		fflush(files(1));
262		loop = 0;
263	}
264	do {
265		/*
266		 * determine if we need to report
267		 *
268		 * NOTE: (q-1)/(2*n) is the k value from 2*k*n + 1.
269		 */
270		if (rept_loop <= ++loop) {
271			/* report this loop */
272			printf("at 2*%d*%d+1, CPU: %f\n",
273				(q-1)/(2*n), n, usertime());
274			fflush(files(1));
275			loop = 0;
276		}
277
278		/*
279		 * skip if divisible by a prime <= 449
280		 *
281		 * The value 281 was determined by timing loops
282		 * which found that 281 was at or near the
283		 * minimum time to factor 2^(2^127-1)-1.
284		 *
285		 * The addition of the do { ... } while (factor(q, 449)>1);
286		 * loop reduced the factoring loop time (36504 k values with
287		 * the hset[] initialization time removed) from 25.69 sec to
288		 * 15.62 sec of CPU time on a 200MHz r4k.
289		 */
290		do {
291			/*
292			 * determine the next factor candidate
293			 */
294			q += hset[++hindx];
295			if (hindx >= hcount) {
296				hindx = 0;
297				/*
298				 * if we cared about g,
299				 * then we wound ++g here too
300				 */
301			}
302		} while (factor(q, 449) > 1);
303	} while (pmod(2,n,q) != 1);
304
305	/*
306	 * return the factor found
307	 *
308	 * q is a factor of (2^n)-1
309	 */
310	return q;
311}
312
313if (config("resource_debug") & 3) {
314	print "mfactor(n [, start_k=1 [, rept_loop=10000 [, p_elim=17]]])"
315}
316