![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to mulmod.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 mulmod.c CVS log | Up to [Research Unix] / researchv10no / cmd / lfactor |
1.1 root 1: #include "mp.h"
2: mint mulmod();
3:
4: /*
5: returns 0 if pseudoprime
6: 1 if composite
7: 2 if composite pseudoprime
8: */
9: int
10: pseudo(nn, arg2)
11: mint nn;
12: mint arg2;
13: {
14: mint nm1, y, z;
15: mint zero, one, two;
16: mint rem;
17: mint mtemp;
18: mint lasty;
19: int i;
20:
21: zero.low = zero.high = 0.0;
22: one = itom(1);
23: two = itom(2);
24:
25: if(mcmp(nn,zero) <= 0) abort(0);
26: if(mcmp(nn,two) <=0) abort(0);
27: mdiv(nn, two, &rem);
28: if(mcmp(rem, zero) == 0.0) abort(0);
29:
30: nm1 = msub(nn, one);
31: y = one;
32: z = arg2;
33:
34: i = 0;
35: while(1){
36: mtemp = mdiv(nm1, two, &rem);
37: if(mcmp(rem,zero) == 0){
38: nm1 = mtemp;
39: i = i + 1;
40: }else break;
41: }
42:
43: while(mcmp(nm1,zero) != 0){
44: nm1 = mdiv(nm1, two, &rem);
45: if(mcmp(rem, one) == 0){
46: y = mulmod(y, z, nn);
47: }
48: z = mulmod(z, z, nn);
49: }
50: if(mcmp(y, one) == 0){
51: return(0);
52: }
53: while(i--){
54: lasty = y;
55: y = mulmod(y,y,nn);
56: if(mcmp(y, one) == 0) break;
57: }
58: if(mcmp(y, one) != 0){
59: return(1);
60: }
61: if(mcmp(lasty, one) == 0){
62: printf("pseudo: error.\n");
63: }
64: if(mcmp(lasty,msub(nn,one)) != 0){
65: return(2);
66: }
67: return(0);
68: }
69: mint
70: mulmod(a,b,c)
71: mint a, b, c;
72: {
73:
74: mint mjunk;
75: mint mtemp1, mtemp2;
76: mint xy0, xy1, xy2, xy3;
77: mint x1, x2, y1, y2;
78: int i;
79: double z0, z1, z2, z3;
80: double s0, s1, s2;
81: double tquot, dtemp1, dtemp2;
82: mint zero;
83:
84:
85: zero.low = 0.0;
86: zero.high = 0.0;
87: if(mcmp(c,zero) == 0){
88: printf("mulmod: divide check\n");
89: abort(0);
90: }
91:
92: if((mcmp(a,zero)<0) || (mcmp(b,zero)<0)){
93: printf("mulmod: negative argument.\n");
94: abort(0);
95: }
96: if(mcmp(c,zero) < 0){
97: printf("mulmod: negative modulus.\n");
98: abort(0);
99: }
100: x1.low = a.low;
101: x1.high = 0;
102: x2.low = a.high;
103: x2.high = 0;
104:
105: y1.low = b.low;
106: y1.high = 0;
107: y2.low = b.high;
108: y2.high = 0;
109:
110: xy0 = mult(x1,y1);
111: xy1 = mult(x1,y2);
112: xy2 = mult(x2,y1);
113: xy3 = mult(x2,y2);
114:
115: z0 = xy0.low;
116:
117: mtemp1.low = xy0.high;
118: mtemp1.high = 0;
119: mtemp2.low = xy1.low;
120: mtemp2.high = 0;
121: mtemp1 = madd(mtemp1, mtemp2);
122: mtemp2.low = xy2.low;
123: mtemp1 = madd(mtemp1, mtemp2);
124: z1 = mtemp1.low;
125:
126: mtemp1.low = mtemp1.high;
127: mtemp1.high = 0;
128: mtemp2.low = xy1.high;
129: mtemp2.high = 0;
130: mtemp1 = madd(mtemp1, mtemp2);
131: mtemp2.low = xy2.high;
132: mtemp1 = madd(mtemp1, mtemp2);
133: mtemp2.low = xy3.low;
134: mtemp1 = madd(mtemp1, mtemp2);
135: z2 = mtemp1.low;
136:
137: mtemp1.low = mtemp1.high;
138: mtemp1.high = 0;
139: mtemp2.low = xy3.high;
140: mtemp2.high = 0;
141: mtemp1 = madd(mtemp1, mtemp2);
142: z3 = mtemp1.low;
143:
144:
145: if(mtemp1.high != 0.0){
146: printf("mulmod: error\n");
147: }
148:
149: if(c.high == 0.0){
150: mtemp1.high = 0.0;
151: mtemp1.low = z3;
152: mjunk = mdiv(mtemp1, c, &mtemp2);
153: z3 = mtemp2.low;
154: if(mtemp2.high != 0.0) trouble(12);
155: mtemp1.high = z3;
156: mtemp1.low = z2;
157: mjunk = mdiv(mtemp1, c, &mtemp2);
158: mtemp1.high = mtemp2.low;
159: mtemp1.low = z1;
160: mjunk = mdiv(mtemp1, c, &mtemp2);
161: mtemp1.high = mtemp2.low;
162: mtemp1.low = z0;
163: mjunk = mdiv(mtemp1, c, &mtemp2);
164: if(mcmp(mtemp2, c) >= 0) trouble(10);
165: if(mcmp(mtemp2, zero) < 0) trouble(11);
166: return(mtemp2);
167: }
168:
169: mtemp1.high = z3;
170: mtemp1.low = z2;
171: if(mcmp(mtemp1, c) >= 0){
172: mjunk = mdiv(mtemp1, c, &mtemp2);
173: z3 = mtemp2.high;
174: z2 = mtemp2.low;
175: }
176: dtemp1 = (z3*e16 + z2) + z1/e16;
177: tquot = dtemp1/(c.high+c.low/e16);
178: modf(tquot, &tquot);
179: y1.low = tquot;
180: y1.high = 0.0;
181: x1.low = c.low;
182: x1.high = 0.0;
183: x2.low = c.high;
184: x2.high = 0;
185: xy0 = mult(x1, y1);
186: xy1 = mult(x2, y1);
187:
188: s0 = xy0.low;
189:
190: mtemp1.low = xy0.high;
191: mtemp1.high = 0.0;
192: mtemp1 = madd(mtemp1, xy1);
193: s1 = mtemp1.low;
194: s2 = mtemp1.high;
195: s0 = z1 - s0;
196: s1 = z2 - s1;
197: s2 = z3 - s2;
198:
199: if(s2 > 0.0) {s2 = s2 - 1; s1 = s1 + e16;}
200: if(s2 < 0.0) {s2 = s2 + 1; s1 = s1 - e16;}
201: if((s1<0.0) && (s0>0.0)){
202: s1 = s1 + 1;
203: s0 = s0 - e16;
204: }else
205: if((s1>0.0) && (s0<0.0)){
206: s1 = s1 - 1;
207: s0 = s0 + e16;
208: }
209: if((s1*e16+s0) < 0.0){
210: mtemp1.low = s0;
211: mtemp1.high = s1;
212: mtemp1 = madd(mtemp1, c);
213: s1 = mtemp1.high;
214: s0 = mtemp1.low;
215: }
216: if(s2 != 0.0) trouble(1);
217: mtemp1.low = s0;
218: mtemp1.high = s1;
219: if(mcmp(mtemp1,zero) < 0) trouble(2);
220: if(mcmp(mtemp1,c) >= 0){
221: mtemp1.high = s1;
222: mtemp1.low = s0;
223: mtemp1 = msub(mtemp1, c);
224: s1 = mtemp1.high;
225: s0 = mtemp1.low;
226: }
227: if(mcmp(mtemp1,c) >0) trouble(3);
228:
229: z3 = s2;
230: z2 = s1;
231: z1 = s0;
232:
233: dtemp1 = (z2*e16 + z1) + z0/e16;
234: tquot = dtemp1/(c.high + c.low/e16);
235: modf(tquot, &tquot);
236: y1.low = tquot;
237: y1.high = 0.0;
238: x1.low = c.low;
239: x1.high = 0.0;
240: x2.low = c.high;
241: x2.high = 0.0;
242: xy0 = mult(x1, y1);
243: xy1 = mult(x2, y1);
244:
245: s0 = xy0.low;
246:
247: mtemp1.low = xy0.high;
248: mtemp1.high = 0.0;
249: mtemp1 = madd(mtemp1, xy1);
250: s1 = mtemp1.low;
251: s2 = mtemp1.high;
252:
253: s0 = z0 - s0;
254: s1 = z1 - s1;
255: s2 = z2 - s2;
256: if(s2 > 0.0) {s2 = s2 - 1; s1 = s1 + e16;}
257: if(s2 < 0.0) {s2 = s2 + 1; s1 = s1 - e16;}
258: if((s1<0.0) && (s0>0.0)){
259: s1 = s1 + 1;
260: s0 = s0 - e16;
261: }else
262: if((s1>0.0) && (s0<0.0)){
263: s1 = s1 - 1;
264: s0 = s0 + e16;
265: }
266: if((s1*e16+s0) < 0.0){
267: mtemp1.low = s0;
268: mtemp1.high = s1;
269: mtemp1 = madd(mtemp1, c);
270: s1 = mtemp1.high;
271: s0 = mtemp1.low;
272: }
273: if(s2 != 0.0) trouble(4);
274: mtemp1.low = s0;
275: mtemp1.high = s1;
276: if(mcmp(mtemp1,zero) < 0) trouble(5);
277: if(mcmp(mtemp1,c) >= 0){
278: mtemp1.high = s1;
279: mtemp1.low = s0;
280: mtemp1 = msub(mtemp1, c);
281: s1 = mtemp1.high;
282: s0 = mtemp1.low;
283: }
284: if(mcmp(mtemp1,c) > 0) trouble(6);
285:
286: z2 = s2;
287: z1 = s1;
288: z0 = s0;
289: mtemp1.high = s1;
290: mtemp1.low = s0;
291: return(mtemp1);
292: }
293:
294: trouble(i)
295: int i;
296: {
297: printf("mulmod: error %d\n", i);
298: abort(0);
299: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.