![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to lfactor.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / lfactor|
Return to lfactor.c CVS log | Up to [Research Unix] / researchv10no / cmd / lfactor |
1.1 root 1: #include <stdio.h>
2: #include "mp.h"
3: #include "form.h"
4: #define testing 0
5: #define fmax 20
6: #define gridmax 30000
7: form stockf[fmax];
8: struct grid { short a; float b;};
9: union{
10: struct grid grida[gridmax+1];
11: form powers[512];
12: }ugly;
13: FILE *ptab;
14:
15: mint factab[100];
16: int facptr;
17: double flimit;
18:
19: main(argc, argp)
20: int argc;
21: char *argp[];
22: {
23: form formt1;
24: form formtt;
25: form bigf, bigfi;
26: form smallf;
27: form formttt;
28: form regen;
29: unsigned giant;
30: struct grid *j1, *j2;
31: int ordind;
32: int prime;
33: mint nn;
34: mint nugget;
35: mint mjunk, mtemp;
36: mint mtemp1;
37: mint zero, one, two, four;
38: double a, b;
39: mint ma, mb, mc;
40: double temp;
41: mint det;
42: double lastpr;
43: double m,n;
44: double prod;
45: double hbar, h;
46: double pi = 3.141592653589793;
47: double ht, hr, it;
48: double logord;
49: double maxord;
50: double order[fmax];
51: double tt, ttt;
52: float ftt, fttt;
53: int baby;
54: int kbaby;
55: register j;
56: int i;
57: double k;
58: int goti, gotj, gotit;
59: int ttn;
60: int junk;
61: int uneq;
62: int jmax;
63: int fcount;
64: int parity;
65: double sqroot, fifth;
66: int verbose = 0;
67: int nofact = 0;
68: int compare();
69: int compar1();
70: double getpr();
71: double modf();
72: double log(), exp(), sqrt();
73:
74: setbuf(stdout,NULL);
75: verbose = 0;
76: if(argc > 1){
77: if(*argp[1] == '+'){
78: verbose = 1;
79: }
80: if(*argp[1] == '-'){
81: nofact = 1;
82: }
83: }
84: zero = itom(0);
85: one = itom(1);
86: two = itom(2);
87: four = itom(4);
88: facptr = 0;
89:
90: if(min(&nn) == EOF) exit(0);
91: if(mcmp(nn,zero) <= 0) exit(0);
92: while(1){
93: mtemp = sdiv(nn,2,&i);
94: if(i==0){
95: printf("Prime factor: 2\n");
96: nn = mtemp;
97: }else break;
98: }
99:
100: ptab = fopen("/usr/lib/ptab","r");
101: if(ptab == NULL){
102: printf("ptab?\n");
103: exit(1);
104: }
105: retry:
106: getpr(-1);
107: lastpr = 2.;
108: if(mcmp(nn,itom(1))==0)
109: exit(0);
110: sqroot = sqrt(nn.high*1e16 + nn.low);
111: fifth = exp(0.2*log(nn.high*1e16+nn.low));
112: if(fifth < 1000.) fifth = 1000.;
113: sdiv(nn,4, &i);
114: if((i== 1) || (i== 2)){
115: parity = 0;
116: prod = 1.;
117: }
118: if(i== 3)
119: parity = 1;
120: sdiv(nn, 8, &junk);
121: if(junk == 3) prod = 2./3.;
122: if(junk == 7) prod = 2.;
123: fcount = 0;
124: while(1){
125: if((m = getpr(0)) < 0)
126: break;
127: lastpr = m;
128: if(lastpr > sqroot){
129: printf("Prime factor: %.0f\n", nn.low);
130: exit(0);
131: }
132: if(lastpr > 10.*fifth)
133: break;
134: modf(nn.high/m, &temp);
135: n = nn.high - m*temp;
136: modf(e16/m, &temp);
137: temp = e16 - m*temp;
138: n = n*temp;
139: modf(nn.low/m, &temp);
140: temp = nn.low - m*temp;
141: n = n + temp;
142: j = jacobi(-n,m);
143: if((j==0) && (nofact==0)){
144: printf("Prime factor: %.0f\n", m);
145: mtemp.high = 0;
146: mtemp.low = m;
147: nn = mdiv(nn, mtemp, &mtemp);
148: if(mcmp(mtemp,zero)!=0) abort();
149: goto retry;
150: }
151: a = m;
152: if((j==1) && (fcount<fmax)){
153: det = mchs(nn);
154: sdiv(det, 4, &i);
155: if((i== -1) || (i== -2))
156: det = smult(det, 4);
157: for(b=parity; b<=a; b = b+2){
158: ma.high = 0.;
159: ma.low = a;
160: mb.high = 0.;
161: mb.low = b;
162: mc.high = 0.;
163: mc.low = b*b;
164: mc = msub(mc, det);
165: mtemp.high = 0.;
166: mtemp.low = 4.*a;
167: mc = mdiv(mc, mtemp, &mtemp);
168: if(mcmp(mtemp,zero) == 0){
169: stockf[fcount].a = ma;
170: stockf[fcount].b = mb;
171: stockf[fcount].c = mc;
172: #if testing
173: if(verbose){
174: formout(stockf[fcount]);
175: printf("\n");
176: }
177: #endif
178: fcount++;
179: }
180: }
181: }
182: prod *= m/(m-j);
183: }
184: printf("Trying to factor: n = ");
185: mout(nn);
186: printf("Factor limit = %.0f\n", a);
187: flimit = a;
188: mtemp1 = mcbrt(nn, &mjunk);
189: if(mjunk.low == 0){
190: mout1(nn);
191: printf(" is a cube.\n");
192: exit(0);
193: }
194: mtemp1 = msqrt(nn, &mjunk);
195: if(mjunk.low == 0){
196: mout1(nn);
197: printf(" is a square.\n");
198: exit(0);
199: }
200:
201: prime = pseudo(nn, itom(2));
202: if(prime == 0) prime = pseudo(nn, itom(3));
203: if(prime == 0) prime = pseudo(nn, itom(5));
204: if(prime == 0) prime = pseudo(nn, itom(7));
205: if(prime == 0){
206: printf("n is pseudoprime (mod 2,3,5,7).\n");
207: }else{
208: printf("n is composite.\n");
209: }
210:
211: if(parity==0) prod *= 2.;
212: modf(prod*mtemp1.low/pi, &hbar);
213: #if testing
214: printf("h~~ %.0f\n", hbar);
215: #endif
216:
217: ugly.grida[0].a = 0;
218: ugly.grida[0].b = 1.;
219: if(parity == 1){
220: baby = 1;
221: kbaby = 0;
222: }else{
223: baby = 2;
224: kbaby = 0;
225: }
226:
227: sdiv(nn, 8, &i);
228: if(i == 1){
229: baby = 4;
230: kbaby = 0;
231: }else if((i == 5) && (prime != 0)){
232: baby = 4;
233: kbaby = 0;
234: }
235: if((parity == 1) && (prime != 0)){
236: baby = 2;
237: kbaby = 0;
238: }
239:
240: modf(hbar/baby, &temp);
241: hbar = baby*temp;
242: hbar = hbar + kbaby;
243:
244: smallf = formpow(stockf[0], (double)baby);
245: ugly.grida[1].a = 1;
246: ugly.grida[1].b = smallf.a.low;
247: formt1 = formpow(stockf[0], hbar);
248:
249: if(mcmp(formt1.a,one) == 0){
250: h = hbar;
251: goto found;
252: }
253: formtt = smallf;
254: ttn = 1;
255: temp = sqrt(0.1*hbar/(sqrt(fifth)*baby));
256: if(temp > gridmax) temp = gridmax;
257: giant = temp;
258: printf("Grid = %d by %d\n", giant, baby);
259: while(ttn <= giant){
260: if(formtt.a.low == formt1.a.low){
261: if(formtt.b.low == formt1.b.low){
262: h = hbar - ttn*baby;
263: goto found;
264: }
265: if(formtt.b.low == -formt1.b.low){
266: h = hbar + ttn*baby;
267: goto found;
268: }
269: }
270: formtt = compos(formtt, smallf);
271: ttn += 1;
272: ugly.grida[ttn].a = ttn;
273: ugly.grida[ttn].b = formtt.a.low;
274: }
275: qsort(&ugly.grida[0].a, giant+1, sizeof(ugly.grida[0]), compare);
276: bigf = formpow(smallf, (double)(2*giant));
277: bigfi = bigf;
278: bigfi.b.low = -bigfi.b.low;
279: formtt = compos(bigf, formt1);
280: formttt = compos(bigfi, formt1);
281: for(k=1; k<fifth/2.; k=k+1){
282: tt = formtt.a.low;
283: ttt = formttt.a.low;
284: ftt = formtt.a.low;
285: fttt = formttt.a.low;
286: if(search(&ugly.grida[0].a, giant+1, sizeof(ugly.grida[0]), compar1, &ftt, &j1, &j2) >=0){
287: b1:
288: j = j1->a;
289: regen = formpow(smallf, (double)j);
290: if(tt == regen.a.low){
291: if(formtt.b.low == regen.b.low){
292: h = hbar + (2.*giant*k - j)*baby;
293: goto found;
294: }else
295: if(formtt.b.low == -regen.b.low){
296: h = hbar + (2.*giant*k + j)*baby;
297: goto found;
298: }
299: }
300: if(j2 > j1+1){
301: j1 = j1 + 1;
302: goto b1;
303: }
304: }
305: if(search(&ugly.grida[0].a, giant+1, sizeof(ugly.grida[0]), compar1, &fttt, &j1, &j2) >=0){
306: b2:
307: j = j1->a;
308: regen = formpow(smallf, (double)j);
309: if(ttt == regen.a.low){
310: if(formttt.b.low == regen.b.low){
311: h = hbar - (2.*giant*k + j)*baby;
312: goto found;
313: }else
314: if(formttt.b.low == -regen.b.low){
315: h = hbar - (2.*giant*k - j)*baby;
316: goto found;
317: }
318: }
319: if(j2 > j1 + 1){
320: j1 = j1 + 1;
321: goto b2;
322: }
323: }
324: formtt = compos(formtt, bigf);
325: formttt = compos(formttt, bigfi);
326: }
327: printf("Search abandoned at %.0f.\n", 2.*giant*k*baby);
328: exit(1);
329: found:
330: printf("h = %.0f\n", h);
331: temp = hbar - h;
332: printf("Search succeeded at %.0f (%.3f) ", temp, 1000.*temp/h);
333: printf("(%.3f)\n",flimit*temp*temp/(h*h));
334:
335: ht = 0;
336: hr =h;
337: while(1){
338: modf(hr/2, &temp);
339: if(hr == 2.*temp){
340: hr = temp;
341: ht += 1;
342: }else{
343: break;
344: }
345: }
346: /*
347: * h = hr * 2^ht
348: */
349: if(ht == 0){
350: printf("Class number is odd!\n");
351: mout1(nn);
352: printf(" is a prime.\n");
353: exit(0);
354: }
355: #if testing
356: printf("Sylow 2-subgroup of order 2^%.0f\n", ht);
357: #endif
358: for(i=0; i<fmax; i++){
359: stockf[i] = formpow(stockf[i], hr);
360: }
361:
362: for(i=0; i<fmax; i++){
363: if(mcmp(stockf[i].b,zero) < 0){
364: stockf[i].b = mchs(stockf[i].b);
365: }
366: }
367:
368: for(i=1,jmax=1; i<fmax; i++){
369: for(j=0; j<jmax; j++){
370: uneq = 0;
371: if(mcmp(stockf[i].a,stockf[j].a) != 0) uneq = 1;
372: if(mcmp(stockf[i].b,stockf[j].b) != 0) uneq = 1;
373: if(mcmp(stockf[i].c,stockf[j].c) != 0) uneq = 1;
374: if(uneq == 0) break;
375: }
376: if(uneq == 1){
377: stockf[jmax++] = stockf[i];
378: }
379: }
380:
381: for(i=0; i<jmax; i++){
382: it = 0;
383: formt1 = stockf[i];
384: if(mcmp(formt1.a,one) == 0) continue;
385: while(1){
386: formtt = formpow(formt1, (double)2);
387: it = it+1;
388: if(it>ht){
389: printf("main: error no. 1\n");
390: abort();
391: }
392: if(mcmp(formtt.a,one) == 0) break;
393: formt1 = formtt;
394: }
395: order[i] = it;
396: /*
397: printf("f%3d: ", i);
398: formout(stockf[i]);
399: printf(" order = 2^%.0f: ", it);
400: formout(formt1);
401: printf("\n");
402: */
403: }
404: logord = 0;
405: ordind = 0;
406: for(i=0; i<jmax; i++){
407: if(order[i] > logord){
408: ordind = i;
409: logord = order[i];
410: }
411: }
412:
413: if(logord == ht){
414: printf("Sylow 2-subgroup is cyclic.\n");
415: formtt = stockf[ordind];
416: for(i=0; i<order[ordind]-1; i++){
417: formtt = formpow(formtt, (double)2);
418: }
419: if(mcmp(formtt.a,two) == 0){
420: mout1(nn);
421: printf(" is a prime.\n");
422: exit(0);
423: }
424: if(mcmp(formtt.b,zero) == 0){
425: sqtest(formtt.a, flimit);
426: sqtest(formtt.c, flimit);
427: exit(0);
428: }
429: if(mcmp(formtt.a, formtt.b) == 0){
430: sqtest(formtt.a, flimit);
431: sqtest(msub(mult(itom(4),formtt.c),formtt.a), flimit);
432: exit(0);
433: }
434: if(mcmp(formtt.a, formtt.c) == 0){
435: sqtest(msub(mult(itom(2),formtt.a),formtt.b), flimit);
436: sqtest(madd(mult(itom(2),formtt.a),formtt.b), flimit);
437: exit(0);
438: }
439: printf("main: error no. 4\n");
440: exit(1);
441: }
442: if(logord > 10){
443: printf("Sylow 2-subgroup too big.\n");
444: exit(1);
445: }
446:
447: for(i=0,j=1;i<logord;i++){
448: j = j*2;
449: }
450: maxord = j;
451: printf("\n");
452:
453: ugly.powers[0] = stockf[ordind];
454: for(i=1;i<maxord/2;i++){
455: ugly.powers[i] = compos(ugly.powers[i-1],stockf[ordind]);
456: }
457:
458: for(i=0; i<jmax; i++){
459: if(i == ordind) continue;
460: gotit = 0;
461: formt1 = stockf[i];
462: for(j=0; j<=order[i]; j++){
463: for(junk=0; junk<maxord/2; junk++){
464: if(mcmp(formt1.a,ugly.powers[junk].a) != 0)
465: continue;
466: if(mcmp(formt1.c,ugly.powers[junk].c) != 0)
467: continue;
468: if(mcmp(formt1.b,ugly.powers[junk].b) != 0){
469: goti = junk + 1;
470: }else{
471: goti = maxord - junk - 1;
472: }
473: for(junk=0,gotj=1; junk<j; junk++){
474: gotj *= 2;
475: }
476: if(goti%gotj != 0){
477: printf("main: error no. 2\n");
478: }
479: goti = goti/gotj;
480: if(goti > maxord/2){
481: goti = maxord - goti;
482: }
483: stockf[i] = compos(stockf[i],ugly.powers[goti-1]);
484: gotit = 1;
485: break;
486: }
487: if(gotit == 1) break;
488: formt1 = compos(formt1,formt1);
489: }
490: }
491:
492: for(i=0; i<jmax; i++){
493: if(mcmp(stockf[i].b,zero) < 0){
494: stockf[i].b = mchs(stockf[i].b);
495: }
496: }
497:
498: junk = jmax;
499: for(i=1,jmax=1; i<junk; i++){
500: for(j=0; j<jmax; j++){
501: uneq = 0;
502: if(mcmp(stockf[i].a,stockf[j].a) != 0) uneq = 1;
503: if(mcmp(stockf[i].b,stockf[j].b) != 0) uneq = 1;
504: if(mcmp(stockf[i].c,stockf[j].c) != 0) uneq = 1;
505: if(uneq == 0) break;
506: }
507: if(uneq == 1){
508: stockf[jmax++] = stockf[i];
509: }
510: }
511:
512: for(i=0; i<jmax; i++){
513: it = 0;
514: formt1 = stockf[i];
515: if(mcmp(formt1.a,one) == 0) continue;
516: while(1){
517: formtt = formpow(formt1, (double)2);
518: it = it+1;
519: if(it>ht){
520: printf("main: error no. 3\n");
521: abort();
522: }
523: if(mcmp(formtt.a,one) == 0) break;
524: formt1 = formtt;
525: }
526: order[i] = it;
527: if(verbose){
528: printf("f%3d: ", i);
529: formout(stockf[i]);
530: printf(" order = 2^%.0f: ", it);
531: formout(formt1);
532: printf("\n");
533: }
534: if(mcmp(formt1.b,zero) == 0){
535: putfact(formt1.a);
536: putfact(formt1.c);
537: }else if(mcmp(formt1.a,formt1.b) == 0){
538: putfact(formt1.a);
539: putfact(msub(mult(four,formt1.c),formt1.a));
540: }else if(mcmp(formt1.a,formt1.c) == 0){
541: putfact(madd(mult(two,formt1.a),formt1.b));
542: putfact(msub(mult(two,formt1.a),formt1.b));
543: }else{
544: printf("main: error no. 5\n");
545: }
546: }
547: nugget = nn;
548: for(i=0;i<facptr;i++){
549: if(mcmp(factab[i],one) == 0) continue;
550: if(mcmp(factab[i],nn) >= 0) continue;
551: if(mcmp(factab[i],nugget) > 0) continue;
552: prtest(factab[i], flimit);
553: nugget = mdiv(nugget, factab[i], &mjunk);
554: if(mcmp(mjunk,zero) != 0){
555: printf("main: error no. 6\n");
556: }
557: }
558: if(mcmp(nugget,one) > 0){
559: prtest(nugget, flimit);
560: }
561: exit(0);
562: }
563:
564: /*
565: * m must be positive and odd.
566: */
567: int
568: jacobi(n,m)
569: double n,m;
570: {
571: int mr2, mrm1, i, j, sign;
572: double temp;
573: double modf();
574:
575: sign = 1;
576: if(n==0) return(0);
577: if(n==1) return(1);
578: if(m==1) return(1);
579:
580: mr2 = 1;
581: mrm1 = 1;
582: modf(m/8, &temp);
583: i = m - 8*temp;
584: switch(i){
585: case 3:
586: mrm1 = -1;
587: case 5:
588: mr2 = -1;
589: case 1:
590: break;
591: case 7:
592: mrm1 = -1;
593: break;
594: default:
595: abort();
596: }
597:
598: modf(n/m, &temp);
599: n = n - m*temp;
600: if(n<0) n += m;
601: if(n==0) return(0);
602: if(n+n>m){
603: n = m-n;
604: sign = sign*mrm1;
605: }
606: while(1){
607: if(modf(.5*n, &temp) != 0) break;
608: n = temp;
609: sign = sign * mr2;
610: }
611: if(m>=32768.)
612: return(sign * jacobi(mrm1*m,n));
613: else{
614: i = mrm1*m;
615: j = n;
616: return(sign * ijac(i,j));
617: }
618: }
619: /*
620: * m must be positive and odd.
621: */
622: int
623: ijac(n,m)
624: int n,m;
625: {
626: int mr2, mrm1, sign;
627:
628: sign = 1;
629: if(n==0) return(0);
630: if(n==1) return(1);
631: if(m==1) return(1);
632:
633: mr2 = mrm1 = 1;
634: switch(m&07){
635: case 3:
636: mrm1 = -1;
637: case 5:
638: mr2 = -1;
639: case 1:
640: break;
641: case 7:
642: mrm1 = -1;
643: break;
644: default:
645: abort();
646: }
647:
648: n = n%m;
649: if(n<0) n += m;
650: if(n==0) return(0);
651: if((n&01)==1){
652: n = m-n;
653: sign = sign*mrm1;
654: }
655: while((n&03) == 0)
656: n = n>>2;
657: if((n&01) == 0){
658: n = n>>1;
659: sign = sign * mr2;
660: }
661: return(sign * ijac(mrm1*m,n));
662: }
663:
664: int sieve[] = {
665: 1, 7, 11, 13, 17, 19, 23, 29,
666: 31, 37, 41, 43, 47, 49, 53, 59,
667: 61, 67, 71, 73, 77, 79, 83, 89,
668: 91, 97,101,103,107,109,113,119,
669: };
670: double
671: getpr(arg)
672: int arg;
673: {
674: static int i;
675: static unsigned word;
676: static double base;
677:
678: if(arg<0){
679: rewind(ptab);
680: word = getw(ptab);
681: if(feof(ptab)) exit(0);
682: if(ferror(ptab)) exit(0);
683: base = 0;
684: i = -2;
685: return(0.);
686: }
687: if((base == 0) && (i < 0)){
688: if(i == -2){
689: i++;
690: return(3);
691: }
692: if(i == -1){
693: i++;
694: return(5);
695: }
696: }
697: while(1){
698: i++;
699: if(i>(8*sizeof(int))){
700: word = getw(ptab);
701: if(ferror(ptab)) abort();
702: if(feof(ptab)) return((double)-1);
703: i -= 8*sizeof(int);
704: base += 30*sizeof(int);
705: }
706: if((word & (1<<(i-1))) != 0)
707: return(base+sieve[i-1]);
708: }
709: }
710:
711: int
712: compare(a,b)
713: struct grid *a, *b;
714: {
715: if(a->b > b->b)
716: return(1);
717: if(a->b < b->b)
718: return(-1);
719: return(0);
720: }
721: int
722: compar1(a,b)
723: float *a;
724: struct grid *b;
725: {
726: if(*a > b->b)
727: return(1);
728: if(*a < b->b)
729: return(-1);
730: return(0);
731: }
732:
733: void
734: sqtest(nn, flimit)
735: mint nn;
736: double flimit;
737: {
738:
739: mint temp1, temp2;
740: double dtemp;
741:
742: temp1.low = flimit;
743: temp1.high = 0;
744: if(mcmp(temp1,nn) > 0){
745: printf("sqtest: error no. 1\n");
746: exit(1);
747: }
748:
749: dtemp = log(nn.high+nn.low)/log(flimit);
750: if(dtemp >= 5.0){
751: printf("sqtest: error no. 2");
752: exit(1);
753: }
754:
755: temp1 = msqrt(nn, &temp2);
756: if(temp2.low == 0){
757: printf("Square factor: ");
758: mout(temp1);
759: return;
760: }
761: temp1 = mcbrt(nn, &temp2);
762: if(temp2.low == 0){
763: printf("Cubic factor: ");
764: mout(temp1);
765: return;
766: }
767: printf("Prime factor: ");
768: mout(nn);
769: return;
770: }
771: putfact(arg)
772: mint arg;
773: {
774: mint mtemp;
775: int i;
776:
777:
778: while(1){
779: mtemp = sdiv(arg,2,&i);
780: if(i==0){
781: arg = mtemp;
782: }
783: else break;
784: }
785: if(facptr == 0){
786: factab[facptr++] = arg;
787: return(0);
788: }
789: for(i=0; i<facptr; i++){
790: if(mcmp(factab[i],arg) < 0) continue;
791: if(mcmp(factab[i],arg) == 0) return(0);
792: mtemp = factab[i];
793: factab[i] = arg;
794: arg = mtemp;
795: continue;
796: }
797: factab[facptr++] = arg;
798: return;
799:
800: }
801:
802: prtest(arg)
803: mint arg;
804: {
805: double dtemp;
806:
807: dtemp = log(arg.high+arg.low)/log(flimit);
808: if(dtemp < 2.0){
809: printf("Prime factor: ");
810: }else{
811: printf("Factor: ");
812: }
813: mout(arg);
814: return;
815: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.