![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to BUG CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / ccom / vax / tests|
Return to BUG CVS log | Up to [Research Unix] / researchv10no / cmd / ccom / vax / tests |
1.1 root 1: /* foop[i] >>= 1; */
2: #include "qs.h"
3: extern char *malloc();
4: #define ison(n, v) (v[(n) >> 5] & (1 << ((n) & 31)))
5: #define on(n, v) v[(n) >> 5] |= (1 << ((n) & 31))
6: #define off(n, v) v[(n) >> 5] &= ~(1 << ((n) & 31))
7:
8: mint *one;
9:
10: vec *vecs; /* a few extra ones to make sure */
11: int nvecs, vsize;
12: unsigned int *matrix;
13: int *foop;
14:
15: int *pivot;
16:
17: svvec(n, a, b, q, orig)
18: int *a, *b;
19: mint *q, *orig;
20: { int i, p;
21: for(i = 0; i < n; i++) {
22: if((b[i] & 1) == 0)
23: continue;
24: p = a[i];
25: on(p, vecs[nvecs].rval);
26: }
27: vecs[nvecs].iam = -1;
28: vecs[nvecs].aval = itom(0);
29: move(q, vecs[nvecs].aval);
30: vecs[nvecs].qval = itom(0);
31: move(orig, vecs[nvecs].qval);
32: vecs[nvecs].ps = (int *)malloc((n + 2) * sizeof(int));
33: vecs[nvecs].exps = (short *)malloc((n + 2) * sizeof(short));
34: vecs[nvecs].vlen = n;
35: for(i = 0; i < n; i++) {
36: vecs[nvecs].ps[i] = a[i];
37: vecs[nvecs].exps[i] = b[i];
38: }
39: #ifdef DEBUG
40: checkvec(nvecs);
41: #endif
42: lu();
43: nvecs++;
44: if(nvecs >= np + 100) {
45: printf("hey, tried 100 extra times?\n");
46: exit(0);
47: }
48: }
49:
50: /* lu decomposition, maybe */
51: lu()
52: { int i;
53: /* too much work */
54: for(i = 0; i < nvecs; i++) {
55: if(vecs[i].iam < 0)
56: continue;
57: if(!ison(vecs[i].iam, vecs[nvecs].rval))
58: continue;
59: elim(vecs[i].iam, i, nvecs);
60: }
61: for(i = np - 1; i >= 0; i--)
62: if(pivot[i] == -1 && ison(i, vecs[nvecs].rval)) {
63: pivot[i] = nvecs;
64: vecs[nvecs].iam = i;
65: break;
66: }
67: #ifdef DEBUG
68: checkrow(nvecs);
69: #endif
70: if(vecs[nvecs].iam != -1)
71: return;
72: for(i = 0; i < np; i++)
73: if(ison(i, vecs[nvecs].rval) && pivot[i] == -1)
74: return;
75: doit(nvecs);
76: }
77:
78: elim(piv, a, b) /* add row a to row b */
79: { int i;
80: for(i = 0; i < vsize; i++)
81: vecs[b].rval[i] ^= vecs[a].rval[i];
82: if(ison(piv, vecs[b].rval))
83: abort("elim");
84: on(piv, vecs[b].rval);
85: }
86:
87:
88: mint *xa, *xb, *xc, *xd, *xe, *xbb;
89: doit(n) /* row n is now a winner, we hope */
90: { static flg = 0;
91: int i, j;
92: #ifdef DEBUG
93: checkrow(n);
94: #endif
95: if(flg++ == 0) {
96: xa = itom(1), xb = itom(0), xc = itom(0), xd = itom(0);
97: xe = itom(0), xbb = itom(0);
98: }
99: for(i = 0; i < np; i++)
100: foop[i] = 0;
101: move(xa, xb);
102: move(xa, xe);
103: move(xa, xbb);
104: for(i = 0; i < np; i++)
105: if(ison(i, vecs[n].rval)) {
106: refactor(pivot[i]);
107: mult(xe, vecs[pivot[i]].aval, xe);
108: mdiv(xe, num, xc, xe);
109: }
110: refactor(n);
111: mult(xe, vecs[n].aval, xe);
112: mdiv(xe, num, xc, xe);
113: move(one, xd);
114: for(i = 0; i < np; i++) {
115: if(foop[i] & 1) {
116: printf("hey, not a square foop[%d] = %d\n", i, foop[i]);
117: fflush(stdout);
118: abort();
119: }
120: foop[i] >>= 1;
121: }
122: for(i = 1; i < np; i++) /* skip -1 */
123: if(foop[i] > 0) {
124: for(j = 0; j < foop[i]; j++)
125: mult(xd, mintp[i], xd);
126: mdiv(xd, num, xc, xd);
127: }
128: msub(xe, xd, xb);
129: gcd(xb, num, xc);
130: printf("gcd(%d) ", nvecs), mout(xc);
131: if(mcmp(xc, num) != 0 && (xc->len != 1 || xc->val[0] != 1)) {
132: newfactor(xc);
133: mdiv(num, xc, xb, xe);
134: newfactor(xb);
135: return;
136: }
137: if(mcmp(xc, num) == 0)
138: return;
139: madd(xe, xd, xb);
140: mult(xe, xe, xe);
141: mult(xd, xd, xd);
142: msub(xe, xd, xbb);
143: mdiv(xbb, num, xc, xbb);
144: mdiv(xb, num, xc, xb);
145: if(xb->len != 0) {
146: for(i = 0; i < (np<100?np:100); i++)
147: if(foop[i] != 0)
148: printf("p %d e %d\n", goodp[i], foop[i]);
149: printf("...\n");
150: for(i = 0; i <= n; i++)
151: checkrow(i);
152: abort();
153: }
154: }
155:
156: refactor(n)
157: { int i;
158: if(n < 0)
159: printf("refactor(%d)?\n", n);
160: for(i = 0; i < vecs[n].vlen; i++) {
161: if(vecs[n].ps[i] >= np || vecs[n].ps[i] < 0)
162: printf("oops n %d vlen %d ps[%d] %d\n", n, vecs[n].vlen, i, vecs[n].ps[i]);
163: foop[vecs[n].ps[i]] += vecs[n].exps[i];
164: }
165: }
166:
167: #ifdef DEBUG
168: checkrow(n)
169: { vec *v = vecs + n;
170: int i, j;
171: static mint *lhs, *scratch, *rhs, *arhs;
172: static flg = 0;
173: if(flg++ == 0)
174: arhs = itom(0), rhs = itom(0), scratch = itom(0), lhs = itom(0);
175: for(i = 0; i < np; i++)
176: foop[i] = 0;
177: move(v->aval, lhs);
178: for(i = 0; i < np; i++)
179: if(ison(i, v->rval) && pivot[i] >= 0 && pivot[i] != n) {
180: mult(lhs, vecs[pivot[i]].aval, lhs);
181: mdiv(lhs, num, scratch, lhs);
182: }
183: mult(lhs, lhs, lhs);
184: mdiv(lhs, num, scratch, lhs);
185: move(v->qval, rhs);
186: refactor(n);
187: for(i = 0; i < np; i++)
188: if(ison(i, v->rval) && pivot[i] >= 0 && pivot[i] != n) {
189: mult(rhs, vecs[pivot[i]].qval, rhs);
190: mdiv(rhs, num, scratch, rhs);
191: refactor(pivot[i]);
192: }
193: if(rhs->len < 0)
194: madd(rhs, num, rhs);
195: move(one, arhs);
196: for(i = 1; i < np; i++)
197: if(foop[i] > 0) {
198: for(j = 0; j < foop[i]; j++)
199: mult(arhs, mintp[i], arhs);
200: mdiv(arhs, num, scratch, arhs);
201: }
202: if(foop[0] & 1)
203: msub(num, arhs, arhs);
204: if(mcmp(lhs, rhs) == 0 && mcmp(rhs, arhs) == 0)
205: return;
206: printf("row %d failed!\n", n);
207: printf("lhs "), mout(lhs);
208: printf("rhs "), mout(rhs);
209: printf("arhs "), mout(arhs);
210: exit(1);
211: }
212:
213: checkvec(n)
214: { static mint *xa, *xb, *xc;
215: static flg = 0;
216: int i, j;
217: vec *v = vecs + n;
218: if(!flg++) {
219: xa = itom(0), xb = itom(0), xc = itom(0);
220: }
221: mult(v->aval, v->aval, xa);
222: mdiv(xa, num, xa, xb);
223: mdiv(v->qval, num, xa, xc);
224: if(mcmp(xb, xc) != 0) {
225: printf("a * a != q mod num\n");
226: exit(1);
227: }
228: move(one, xa);
229: for(i = 0; i < v->vlen; i++) {
230: for(j = 0; j < v->exps[i]; j++) {
231: mult(xa, mintp[v->ps[i]], xa);
232: mdiv(xa, num, xb, xa);
233: }
234: }
235: move(v->qval, xb);
236: if(xb->len < 0)
237: madd(xb, num, xb);
238: if(mcmp(xa, xb) != 0) {
239: printf("vec not q! "), mout(xa), mout(v->qval);
240: exit(1);
241: }
242: }
243: #endif
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.