[BACK]Return to gmpq.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Annotation of OpenXM_contrib2/asir2000/engine/gmpq.c, Revision 1.3

1.1       noro        1: #include "ca.h"
                      2: #include "gmp.h"
                      3: #include "base.h"
                      4: #include "inline.h"
                      5:
                      6: mpz_t ONEMPZ;
                      7: GZ ONEGZ;
                      8:
                      9: void isqrtgz(GZ a,GZ *r);
                     10: void bshiftgz(GZ a,int n,GZ *r);
                     11:
                     12: void *gc_realloc(void *p,size_t osize,size_t nsize)
                     13: {
1.2       noro       14:        return (void *)Risa_GC_realloc(p,nsize);
1.1       noro       15: }
                     16:
                     17: void gc_free(void *p,size_t size)
                     18: {
1.2       noro       19:        Risa_GC_free(p);
1.1       noro       20: }
                     21:
                     22: void init_gmpq()
                     23: {
1.2       noro       24:        mp_set_memory_functions(Risa_GC_malloc_atomic,gc_realloc,gc_free);
1.1       noro       25:
                     26:        mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOGZ(ONEMPZ,ONEGZ);
                     27: }
                     28:
                     29: GZ utogz(unsigned int u)
                     30: {
                     31:        mpz_t z;
                     32:        GZ r;
                     33:
                     34:        if ( !u ) return 0;
                     35:        mpz_init(z); mpz_set_ui(z,u); MPZTOGZ(z,r); return r;
                     36: }
                     37:
                     38: GZ stogz(int s)
                     39: {
                     40:        mpz_t z;
                     41:        GZ r;
                     42:
                     43:        if ( !s ) return 0;
                     44:        mpz_init(z); mpz_set_si(z,s); MPZTOGZ(z,r); return r;
                     45: }
                     46:
                     47: GQ mpqtogzq(mpq_t a)
                     48: {
                     49:        GZ z;
                     50:        GQ q;
                     51:
                     52:        if ( INTMPQ(a) ) {
                     53:                MPZTOGZ(mpq_numref(a),z); return (GQ)z;
                     54:        } else {
                     55:                MPQTOGQ(a,q); return q;
                     56:        }
                     57: }
                     58:
                     59: GZ ztogz(Q a)
                     60: {
                     61:        mpz_t z;
                     62:        mpq_t b;
                     63:        N nm;
                     64:        GZ s;
                     65:
                     66:        if ( !a ) return 0;
                     67:        nm = NM(a);
                     68:        mpz_init(z);
                     69:        mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
                     70:        if ( SGN(a)<0 ) mpz_neg(z,z);
                     71:        MPZTOGZ(z,s);
                     72:        return s;
                     73: }
                     74:
                     75: Q gztoz(GZ a)
                     76: {
                     77:        N nm;
                     78:        Q q;
                     79:        int sgn;
                     80:        size_t len;
                     81:
                     82:        if ( !a ) return 0;
                     83:        len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len);
                     84:        mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a));
                     85:        PL(nm) = len;
                     86:        sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q);
                     87:        return q;
                     88: }
                     89:
                     90: int n_bits_gz(GZ a)
                     91: {
                     92:        return a ? mpz_sizeinbase(BDY(a),2) : 0;
                     93: }
                     94:
                     95: GQ qtogq(Q a)
                     96: {
                     97:        mpz_t z;
                     98:        mpq_t b;
                     99:        N nm,dn;
                    100:        GZ s;
                    101:        GQ r;
                    102:
                    103:        if ( !a ) return 0;
                    104:        if ( INT(a) ) {
                    105:                nm = NM(a);
                    106:                mpz_init(z);
                    107:                mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
                    108:                if ( SGN(a)<0 ) mpz_neg(z,z);
                    109:                MPZTOGZ(z,s);
                    110:                return (GQ)s;
                    111:        } else {
                    112:                nm = NM(a); dn = DN(a);
                    113:                mpq_init(b);
                    114:                mpz_import(mpq_numref(b),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
                    115:                mpz_import(mpq_denref(b),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn));
                    116:                if ( SGN(a)<0 ) mpq_neg(b,b);
                    117:                MPQTOGQ(b,r);
                    118:                return r;
                    119:        }
                    120: }
                    121:
                    122: Q gqtoq(GQ a)
                    123: {
                    124:        N nm,dn;
                    125:        Q q;
                    126:        int sgn;
                    127:        size_t len;
                    128:
                    129:        if ( !a ) return 0;
                    130:        if ( NID(a) == N_GZ ) {
                    131:                len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len);
                    132:                mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a));
                    133:                PL(nm) = len;
                    134:                sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q);
                    135:        } else {
                    136:                len = WORDSIZE_IN_N(mpq_numref(BDY(a))); nm = NALLOC(len);
                    137:                mpz_export(BD(nm),&len,-1,sizeof(int),0,0,mpq_numref(BDY(a)));
                    138:                PL(nm) = len;
                    139:                len = WORDSIZE_IN_N(mpq_denref(BDY(a))); dn = NALLOC(len);
                    140:                mpz_export(BD(dn),&len,-1,sizeof(int),0,0,mpq_denref(BDY(a)));
                    141:                PL(dn) = len;
                    142:                sgn = mpz_sgn(mpq_numref(BDY(a))); NDTOQ(nm,dn,sgn,q);
                    143:        }
                    144:        return q;
                    145: }
                    146:
                    147: P ptogp(P a)
                    148: {
                    149:        DCP dc,dcr,dcr0;
                    150:        P b;
                    151:
                    152:        if ( !a ) return 0;
                    153:        if ( NUM(a) ) return (P)qtogq((Q)a);
                    154:        for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) {
                    155:                NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = (P)ptogp(COEF(dc));
                    156:        }
                    157:        NEXT(dcr) = 0; MKP(VR(a),dcr0,b);
                    158:        return b;
                    159: }
                    160:
                    161: P gptop(P a)
                    162: {
                    163:        DCP dc,dcr,dcr0;
                    164:        P b;
                    165:
                    166:        if ( !a ) return 0;
                    167:        if ( NUM(a) ) b = (P)gqtoq((GQ)a);
                    168:        else {
                    169:                for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) {
                    170:                        NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
                    171:                        COEF(dcr) = (P)gptop(COEF(dc));
                    172:                }
                    173:                NEXT(dcr) = 0; MKP(VR(a),dcr0,b);
                    174:        }
                    175:        return b;
                    176: }
                    177:
                    178: void addgz(GZ n1,GZ n2,GZ *nr)
                    179: {
                    180:        mpz_t t;
                    181:        int s1,s2;
                    182:
                    183:        if ( !n1 ) *nr = n2;
                    184:        else if ( !n2 ) *nr = n1;
                    185:        else {
                    186:                mpz_init(t); mpz_add(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
                    187:        }
                    188: }
                    189:
                    190: void subgz(GZ n1,GZ n2,GZ *nr)
                    191: {
                    192:        mpz_t t;
                    193:
                    194:        if ( !n1 )
                    195:                if ( !n2 )
                    196:                        *nr = 0;
                    197:                else {
                    198:                        t[0] = BDY(n2)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr);
                    199:                }
                    200:        else if ( !n2 )
                    201:                *nr = n1;
                    202:        else if ( n1 == n2 )
                    203:                *nr = 0;
                    204:        else {
                    205:                mpz_init(t); mpz_sub(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
                    206:        }
                    207: }
                    208:
                    209: void mulgz(GZ n1,GZ n2,GZ *nr)
                    210: {
                    211:        mpz_t t;
                    212:
                    213:        if ( !n1 || !n2 ) *nr = 0;
                    214:        else if ( UNIGZ(n1) ) *nr = n2;
                    215:        else if ( UNIGZ(n2) ) *nr = n1;
                    216:        else if ( MUNIGZ(n1) ) chsgngz(n2,nr);
                    217:        else if ( MUNIGZ(n2) ) chsgngz(n1,nr);
                    218:        else {
                    219:                mpz_init(t); mpz_mul(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
                    220:        }
                    221: }
                    222:
1.3     ! noro      223: void mul1gz(GZ n1,int n2,GZ *nr)
        !           224: {
        !           225:        mpz_t t;
        !           226:
        !           227:        if ( !n1 || !n2 ) *nr = 0;
        !           228:        else {
        !           229:                mpz_init(t); mpz_mul_ui(t,BDY(n1),(long)n2); MPZTOGZ(t,*nr);
        !           230:        }
        !           231: }
        !           232:
1.1       noro      233: void divgz(GZ n1,GZ n2,GZ *nq)
                    234: {
                    235:        mpz_t t;
                    236:        mpq_t a, b, q;
                    237:
                    238:        if ( !n2 ) {
                    239:                error("division by 0");
                    240:                *nq = 0;
                    241:        } else if ( !n1 )
                    242:                *nq = 0;
                    243:        else if ( n1 == n2 ) {
                    244:                mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq);
                    245:        } else {
                    246:                MPZTOMPQ(BDY(n1),a); MPZTOMPQ(BDY(n2),b);
                    247:                mpq_init(q); mpq_div(q,a,b); *nq = (GZ)mpqtogzq(q);
                    248:        }
                    249: }
                    250:
                    251: void divqrgz(GZ n1,GZ n2,GZ *nq,GZ *nr)
                    252: {
                    253:        mpz_t t, a, b, q, r;
                    254:
                    255:        if ( !n2 ) {
                    256:                error("division by 0");
                    257:                *nq = 0; *nr = 0;
                    258:        } else if ( !n1 ) {
                    259:                *nq = 0; *nr = 0;
                    260:        } else if ( n1 == n2 ) {
                    261:                mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq); *nr = 0;
                    262:        } else {
                    263:                mpz_init(q); mpz_init(r);
                    264:                mpz_fdiv_qr(q,r,BDY(n1),BDY(n2));
                    265:                if ( !mpz_sgn(q) ) *nq = 0;
                    266:                else MPZTOGZ(q,*nq);
                    267:                if ( !mpz_sgn(r) ) *nr = 0;
                    268:                else MPZTOGZ(r,*nr);
                    269:        }
                    270: }
                    271:
                    272: void divsgz(GZ n1,GZ n2,GZ *nq)
                    273: {
                    274:        mpz_t t;
                    275:        mpq_t a, b, q;
                    276:
                    277:        if ( !n2 ) {
                    278:                error("division by 0");
                    279:                *nq = 0;
                    280:        } else if ( !n1 )
                    281:                *nq = 0;
                    282:        else if ( n1 == n2 ) {
                    283:                mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq);
                    284:        } else {
                    285:                mpz_init(t); mpz_divexact(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nq);
                    286:        }
                    287: }
                    288:
                    289: void chsgngz(GZ n,GZ *nr)
                    290: {
                    291:        mpz_t t;
                    292:
                    293:        if ( !n )
                    294:                *nr = 0;
                    295:        else {
                    296:                t[0] = BDY(n)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr);
                    297:        }
                    298: }
                    299:
                    300: void pwrgz(GZ n1,Q n,GZ *nr)
                    301: {
                    302:        mpq_t t,q;
                    303:        mpz_t z;
                    304:        GQ p,r;
                    305:
                    306:        if ( !n || UNIGZ(n1) ) *nr = ONEGZ;
                    307:        else if ( !n1 ) *nr = 0;
                    308:        else if ( !INT(n) ) {
                    309:                error("can't calculate fractional power."); *nr = 0;
                    310:        } else if ( MUNIGZ(n1) ) *nr = BD(NM(n))[0]%2 ? n1 : ONEGZ;
                    311:        else if ( PL(NM(n)) > 1 ) {
                    312:                error("exponent too big."); *nr = 0;
                    313:        } else if ( NID(n1)==N_GZ && SGN(n)>0 ) {
                    314:                mpz_init(z); mpz_pow_ui(z,BDY(n1),QTOS(n)); MPZTOGZ(z,*nr);
                    315:        } else {
                    316:                MPZTOMPQ(BDY(n1),q); MPQTOGQ(q,r);
                    317:                pwrgq(r,n,&p); *nr = (GZ)p;
                    318:        }
                    319: }
                    320:
                    321: int cmpgz(GZ q1,GZ q2)
                    322: {
                    323:        int sgn;
                    324:
                    325:        if ( !q1 )
                    326:                if ( !q2 )
                    327:                        return 0;
                    328:                else
                    329:                        return -mpz_sgn(BDY(q2));
                    330:        else if ( !q2 )
                    331:                return mpz_sgn(BDY(q1));
                    332:        else if ( (sgn = mpz_sgn(BDY(q1))) != mpz_sgn(BDY(q2)) )
                    333:                        return sgn;
                    334:        else {
                    335:                sgn = mpz_cmp(BDY(q1),BDY(q2));
                    336:                if ( sgn > 0 ) return 1;
                    337:                else if ( sgn < 0 ) return -1;
                    338:                else return 0;
                    339:        }
                    340: }
                    341:
                    342: void gcdgz(GZ n1,GZ n2,GZ *nq)
                    343: {
                    344:        mpz_t t;
                    345:
                    346:        if ( !n1 ) *nq = n2;
                    347:        else if ( !n2 ) *nq = n1;
                    348:        else {
                    349:                mpz_init(t); mpz_gcd(t,BDY(n1),BDY(n2));
                    350:                MPZTOGZ(t,*nq);
                    351:        }
                    352: }
                    353:
                    354: void lcmgz(GZ n1,GZ n2,GZ *nq)
                    355: {
                    356:        GZ g,t;
                    357:
                    358:        if ( !n1 || !n2 ) *nq = 0;
                    359:        else {
                    360:                gcdgz(n1,n2,&g); divsgz(n1,g,&t);
                    361:                mulgz(n2,t,nq);
                    362:        }
                    363: }
                    364:
                    365: void gcdvgz(VECT v,GZ *q)
                    366: {
                    367:        int n,i;
                    368:        GZ *b;
                    369:        GZ g,g1;
                    370:
                    371:        n = v->len;
                    372:        b = (GZ *)v->body;
                    373:        g = b[0];
                    374:        for ( i = 1; i < n; i++ ) {
                    375:                gcdgz(g,b[i],&g1); g = g1;
                    376:        }
                    377:        *q = g;
                    378: }
                    379:
                    380: void gcdvgz_estimate(VECT v,GZ *q)
                    381: {
                    382:        int n,m,i;
                    383:        GZ s,t,u;
                    384:        GZ *b;
                    385:
                    386:        n = v->len;
                    387:        b = (GZ *)v->body;
                    388:        if ( n == 1 ) {
                    389:                if ( mpz_sgn(BDY(b[0]))<0 ) chsgngz(b[0],q);
                    390:                else *q = b[0];
                    391:        }
                    392:        m = n/2;
                    393:        for ( i = 0, s = 0; i < m; i++ ) {
                    394:                if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(s,b[i],&u);
                    395:                else addgz(s,b[i],&u);
                    396:                s = u;
                    397:        }
                    398:        for ( i = 0, t = 0; i < n; i++ ) {
                    399:                if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(t,b[i],&u);
                    400:                else addgz(t,b[i],&u);
                    401:                t = u;
                    402:        }
                    403:        gcdgz(s,t,q);
                    404: }
                    405:
                    406: void addgq(GQ n1,GQ n2,GQ *nr)
                    407: {
                    408:        mpq_t q1,q2,t;
                    409:
                    410:        if ( !n1 ) *nr = n2;
                    411:        else if ( !n2 ) *nr = n1;
                    412:        else {
                    413:                if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
                    414:                else q1[0] = BDY(n1)[0];
                    415:                if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
                    416:                else q2[0] = BDY(n2)[0];
                    417:                mpq_init(t); mpq_add(t,q1,q2); *nr = mpqtogzq(t);
                    418:        }
                    419: }
                    420:
                    421: void subgq(GQ n1,GQ n2,GQ *nr)
                    422: {
                    423:        mpq_t q1,q2,t;
                    424:
                    425:        if ( !n1 )
                    426:                if ( !n2 ) *nr = 0;
                    427:                else {
                    428:                        if ( NID(n1) == N_GZ ) chsgngz((GZ)n1,(GZ *)nr);
                    429:                        else {
                    430:                                mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOGQ(t,*nr);
                    431:                        }
                    432:                }
                    433:        else if ( !n2 ) *nr = n1;
                    434:        else if ( n1 == n2 ) *nr = 0;
                    435:        else {
                    436:                if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
                    437:                else q1[0] = BDY(n1)[0];
                    438:                if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
                    439:                else q2[0] = BDY(n2)[0];
                    440:                mpq_init(t); mpq_sub(t,q1,q2); *nr = mpqtogzq(t);
                    441:        }
                    442: }
                    443:
                    444: void mulgq(GQ n1,GQ n2,GQ *nr)
                    445: {
                    446:        mpq_t t,q1,q2;
                    447:
                    448:        if ( !n1 || !n2 ) *nr = 0;
                    449:        else {
                    450:                if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
                    451:                else q1[0] = BDY(n1)[0];
                    452:                if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
                    453:                else q2[0] = BDY(n2)[0];
                    454:                mpq_init(t); mpq_mul(t,q1,q2); *nr = mpqtogzq(t);
                    455:        }
                    456: }
                    457:
                    458: void divgq(GQ n1,GQ n2,GQ *nq)
                    459: {
                    460:        mpq_t t,q1,q2;
                    461:
                    462:        if ( !n2 ) {
                    463:                error("division by 0");
                    464:                *nq = 0;
                    465:                return;
                    466:        } else if ( !n1 ) *nq = 0;
                    467:        else if ( n1 == n2 ) *nq = (GQ)ONEGZ;
                    468:        else {
                    469:                if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
                    470:                else q1[0] = BDY(n1)[0];
                    471:                if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
                    472:                else q2[0] = BDY(n2)[0];
                    473:                mpq_init(t); mpq_div(t,q1,q2); *nq = mpqtogzq(t);
                    474:        }
                    475: }
                    476:
                    477: void chsgngq(GQ n,GQ *nr)
                    478: {
                    479:        mpq_t t;
                    480:
                    481:        if ( !n ) *nr = 0;
                    482:        else if ( NID(n) == N_GZ ) chsgngz((GZ)n,(GZ *)nr);
                    483:        else {
                    484:                mpq_init(t); mpq_neg(t,BDY(n)); MPQTOGQ(t,*nr);
                    485:        }
                    486: }
                    487:
                    488: void pwrgq(GQ n1,Q n,GQ *nr)
                    489: {
                    490:        int e;
                    491:        mpz_t nm,dn;
                    492:        mpq_t t;
                    493:
                    494:        if ( !n || UNIGZ((GZ)n1) || UNIGQ(n1) ) *nr = (GQ)ONEGZ;
                    495:        else if ( !n1 ) *nr = 0;
                    496:        else if ( !INT(n) ) {
                    497:                error("can't calculate fractional power."); *nr = 0;
                    498:        } else if ( PL(NM(n)) > 1 ) {
                    499:                error("exponent too big."); *nr = 0;
                    500:        } else {
                    501:                e = QTOS(n);
                    502:                if ( e < 0 ) {
                    503:                        e = -e;
                    504:                        if ( NID(n1)==N_GZ ) {
                    505:                                nm[0] = ONEMPZ[0];
                    506:                                dn[0] = BDY((GZ)n1)[0];
                    507:                        } else {
                    508:                                nm[0] = mpq_denref(BDY(n1))[0]; dn[0] = mpq_numref(BDY(n1))[0];
                    509:                        }
                    510:                } else {
                    511:                        if ( NID(n1)==N_GZ ) {
                    512:                                nm[0] = BDY((GZ)n1)[0]; dn[0] = ONEMPZ[0];
                    513:                        } else {
                    514:                                nm[0] = mpq_numref(BDY(n1))[0]; dn[0] = mpq_denref(BDY(n1))[0];
                    515:                        }
                    516:                }
                    517:                mpq_init(t);
                    518:                mpz_pow_ui(mpq_numref(t),nm,e); mpz_pow_ui(mpq_denref(t),dn,e);
                    519:                *nr = mpqtogzq(t);
                    520:        }
                    521: }
                    522:
                    523: int cmpgq(GQ n1,GQ n2)
                    524: {
                    525:        mpq_t q1,q2;
                    526:        int sgn;
                    527:
                    528:        if ( !n1 )
                    529:                if ( !n2 ) return 0;
                    530:                else return (NID(n2)==N_GZ) ? -mpz_sgn(BDY((GZ)n2)) : -mpq_sgn(BDY(n2));
                    531:        if ( !n2 ) return (NID(n1)==N_GZ) ? mpz_sgn(BDY((GZ)n1)) : mpq_sgn(BDY(n1));
                    532:        else if ( (sgn = mpq_sgn(BDY(n1))) != mpq_sgn(BDY(n2)) ) return sgn;
                    533:        else {
                    534:                if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
                    535:                else q1[0] = BDY(n1)[0];
                    536:                if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
                    537:                else q2[0] = BDY(n2)[0];
                    538:                sgn = mpq_cmp(q1,q2);
                    539:                if ( sgn > 0 ) return 1;
                    540:                else if ( sgn < 0 ) return -1;
                    541:                else return 0;
                    542:        }
                    543: }
                    544:
                    545: void mkgwc(int k,int l,GZ *t)
                    546: {
                    547:        mpz_t a,b,q,nm,z,u;
                    548:        int i,n;
                    549:
                    550:        n = MIN(k,l);
                    551:        mpz_init_set_ui(z,1);
                    552:        mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[0]);
                    553:        mpz_init(a); mpz_init(b); mpz_init(nm);
                    554:        for ( i = 1; i <= n; i++ ) {
                    555:                mpz_set_ui(a,k-i+1); mpz_set_ui(b,l-i+1); mpz_mul(nm,a,b);
                    556:                mpz_mul(z,BDY(t[i-1]),nm); mpz_fdiv_q_ui(z,z,i);
                    557:                mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[i]);
                    558:        }
                    559: }
                    560:
                    561: void gz_lgp(P p,GZ *g,GZ *l);
                    562:
                    563: void gz_ptozp(P p,int sgn,GQ *c,P *pr)
                    564: {
                    565:        GZ nm,dn;
                    566:
                    567:        if ( !p ) {
                    568:                *c = 0; *pr = 0;
                    569:        } else {
                    570:                gz_lgp(p,&nm,&dn);
                    571:                divgz(nm,dn,(GZ *)c);
                    572:                divsp(CO,p,(P)*c,pr);
                    573:        }
                    574: }
                    575:
                    576: void gz_lgp(P p,GZ *g,GZ *l)
                    577: {
                    578:        DCP dc;
                    579:        GZ g1,g2,l1,l2,l3,l4;
                    580:
                    581:        if ( NUM(p) ) {
                    582:                if ( NID((GZ)p)==N_GZ ) {
                    583:                        MPZTOGZ(BDY((GZ)p),*g);
                    584:                        *l = ONEGZ;
                    585:                } else {
                    586:                        MPZTOGZ(mpq_numref(BDY((GQ)p)),*g);
                    587:                        MPZTOGZ(mpq_denref(BDY((GQ)p)),*l);
                    588:                }
                    589:        } else {
                    590:                dc = DC(p); gz_lgp(COEF(dc),g,l);
                    591:                for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
                    592:                        gz_lgp(COEF(dc),&g1,&l1); gcdgz(*g,g1,&g2); *g = g2;
                    593:                        gcdgz(*l,l1,&l2); mulgz(*l,l1,&l3); divgz(l3,l2,l);
                    594:                }
                    595:        }
                    596: }
                    597:
                    598: void gz_qltozl(GQ *w,int n,GZ *dvr)
                    599: {
                    600:        GZ nm,dn;
                    601:        GZ g,g1,l1,l2,l3;
                    602:        GQ c;
                    603:        int i;
                    604:        struct oVECT v;
                    605:
                    606:        for ( i = 0; i < n; i++ )
                    607:                if ( w[i] && NID(w[i])==N_GQ )
                    608:                        break;
                    609:        if ( i == n ) {
                    610:                v.id = O_VECT; v.len = n; v.body = (pointer *)w;
                    611:                gcdvgz(&v,dvr); return;
                    612:        }
                    613:        for ( i = 0; !w[i]; i++ );
                    614:        c = w[i];
                    615:        if ( NID(c)==N_GQ ) {
                    616:                MPZTOGZ(mpq_numref(BDY(c)),nm); MPZTOGZ(mpq_denref(BDY(c)),dn);
                    617:        } else {
                    618:                MPZTOGZ(BDY((GZ)c),nm); dn = ONEGZ;
                    619:        }
                    620:        for ( i++; i < n; i++ ) {
                    621:                c = w[i];
                    622:                if ( !c ) continue;
                    623:                if ( NID(c)==N_GQ ) {
                    624:                        MPZTOGZ(mpq_numref(BDY(c)),g1); MPZTOGZ(mpq_denref(BDY(c)),l1);
                    625:                } else {
                    626:                        MPZTOGZ(BDY((GZ)c),g1); l1 = ONEGZ;
                    627:                }
                    628:                gcdgz(nm,g1,&g); nm = g;
                    629:                gcdgz(dn,l1,&l2); mulgz(dn,l1,&l3); divgz(l3,l2,&dn);
                    630:        }
                    631:        divgz(nm,dn,dvr);
                    632: }
                    633:
                    634: int gz_bits(GQ q)
                    635: {
                    636:        if ( !q ) return 0;
                    637:        else if ( NID(q)==N_Q )
                    638:                return n_bits(NM((Q)q))+(INT((Q)q)?0:n_bits(DN((Q)q)));
                    639:        else if ( NID(q)==N_GZ ) return mpz_sizeinbase(BDY((GZ)q),2);
                    640:        else
                    641:                return mpz_sizeinbase(mpq_numref(BDY(q)),2)
                    642:                        + mpz_sizeinbase(mpq_denref(BDY(q)),2);
                    643: }
                    644:
                    645: int gzp_mag(P p)
                    646: {
                    647:        int s;
                    648:        DCP dc;
                    649:
                    650:        if ( !p ) return 0;
                    651:        else if ( OID(p) == O_N ) return gz_bits((GQ)p);
                    652:        else {
                    653:                for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) s += gzp_mag(COEF(dc));
                    654:                return s;
                    655:        }
                    656: }
                    657:
                    658: void makesubstgz(VL v,NODE *s)
                    659: {
                    660:        NODE r,r0;
                    661:        GZ q;
                    662:        unsigned int n;
                    663:
                    664:        for ( r0 = 0; v; v = NEXT(v) ) {
                    665:                NEXTNODE(r0,r); BDY(r) = (pointer)v->v;
                    666: #if defined(_PA_RISC1_1)
                    667:                n = mrand48()&BMASK; q = utogz(n);
                    668: #else
                    669:                n = random(); q = utogz(n);
                    670: #endif
                    671:                NEXTNODE(r0,r); BDY(r) = (pointer)q;
                    672:        }
                    673:        if ( r0 ) NEXT(r) = 0;
                    674:        *s = r0;
                    675: }
                    676:
                    677: unsigned int remgq(GQ a,unsigned int mod)
                    678: {
                    679:        unsigned int c,nm,dn;
                    680:        mpz_t r;
                    681:
                    682:        if ( !a ) return 0;
                    683:        else if ( NID(a)==N_GZ ) {
                    684:                mpz_init(r);
                    685:                c = mpz_fdiv_r_ui(r,BDY((GZ)a),mod);
                    686:        } else {
                    687:                mpz_init(r);
                    688:                nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod);
                    689:                dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod);
                    690:                dn = invm(dn,mod);
                    691:                DMAR(nm,dn,0,mod,c);
                    692:        }
                    693:        return c;
                    694: }
                    695:
                    696: extern int DP_Print;
                    697:
                    698: #define GZ_F4_INTRAT_PERIOD 8
                    699:
                    700: int gz_generic_gauss_elim(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
                    701: {
                    702:        int **wmat;
                    703:        GZ **bmat,**tmat,*bmi,*tmi;
                    704:        GZ q,m1,m2,m3,s,u;
                    705:        int *wmi,*colstat,*wcolstat,*rind,*cind;
                    706:        int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
                    707:        MAT r,crmat;
                    708:        int ret;
                    709:
                    710:        bmat = (GZ **)mat->body;
                    711:        row = mat->row; col = mat->col;
                    712:        wmat = (int **)almat(row,col);
                    713:        colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                    714:        wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                    715:        for ( ind = 0; ; ind++ ) {
                    716:                if ( DP_Print ) {
                    717:                        fprintf(asir_out,"."); fflush(asir_out);
                    718:                }
                    719:                md = get_lprime(ind);
                    720:                for ( i = 0; i < row; i++ )
                    721:                        for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
                    722:                                wmi[j] = remgq((GQ)bmi[j],md);
                    723:                rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
                    724:                if ( !ind ) {
                    725: RESET:
                    726:                        m1 = utogz(md);
                    727:                        rank0 = rank;
                    728:                        bcopy(wcolstat,colstat,col*sizeof(int));
                    729:                        MKMAT(crmat,rank,col-rank);
                    730:                        MKMAT(r,rank,col-rank); *nm = r;
                    731:                        tmat = (GZ **)crmat->body;
                    732:                        for ( i = 0; i < rank; i++ )
                    733:                                for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                    734:                                        if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]);
                    735:                } else {
                    736:                        if ( rank < rank0 ) {
                    737:                                if ( DP_Print ) {
                    738:                                        fprintf(asir_out,"lower rank matrix; continuing...\n");
                    739:                                        fflush(asir_out);
                    740:                                }
                    741:                                continue;
                    742:                        } else if ( rank > rank0 ) {
                    743:                                if ( DP_Print ) {
                    744:                                        fprintf(asir_out,"higher rank matrix; resetting...\n");
                    745:                                        fflush(asir_out);
                    746:                                }
                    747:                                goto RESET;
                    748:                        } else {
                    749:                                for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
                    750:                                if ( j < col ) {
                    751:                                        if ( DP_Print ) {
                    752:                                                fprintf(asir_out,"inconsitent colstat; resetting...\n");
                    753:                                                fflush(asir_out);
                    754:                                        }
                    755:                                        goto RESET;
                    756:                                }
                    757:                        }
                    758:
                    759:                        inv = invm(remgq((GQ)m1,md),md);
                    760:                        m2 = utogz(md); mulgz(m1,m2,&m3);
                    761:                        for ( i = 0; i < rank; i++ )
                    762:                                for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                    763:                                        if ( !colstat[j] ) {
                    764:                                                if ( tmi[k] ) {
                    765:                                                /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                    766:                                                        t = remgq((GQ)tmi[k],md);
                    767:                                                        if ( wmi[j] >= t ) t = wmi[j]-t;
                    768:                                                        else t = md-(t-wmi[j]);
                    769:                                                        DMAR(t,inv,0,md,t1)
                    770:                                                        u = utogz(t1); mulgz(m1,u,&s);
                    771:                                                        addgz(tmi[k],s,&u); tmi[k] = u;
                    772:                                                } else if ( wmi[j] ) {
                    773:                                                /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                    774:                                                        DMAR(wmi[j],inv,0,md,t)
                    775:                                                        u = utogz(t); mulgz(m1,u,&s); tmi[k] = s;
                    776:                                                }
                    777:                                                k++;
                    778:                                        }
                    779:                        m1 = m3;
                    780:                        if ( ind % GZ_F4_INTRAT_PERIOD )
                    781:                                ret = 0;
                    782:                        else
                    783:                                ret = gz_intmtoratm(crmat,m1,*nm,dn);
                    784:                        if ( ret ) {
                    785:                                *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                    786:                                *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
                    787:                                for ( j = k = l = 0; j < col; j++ )
                    788:                                        if ( colstat[j] ) rind[k++] = j;
                    789:                                        else cind[l++] = j;
                    790:                                if ( gz_gensolve_check(mat,*nm,*dn,rind,cind) )
                    791:                                        return rank;
                    792:                        }
                    793:                }
                    794:        }
                    795: }
                    796:
                    797: int gz_generic_gauss_elim2(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
                    798: {
                    799:
                    800:        MAT full;
                    801:        GZ **bmat,**b;
                    802:        GZ *bmi;
                    803:        GZ dn0;
                    804:        int row,col,md,i,j,rank,ret;
                    805:        int **wmat;
                    806:        int *wmi;
                    807:        int *colstat,*rowstat;
                    808:
                    809:        bmat = (GZ **)mat->body;
                    810:        row = mat->row; col = mat->col;
                    811:        wmat = (int **)almat(row,col);
                    812:        colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                    813:        rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
                    814:        /* XXX */
                    815:        md = get_lprime(0);
                    816:        for ( i = 0; i < row; i++ )
                    817:                for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
                    818:                        wmi[j] = remgq((GQ)bmi[j],md);
                    819:        rank = generic_gauss_elim_mod2(wmat,row,col,md,colstat,rowstat);
                    820:        b = (GZ **)MALLOC(rank*sizeof(GZ));
                    821:        for ( i = 0; i < rank; i++ ) b[i] = bmat[rowstat[i]];
                    822:        NEWMAT(full); full->row = rank; full->col = col; full->body = (pointer **)b;
                    823:        ret = gz_generic_gauss_elim_full(full,nm,dn,rindp,cindp);
                    824:        if ( !ret ) {
                    825:                rank = gz_generic_gauss_elim(mat,nm,&dn0,rindp,cindp);
                    826:                for ( i = 0; i < rank; i++ ) dn[i] = dn0;
                    827:        }
                    828:        return rank;
                    829: }
                    830:
                    831: int gz_generic_gauss_elim_full(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
                    832: {
                    833:        int **wmat;
                    834:        int *stat;
                    835:        GZ **bmat,**tmat,*bmi,*tmi,*ri;
                    836:        GZ q,m1,m2,m3,s,u;
                    837:        int *wmi,*colstat,*wcolstat,*rind,*cind;
                    838:        int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv,h;
                    839:        MAT r,crmat;
                    840:        int ret,initialized,done;
                    841:
                    842:        initialized = 0;
                    843:        bmat = (GZ **)mat->body;
                    844:        row = mat->row; col = mat->col;
                    845:        wmat = (int **)almat(row,col);
                    846:        stat = (int *)MALLOC_ATOMIC(row*sizeof(int));
                    847:        for ( i = 0; i < row; i++ ) stat[i] = 0;
                    848:        colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                    849:        wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                    850:        for ( ind = 0; ; ind++ ) {
                    851:                if ( DP_Print ) {
                    852:                        fprintf(asir_out,"."); fflush(asir_out);
                    853:                }
                    854:                md = get_lprime(ind);
                    855:                for ( i = 0; i < row; i++ )
                    856:                        for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
                    857:                                wmi[j] = remgq((GQ)bmi[j],md);
                    858:                rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
                    859:                if ( rank < row ) continue;
                    860:                if ( !initialized ) {
                    861:                        m1 = utogz(md);
                    862:                        bcopy(wcolstat,colstat,col*sizeof(int));
                    863:                        MKMAT(crmat,row,col-row);
                    864:                        MKMAT(r,row,col-row); *nm = r;
                    865:                        tmat = (GZ **)crmat->body;
                    866:                        for ( i = 0; i < row; i++ )
                    867:                                for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                    868:                                        if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]);
                    869:                        initialized = 1;
                    870:                } else {
                    871:                        for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
                    872:                        if ( j < col ) continue;
                    873:
                    874:                        inv = invm(remgq((GQ)m1,md),md);
                    875:                        m2 = utogz(md); mulgz(m1,m2,&m3);
                    876:                        for ( i = 0; i < row; i++ )
                    877:                                switch ( stat[i] ) {
                    878:                                        case 1:
                    879:                                                /* consistency check */
                    880:                                                ri = (GZ *)BDY(r)[i]; wmi = wmat[i];
                    881:                                                for ( j = 0; j < col; j++ ) if ( colstat[j] ) break;
                    882:                                                h = md-remgq((GQ)dn[i],md);
                    883:                                                for ( j++, k = 0; j < col; j++ )
                    884:                                                        if ( !colstat[j] ) {
                    885:                                                                t = remgq((GQ)ri[k],md);
                    886:                                                                DMAR(wmi[i],h,t,md,t1);
                    887:                                                                if ( t1 ) break;
                    888:                                                        }
                    889:                                                if ( j == col ) { stat[i]++; break; }
                    890:                                                else {
                    891:                                                        /* fall to the case 0 */
                    892:                                                        stat[i] = 0;
                    893:                                                }
                    894:                                        case 0:
                    895:                                                tmi = tmat[i]; wmi = wmat[i];
                    896:                                                for ( j = k = 0; j < col; j++ )
                    897:                                                        if ( !colstat[j] ) {
                    898:                                                                if ( tmi[k] ) {
                    899:                                                                /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                    900:                                                                        t = remgq((GQ)tmi[k],md);
                    901:                                                                        if ( wmi[j] >= t ) t = wmi[j]-t;
                    902:                                                                        else t = md-(t-wmi[j]);
                    903:                                                                        DMAR(t,inv,0,md,t1)
                    904:                                                                        u = utogz(t1); mulgz(m1,u,&s);
                    905:                                                                        addgz(tmi[k],s,&u); tmi[k] = u;
                    906:                                                                } else if ( wmi[j] ) {
                    907:                                                                /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                    908:                                                                        DMAR(wmi[j],inv,0,md,t)
                    909:                                                                        u = utogz(t); mulgz(m1,u,&s); tmi[k] = s;
                    910:                                                                }
                    911:                                                                k++;
                    912:                                                        }
                    913:                                                break;
                    914:                                        case 2: default:
                    915:                                                break;
                    916:                                }
                    917:                        m1 = m3;
                    918:                        if ( ind % 4 )
                    919:                                ret = 0;
                    920:                        else
                    921:                                ret = gz_intmtoratm2(crmat,m1,*nm,dn,stat);
                    922:                        if ( ret ) {
                    923:                                *rindp = rind = (int *)MALLOC_ATOMIC(row*sizeof(int));
                    924:                                *cindp = cind = (int *)MALLOC_ATOMIC((col-row)*sizeof(int));
                    925:                                for ( j = k = l = 0; j < col; j++ )
                    926:                                        if ( colstat[j] ) rind[k++] = j;
                    927:                                        else cind[l++] = j;
                    928:                                return gz_gensolve_check2(mat,*nm,dn,rind,cind);
                    929:                        }
                    930:                }
                    931:        }
                    932: }
                    933:
                    934: int gz_generic_gauss_elim_direct(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp){
                    935:        GZ **bmat,*s;
                    936:        GZ u,v,w,x,d,t,y;
                    937:        int row,col,i,j,k,l,m,rank;
                    938:        int *colstat,*colpos,*cind;
                    939:        MAT r,in;
                    940:
                    941:        row = mat->row; col = mat->col;
                    942:        MKMAT(in,row,col);
                    943:        for ( i = 0; i < row; i++ )
                    944:                for ( j = 0; j < col; j++ ) in->body[i][j] = mat->body[i][j];
                    945:        bmat = (GZ **)in->body;
                    946:        colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                    947:        *rindp = colpos = (int *)MALLOC_ATOMIC(row*sizeof(int));
                    948:        for ( j = 0, rank = 0, d = ONEGZ; j < col; j++ ) {
                    949:                for ( i = rank; i < row && !bmat[i][j]; i++  );
                    950:                if ( i == row ) { colstat[j] = 0; continue; }
                    951:                else { colstat[j] = 1; colpos[rank] = j; }
                    952:                if ( i != rank )
                    953:                        for ( k = j; k < col; k++ ) {
                    954:                                t = bmat[i][k]; bmat[i][k] = bmat[rank][k]; bmat[rank][k] = t;
                    955:                        }
                    956:                for ( i = rank+1, v = bmat[rank][j]; i < row; i++ )
                    957:                        for ( k = j, u = bmat[i][j]; k < col; k++ ) {
                    958:                                mulgz(bmat[i][k],v,&w); mulgz(bmat[rank][k],u,&x);
                    959:                                subgz(w,x,&y); divsgz(y,d,&bmat[i][k]);
                    960:                        }
                    961:                d = v; rank++;
                    962:        }
                    963:        *dn = d;
                    964:        s = (GZ *)MALLOC(col*sizeof(GZ));
                    965:        for ( i = rank-1; i >= 0; i-- ) {
                    966:                for ( k = colpos[i]; k < col; k++ ) mulgz(bmat[i][k],d,&s[k]);
                    967:                for ( m = rank-1; m > i; m-- ) {
                    968:                        for ( k = colpos[m], u = bmat[i][k]; k < col; k++ ) {
                    969:                                mulgz(bmat[m][k],u,&w); subgz(s[k],w,&x); s[k] = x;
                    970:                        }
                    971:                }
                    972:                for ( k = colpos[i], u = bmat[i][k]; k < col; k++ )
                    973:                        divgz(s[k],u,&bmat[i][k]);
                    974:        }
                    975:        *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
                    976:        MKMAT(r,rank,col-rank); *nm = r;
                    977:        for ( j = 0, k = 0; j < col; j++ )
                    978:                if ( !colstat[j] ) {
                    979:                        cind[k] = j;
                    980:                        for ( i = 0; i < rank; i++ ) r->body[i][k] = bmat[i][j];
                    981:                        k++;
                    982:                }
                    983:        return rank;
                    984: }
                    985:
                    986: int gz_intmtoratm(MAT mat,GZ md,MAT nm,GZ *dn)
                    987: {
                    988:        GZ t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy;
                    989:        int i,j,k,l,row,col,sgn,ret;
                    990:        GZ **rmat,**tmat,*tmi,*nmk;
                    991:
                    992:        if ( UNIGZ(md) )
                    993:                return 0;
                    994:        row = mat->row; col = mat->col;
                    995:        bshiftgz(md,1,&t);
                    996:        isqrtgz(t,&s);
                    997:        bshiftgz(s,64,&b);
                    998:        if ( !b ) b = ONEGZ;
                    999:        dn0 = ONEGZ;
                   1000:        tmat = (GZ **)mat->body;
                   1001:        rmat = (GZ **)nm->body;
                   1002:        for ( i = 0; i < row; i++ )
                   1003:                for ( j = 0, tmi = tmat[i]; j < col; j++ )
                   1004:                        if ( tmi[j] ) {
                   1005:                                mulgz(tmi[j],dn0,&s);
                   1006:                                divqrgz(s,md,&dmy,&u);
                   1007:                                ret = gz_inttorat(u,md,b,&sgn,&unm,&udn);
                   1008:                                if ( !ret ) return 0;
                   1009:                                else {
                   1010:                                        if ( sgn < 0 ) chsgngz(unm,&nm1);
                   1011:                                        else nm1 = unm;
                   1012:                                        dn1 = udn;
                   1013:                                        if ( !UNIGZ(dn1) ) {
                   1014:                                                for ( k = 0; k < i; k++ )
                   1015:                                                        for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
                   1016:                                                                mulgz(nmk[l],dn1,&q); nmk[l] = q;
                   1017:                                                        }
                   1018:                                                for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
                   1019:                                                        mulgz(nmk[l],dn1,&q); nmk[l] = q;
                   1020:                                                }
                   1021:                                        }
                   1022:                                        rmat[i][j] = nm1;
                   1023:                                        mulgz(dn0,dn1,&q); dn0 = q;
                   1024:                                }
                   1025:                        }
                   1026:        *dn = dn0;
                   1027:        return 1;
                   1028: }
                   1029:
                   1030: int gz_intmtoratm2(MAT mat,GZ md,MAT nm,GZ *dn,int *stat)
                   1031: {
                   1032:        int row,col,i,j,ret;
                   1033:        GZ dn0,dn1,t,s,b;
                   1034:        GZ *w,*tmi;
                   1035:        GZ **tmat;
                   1036:
                   1037:        bshiftgz(md,1,&t);
                   1038:        isqrtgz(t,&s);
                   1039:        bshiftgz(s,64,&b);
                   1040:        tmat = (GZ **)mat->body;
                   1041:        if ( UNIGZ(md) ) return 0;
                   1042:        row = mat->row; col = mat->col;
                   1043:        dn0 = ONEGZ;
                   1044:        for ( i = 0; i < row; i++ )
                   1045:                if ( cmpgz(dn[i],dn0) > 0 ) dn0 = dn[i];
                   1046:        w = (GZ *)MALLOC(col*sizeof(GZ));
                   1047:        for ( i = 0; i < row; i++ )
                   1048:                if ( stat[i] == 0 ) {
                   1049:                        for ( j = 0, tmi = tmat[i]; j < col; j++ )
                   1050:                                        mulgz(tmi[j],dn0,&w[j]);
                   1051:                        ret = gz_intvtoratv(w,col,md,b,BDY(nm)[i],&dn[i]);
                   1052:                        if ( ret ) {
                   1053:                                stat[i] = 1;
                   1054:                                mulgz(dn0,dn[i],&t); dn[i] = t; dn0 = t;
                   1055:                        }
                   1056:                }
                   1057:        for ( i = 0; i < row; i++ ) if ( !stat[i] ) break;
                   1058:        if ( i == row ) return 1;
                   1059:        else return 0;
                   1060: }
                   1061:
                   1062: int gz_intvtoratv(GZ *v,int n,GZ md,GZ b,GZ *nm,GZ *dn)
                   1063: {
                   1064:        GZ dn0,dn1,q,s,u,nm1,unm,udn,dmy;
                   1065:        GZ *nmk;
                   1066:        int j,l,col,ret,sgn;
                   1067:
                   1068:        for ( j = 0; j < n; j++ ) nm[j] = 0;
                   1069:        dn0 = ONEGZ;
                   1070:        for ( j = 0; j < n; j++ ) {
                   1071:                if ( !v[j] ) continue;
                   1072:                mulgz(v[j],dn0,&s);
                   1073:                divqrgz(s,md,&dmy,&u);
                   1074:                ret = gz_inttorat(u,md,b,&sgn,&unm,&udn);
                   1075:                if ( !ret ) return 0;
                   1076:                if ( sgn < 0 ) chsgngz(unm,&nm1);
                   1077:                else nm1 = unm;
                   1078:                dn1 = udn;
                   1079:                if ( !UNIGZ(dn1) )
                   1080:                        for ( l = 0; l < j; l++ ) {
                   1081:                                mulgz(nm[l],dn1,&q); nm[l] = q;
                   1082:                        }
                   1083:                nm[j] = nm1;
                   1084:                mulgz(dn0,dn1,&q); dn0 = q;
                   1085:        }
                   1086:        *dn = dn0;
                   1087:        return 1;
                   1088: }
                   1089:
                   1090: /* assuming 0 < c < m */
                   1091:
                   1092: int gz_inttorat(GZ c,GZ m,GZ b,int *sgnp,GZ *nmp,GZ *dnp)
                   1093: {
                   1094:        GZ qq,t,u1,v1,r1;
                   1095:        GZ q,u2,v2,r2;
                   1096:
                   1097:        u1 = 0; v1 = ONEGZ; u2 = m; v2 = c;
                   1098:        while ( cmpgz(v2,b) >= 0 ) {
                   1099:                divqrgz(u2,v2,&q,&r2); u2 = v2; v2 = r2;
                   1100:                mulgz(q,v1,&t); subgz(u1,t,&r1); u1 = v1; v1 = r1;
                   1101:        }
                   1102:        if ( cmpgz(v1,b) >= 0 ) return 0;
                   1103:        else {
                   1104:                *nmp = v2;
                   1105:                if ( mpz_sgn(BDY(v1))<0  ) {
                   1106:                        *sgnp = -1; chsgngz(v1,dnp);
                   1107:                } else {
                   1108:                        *sgnp = 1; *dnp = v1;
                   1109:                }
                   1110:                return 1;
                   1111:        }
                   1112: }
                   1113:
                   1114: extern int f4_nocheck;
                   1115:
                   1116: int gz_gensolve_check(MAT mat,MAT nm,GZ dn,int *rind,int *cind)
                   1117: {
                   1118:        int row,col,rank,clen,i,j,k,l;
                   1119:        GZ s,t;
                   1120:        GZ *w;
                   1121:        GZ *mati,*nmk;
                   1122:
                   1123:        if ( f4_nocheck ) return 1;
                   1124:        row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
                   1125:        w = (GZ *)MALLOC(clen*sizeof(GZ));
                   1126:        for ( i = 0; i < row; i++ ) {
                   1127:                mati = (GZ *)mat->body[i];
                   1128:                bzero(w,clen*sizeof(GZ));
                   1129:                for ( k = 0; k < rank; k++ )
                   1130:                        for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) {
                   1131:                                mulgz(mati[rind[k]],nmk[l],&t); addgz(w[l],t,&s); w[l] = s;
                   1132:                        }
                   1133:                for ( j = 0; j < clen; j++ ) {
                   1134:                        mulgz(dn,mati[cind[j]],&t);
                   1135:                        if ( cmpgz(w[j],t) ) break;
                   1136:                }
                   1137:                if ( j != clen ) break;
                   1138:        }
                   1139:        if ( i != row ) return 0;
                   1140:        else return 1;
                   1141: }
                   1142:
                   1143: int gz_gensolve_check2(MAT mat,MAT nm,GZ *dn,int *rind,int *cind)
                   1144: {
                   1145:        int row,col,rank,clen,i,j,k,l;
                   1146:        GZ s,t,u,d;
                   1147:        GZ *w,*m;
                   1148:        GZ *mati,*nmk;
                   1149:
                   1150:        if ( f4_nocheck ) return 1;
                   1151:        row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
                   1152:        w = (GZ *)MALLOC(clen*sizeof(GZ));
                   1153:        m = (GZ *)MALLOC(clen*sizeof(GZ));
                   1154:        for ( d = dn[0], i = 1; i < rank; i++ ) {
                   1155:                lcmgz(d,dn[i],&t); d = t;
                   1156:        }
                   1157:        for ( i = 0; i < rank; i++ ) divsgz(d,dn[i],&m[i]);
                   1158:        for ( i = 0; i < row; i++ ) {
                   1159:                mati = (GZ *)mat->body[i];
                   1160:                bzero(w,clen*sizeof(GZ));
                   1161:                for ( k = 0; k < rank; k++ ) {
                   1162:                        mulgz(mati[rind[k]],m[k],&u);
                   1163:                        for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) {
                   1164:                                mulgz(u,nmk[l],&t); addgz(w[l],t,&s); w[l] = s;
                   1165:                        }
                   1166:                }
                   1167:                for ( j = 0; j < clen; j++ ) {
                   1168:                        mulgz(d,mati[cind[j]],&t);
                   1169:                        if ( cmpgz(w[j],t) ) break;
                   1170:                }
                   1171:                if ( j != clen ) break;
                   1172:        }
                   1173:        if ( i != row ) return 0;
                   1174:        else return 1;
                   1175: }
                   1176:
                   1177: void isqrtgz(GZ a,GZ *r)
                   1178: {
                   1179:        int k;
                   1180:        GZ x,t,x2,xh,quo,rem;
                   1181:        Q two;
                   1182:
                   1183:        if ( !a ) *r = 0;
                   1184:        else if ( UNIGZ(a) ) *r = ONEGZ;
                   1185:        else {
                   1186:                k = gz_bits((GQ)a); /* a <= 2^k-1 */
                   1187:                bshiftgz(ONEGZ,-((k>>1)+(k&1)),&x); /* a <= x^2 */
                   1188:                STOQ(2,two);
                   1189:                while ( 1 ) {
                   1190:                        pwrgz(x,two,&t);
                   1191:                        if ( cmpgz(t,a) <= 0 ) {
                   1192:                                *r = x; return;
                   1193:                        } else {
                   1194:                                if ( mpz_tstbit(BDY(x),0) ) addgz(x,a,&t);
                   1195:                                else t = a;
                   1196:                                bshiftgz(x,-1,&x2); divqrgz(t,x2,&quo,&rem);
                   1197:                                bshiftgz(x,1,&xh); addgz(quo,xh,&x);
                   1198:                        }
                   1199:                }
                   1200:        }
                   1201: }
                   1202:
                   1203: void bshiftgz(GZ a,int n,GZ *r)
                   1204: {
                   1205:        mpz_t t;
                   1206:
                   1207:        if ( !a ) *r = 0;
                   1208:        else if ( n == 0 ) *r = a;
                   1209:        else if ( n < 0 ) {
                   1210:                mpz_init(t); mpz_mul_2exp(t,BDY(a),-n); MPZTOGZ(t,*r);
                   1211:        } else {
                   1212:                mpz_init(t); mpz_fdiv_q_2exp(t,BDY(a),n);
                   1213:                if ( !mpz_sgn(t) ) *r = 0;
                   1214:                else MPZTOGZ(t,*r);
                   1215:        }
                   1216: }
                   1217:

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>