1/* 2 * ellip - attempt to factor numbers using elliptic functions 3 * 4 * Copyright (C) 1999 David I. Bell 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: 1990/02/15 01:50:33 21 * File existed as early as: before 1990 22 * 23 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ 24 */ 25 26/* 27 * Attempt to factor numbers using elliptic functions: 28 * 29 * y^2 = x^3 + a*x + b (mod ellip_N). 30 * 31 * Many points (x,y) (mod ellip_N) are found that solve the above equation, 32 * starting from a trivial solution and 'multiplying' that point together 33 * to generate high powers of the point, looking for such a point whose 34 * order contains a common factor with ellip_N. The order of the group of 35 * points varies almost randomly within a certain interval for each choice of 36 * a and b, and thus each choice provides an independent opportunity to 37 * factor ellip_N. To generate a trivial solution, a is chosen and then b is 38 * selected so that (1,1) is a solution. The multiplication is done using 39 * the basic fact that the equation is a cubic, and so if a line hits the 40 * curve in two rational points, then the third intersection point must 41 * also be rational. Thus by drawing lines between known rational points 42 * the number of rational solutions can be made very large. When modular 43 * arithmetic is used, solving for the third point requires the taking of a 44 * modular inverse (instead of division), and if this fails, then the GCD 45 * of the failing value and ellip_N provides a factor of ellip_N. 46 * This description is only an approximation, read "A Course in Number 47 * Theory and Cryptography" by Neal Koblitz for a good explanation. 48 * 49 * efactor(iN, ia, B, force) 50 * iN is the number to be factored. 51 * ia is the initial value of a in the equation, and each successive 52 * value of a is an independent attempt at factoring (default 1). 53 * B is the limit of the primes that make up the high power that the 54 * point is raised to for each factoring attempt (default 100). 55 * force is a flag to attempt to factor numbers even if they are 56 * thought to already be prime (default FALSE). 57 * 58 * Making B larger makes the power the point being raised to contain more 59 * prime factors, thus increasing the chance that the order of the point 60 * will be made up of those factors. The higher B is then, the greater 61 * the chance that any individual attempt will find a factor. However, 62 * a higher B also slows down the number of independent functions being 63 * examined. The order of the point for any particular function might 64 * contain a large prime and so won't succeed even for a really large B, 65 * whereas the next function might have an order which is quickly found. 66 * So you want to trade off the depth of a particular search with the 67 * number of searches made. For example, for factoring 30 digits, I make 68 * B be about 1000 (probably still too small). 69 * 70 * If you have lots of machines available, then you can run parallel 71 * factoring attempts for the same number by giving different starting 72 * values of ia for each machine (e.g. 1000, 2000, 3000). 73 * 74 * The output as the function is running is (occasionally) the value of a 75 * when a new function is started, the prime that is being included in the 76 * high power being calculated, and the current point which is the result 77 * of the powers so far. 78 * 79 * If a factor is found, it is returned and is also saved in the global 80 * variable f. The number being factored is also saved in the global 81 * variable ellip_N. 82 */ 83 84 85obj point {x, y}; 86global ellip_N; /* number to factor */ 87global ellip_a; /* first coefficient */ 88global ellip_b; /* second coefficient */ 89global ellip_f; /* found factor */ 90 91 92define efactor(iN, ia, B, force) 93{ 94 local C, x, p; 95 96 if (!force && ptest(iN, 50)) 97 return 1; 98 if (isnull(B)) 99 B = 100; 100 if (isnull(ia)) 101 ia = 1; 102 obj point x; 103 ellip_a = ia; 104 ellip_b = -ia; 105 ellip_N = iN; 106 C = isqrt(ellip_N); 107 C = 2 * C + 2 * isqrt(C) + 1; 108 ellip_f = 0; 109 while (ellip_f == 0) { 110 print "A =", ellip_a; 111 x.x = 1; 112 x.y = 1; 113 print 2, x; 114 x = x ^ (2 ^ (highbit(C) + 1)); 115 for (p = 3; ((p < B) && (ellip_f == 0)); p += 2) { 116 if (!ptest(p, 1)) 117 continue; 118 print p, x; 119 x = x ^ (p ^ ((highbit(C) // highbit(p)) + 1)); 120 } 121 ellip_a++; 122 ellip_b--; 123 } 124 return ellip_f; 125} 126 127 128define point_print(p) 129{ 130 print "(" : p.x : "," : p.y : ")" :; 131} 132 133 134define point_mul(p1, p2) 135{ 136 local r, m; 137 138 if (p2 == 1) 139 return p1; 140 if (p1 == p2) 141 return point_square(`p1); 142 obj point r; 143 m = (minv(p2.x - p1.x, ellip_N) * (p2.y - p1.y)) % ellip_N; 144 if (m == 0) { 145 if (ellip_f == 0) 146 ellip_f = gcd(p2.x - p1.x, ellip_N); 147 r.x = 1; 148 r.y = 1; 149 return r; 150 } 151 r.x = (m^2 - p1.x - p2.x) % ellip_N; 152 r.y = ((m * (p1.x - r.x)) - p1.y) % ellip_N; 153 return r; 154} 155 156 157define point_square(p) 158{ 159 local r, m; 160 161 obj point r; 162 m = ((3 * p.x^2 + ellip_a) * minv(p.y << 1, ellip_N)) % ellip_N; 163 if (m == 0) { 164 if (ellip_f == 0) 165 ellip_f = gcd(p.y << 1, ellip_N); 166 r.x = 1; 167 r.y = 1; 168 return r; 169 } 170 r.x = (m^2 - p.x - p.x) % ellip_N; 171 r.y = ((m * (p.x - r.x)) - p.y) % ellip_N; 172 return r; 173} 174 175 176define point_pow(p, pow) 177{ 178 local bit, r, t; 179 180 r = 1; 181 if (isodd(pow)) 182 r = p; 183 t = p; 184 for (bit = 2; ((bit <= pow) && (ellip_f == 0)); bit <<= 1) { 185 t = point_square(`t); 186 if (bit & pow) 187 r = point_mul(`t, `r); 188 } 189 return r; 190} 191