Return to BUG CVS log

Up to [Research Unix] / researchv10no / cmd / ccom / vax / tests

Annotation of researchv10no/cmd/ccom/vax/tests/BUG, revision 1.1

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.