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