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

Annotation of OpenXM_contrib2/asir2018/engine/Q.c, Revision 1.7

1.7     ! noro        1: /* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.6 2018/10/01 05:49:06 noro Exp $ */
1.1       noro        2: #include "ca.h"
                      3: #include "gmp.h"
                      4: #include "base.h"
                      5: #include "inline.h"
                      6:
                      7: mpz_t ONEMPZ;
                      8: Z ONE;
                      9: int lf_lazy;
                     10: Z current_mod_lf;
                     11: int current_mod_lf_size;
                     12: gmp_randstate_t GMP_RAND;
                     13:
1.6       noro       14: #define F4_INTRAT_PERIOD 8
                     15:
                     16: extern int DP_Print;
                     17:
1.1       noro       18: void isqrtz(Z a,Z *r);
                     19: void bshiftz(Z a,int n,Z *r);
                     20:
                     21: void *gc_realloc(void *p,size_t osize,size_t nsize)
                     22: {
                     23:   return (void *)Risa_GC_realloc(p,nsize);
                     24: }
                     25:
                     26: void gc_free(void *p,size_t size)
                     27: {
                     28:   Risa_GC_free(p);
                     29: }
                     30:
                     31: void init_gmpq()
                     32: {
1.3       noro       33:   mp_set_memory_functions(Risa_GC_malloc_atomic,gc_realloc,gc_free);
1.1       noro       34:
                     35:   mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOZ(ONEMPZ,ONE);
                     36:   gmp_randinit_default(GMP_RAND);
                     37: }
                     38:
1.7     ! noro       39: void printexpr(VL,Obj);
        !            40:
1.3       noro       41: void pmat(Z **a,int row,int col)
                     42: {
                     43:   int i,j;
                     44:
                     45:   for ( i = 0; i < row; i++, printf("\n") )
                     46:     for ( j = 0; j < col; j++, printf(" ") )
1.7     ! noro       47:       printexpr(CO,(Obj)a[i][j]);
1.3       noro       48:   printf("\n");
                     49: }
                     50:
1.1       noro       51: Z utoz(unsigned int u)
                     52: {
                     53:   mpz_t z;
                     54:   Z r;
                     55:
                     56:   if ( !u ) return 0;
                     57:   mpz_init(z); mpz_set_ui(z,u); MPZTOZ(z,r); return r;
                     58: }
                     59:
                     60: Z stoz(int s)
                     61: {
                     62:   mpz_t z;
                     63:   Z r;
                     64:
                     65:   if ( !s ) return 0;
                     66:   mpz_init(z); mpz_set_si(z,s); MPZTOZ(z,r); return r;
                     67: }
                     68:
                     69: int sgnz(Z z)
                     70: {
                     71:   if ( !z ) return 0;
                     72:   else return mpz_sgn(BDY(z));
                     73: }
                     74:
                     75: void nmq(Q q,Z *r)
                     76: {
                     77:   if ( !q ) *r = 0;
                     78:   else if ( INT(q) ) *r = (Z)q;
                     79:   else {
                     80:     MPZTOZ(mpq_numref(BDY(q)),*r);
                     81:   }
                     82: }
                     83:
                     84: void dnq(Q q,Z *r)
                     85: {
                     86:   if ( !q ) *r = 0;
                     87:   else if ( INT(q) ) *r = ONE;
                     88:   else {
                     89:     MPZTOZ(mpq_denref(BDY(q)),*r);
                     90:   }
                     91: }
                     92:
                     93: int sgnq(Q q)
                     94: {
                     95:   if ( !q ) return 0;
                     96:   else if ( q->z ) return mpz_sgn(BDY((Z)q));
                     97:   else return mpz_sgn(mpq_numref(BDY(q)));
                     98: }
                     99:
                    100: Q mpqtozq(mpq_t a)
                    101: {
                    102:   Z z;
                    103:   Q q;
                    104:
                    105:   if ( INTMPQ(a) ) {
                    106:     MPZTOZ(mpq_numref(a),z); return (Q)z;
                    107:   } else {
                    108:     MPQTOQ(a,q); return q;
                    109:   }
                    110: }
                    111:
                    112: void dupz(Z a,Z *b)
                    113: {
                    114:   mpz_t t;
                    115:
                    116:   if ( !a ) *b = a;
                    117:   else {
                    118:     mpz_init(t); mpz_set(t,BDY(a)); MPZTOZ(t,*b);
                    119:   }
                    120: }
                    121:
                    122: int n_bits_z(Z a)
                    123: {
                    124:   return a ? mpz_sizeinbase(BDY(a),2) : 0;
                    125: }
                    126:
                    127: void addz(Z n1,Z n2,Z *nr)
                    128: {
                    129:   mpz_t t;
                    130:   int s1,s2;
                    131:
                    132:   if ( !n1 ) *nr = n2;
                    133:   else if ( !n2 ) *nr = n1;
                    134:   else if ( !n1->z || !n2->z )
                    135:     error("addz : invalid argument");
                    136:   else {
                    137:     mpz_init(t); mpz_add(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nr);
                    138:   }
                    139: }
                    140:
                    141: void subz(Z n1,Z n2,Z *nr)
                    142: {
                    143:   mpz_t t;
                    144:
                    145:   if ( !n1 ) {
                    146:     if ( !n2 )
                    147:       *nr = 0;
                    148:     else
                    149:       chsgnz(n2,nr);
                    150:   } else if ( !n2 )
                    151:     *nr = n1;
                    152:   else if ( n1 == n2 )
                    153:     *nr = 0;
                    154:   else if ( !n1->z || !n2->z )
                    155:     error("subz : invalid argument");
                    156:   else {
                    157:     mpz_init(t); mpz_sub(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nr);
                    158:   }
                    159: }
                    160:
                    161: void mulz(Z n1,Z n2,Z *nr)
                    162: {
                    163:   mpz_t t;
                    164:
                    165:   if ( !n1 || !n2 ) *nr = 0;
                    166:   else if ( !n1->z || !n2->z )
                    167:     error("mulz : invalid argument");
                    168:   else if ( UNIQ(n1) ) *nr = n2;
                    169:   else if ( UNIQ(n2) ) *nr = n1;
                    170:   else if ( MUNIQ(n1) ) chsgnz(n2,nr);
                    171:   else if ( MUNIQ(n2) ) chsgnz(n1,nr);
                    172:   else {
                    173:     mpz_init(t); mpz_mul(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nr);
                    174:   }
                    175: }
                    176:
                    177: /* nr += n1*n2 */
                    178:
                    179: void muladdtoz(Z n1,Z n2,Z *nr)
                    180: {
1.3       noro      181: #if 0
1.1       noro      182:   Z t;
                    183:
                    184:   if ( n1 && n2 ) {
                    185:         if ( !(*nr) ) {
                    186:           NEWZ(t); mpz_init(BDY(t)); *nr = t;
                    187:         }
                    188:         mpz_addmul(BDY(*nr),BDY(n1),BDY(n2));
1.2       noro      189:         if ( !mpz_sgn(BDY(*nr)) )
                    190:           *nr = 0;
1.3       noro      191:   }
1.2       noro      192: #else
                    193:   Z t,s;
                    194:
                    195:   mulz(n1,n2,&t); addz(*nr,t,&s); *nr = s;
                    196: #endif
1.1       noro      197: }
                    198:
                    199: /* nr += n1*u */
                    200:
                    201: void mul1addtoz(Z n1,long u,Z *nr)
                    202: {
1.3       noro      203: #if 0
1.1       noro      204:   Z t;
                    205:
                    206:   if ( n1 && u ) {
                    207:         if ( !(*nr) ) {
                    208:           NEWZ(t); mpz_init(BDY(t)); *nr = t;
                    209:         }
                    210:         if ( u >= 0 )
                    211:           mpz_addmul_ui(BDY(*nr),BDY(n1),(unsigned long)u);
                    212:         else
                    213:           mpz_submul_ui(BDY(*nr),BDY(n1),(unsigned long)(-u));
1.2       noro      214:         if ( !mpz_sgn(BDY(*nr)) )
                    215:           *nr = 0;
1.1       noro      216:     }
1.3       noro      217: #else
                    218:   Z t,s;
                    219:
                    220:   mul1z(n1,u,&t); addz(*nr,t,&s); *nr = s;
                    221: #endif
1.1       noro      222: }
                    223:
                    224: void mul1z(Z n1,long n2,Z *nr)
                    225: {
                    226:   mpz_t t;
                    227:
                    228:   if ( !n1 || !n2 ) *nr = 0;
                    229:   else {
                    230:     mpz_init(t); mpz_mul_si(t,BDY(n1),n2); MPZTOZ(t,*nr);
                    231:   }
                    232: }
                    233:
                    234: void divz(Z n1,Z n2,Z *nq)
                    235: {
                    236:   mpz_t t;
                    237:   mpq_t a, b, q;
                    238:
                    239:   if ( !n2 ) {
                    240:     error("division by 0");
                    241:     *nq = 0;
                    242:   } else if ( !n1 )
                    243:     *nq = 0;
                    244:   else if ( n1 == n2 ) {
                    245:     mpz_init(t); mpz_set_ui(t,1); MPZTOZ(t,*nq);
                    246:   } else {
                    247:     MPZTOMPQ(BDY(n1),a); MPZTOMPQ(BDY(n2),b);
                    248:     mpq_init(q); mpq_div(q,a,b); *nq = (Z)mpqtozq(q);
                    249:   }
                    250: }
                    251:
                    252: void remz(Z n1,Z n2,Z *nr)
                    253: {
                    254:   mpz_t r;
                    255:
                    256:   if ( !n2 ) {
                    257:     error("division by 0");
                    258:     *nr = 0;
                    259:   } else if ( !n1 || n1 == n2 )
                    260:     *nr = 0;
                    261:   else if ( !n1->z || !n2->z )
                    262:     error("remz : invalid argument");
                    263:   else {
                    264:     mpz_init(r);
                    265:     mpz_mod(r,BDY(n1),BDY(n2));
                    266:     if ( !mpz_sgn(r) ) *nr = 0;
                    267:     else MPZTOZ(r,*nr);
                    268:   }
                    269: }
                    270:
                    271: void divqrz(Z n1,Z n2,Z *nq,Z *nr)
                    272: {
                    273:   mpz_t t, a, b, q, r;
                    274:
                    275:   if ( !n2 ) {
                    276:     error("division by 0");
                    277:     *nq = 0; *nr = 0;
                    278:   } else if ( !n1 ) {
                    279:     *nq = 0; *nr = 0;
                    280:   } else if ( !n1->z || !n2->z )
                    281:     error("divqrz : invalid argument");
                    282:   else if ( n1 == n2 ) {
                    283:     mpz_init(t); mpz_set_ui(t,1); MPZTOZ(t,*nq); *nr = 0;
                    284:   } else {
                    285:     mpz_init(q); mpz_init(r);
                    286:     mpz_fdiv_qr(q,r,BDY(n1),BDY(n2));
                    287:     if ( !mpz_sgn(q) ) *nq = 0;
                    288:     else MPZTOZ(q,*nq);
                    289:     if ( !mpz_sgn(r) ) *nr = 0;
                    290:     else MPZTOZ(r,*nr);
                    291:   }
                    292: }
                    293:
                    294: void divsz(Z n1,Z n2,Z *nq)
                    295: {
                    296:   mpz_t t;
                    297:   mpq_t a, b, q;
                    298:
                    299:   if ( !n2 ) {
                    300:     error("division by 0");
                    301:     *nq = 0;
                    302:   } else if ( !n1 )
                    303:     *nq = 0;
                    304:   else if ( !n1->z || !n2->z )
                    305:     error("divsz : invalid argument");
                    306:   else if ( n1 == n2 ) {
                    307:     mpz_init(t); mpz_set_ui(t,1); MPZTOZ(t,*nq);
                    308:   } else {
                    309:     mpz_init(t); mpz_divexact(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nq);
                    310:   }
                    311: }
                    312:
                    313: void chsgnz(Z n,Z *nr)
                    314: {
                    315:   mpz_t t;
                    316:
                    317:   if ( !n )
                    318:     *nr = 0;
                    319:   else if ( !n->z )
                    320:     error("chsgnz : invalid argument");
                    321:   else {
                    322:     t[0] = BDY(n)[0]; mpz_neg(t,t); MPZTOZ(t,*nr);
                    323:   }
                    324: }
                    325:
                    326: void absz(Z n,Z *nr)
                    327: {
                    328:   if ( !n ) *nr = 0;
                    329:   else if ( !n->z )
                    330:     error("absz : invalid argument");
                    331:   else if ( sgnz(n) < 0 ) chsgnz(n,nr);
                    332:   else *nr = n;
                    333: }
                    334:
                    335: int evenz(Z n)
                    336: {
                    337:   return !n ? 1 : mpz_even_p(BDY(n));
                    338: }
                    339:
                    340: int smallz(Z n)
                    341: {
                    342:   if ( !n ) return 1;
                    343:   else if ( INT(n) && mpz_fits_sint_p(BDY(n)) ) return 1;
                    344:   else return 0;
                    345: }
                    346:
                    347: void pwrz(Z n1,Z n,Z *nr)
                    348: {
                    349:   mpq_t t,q;
                    350:   mpz_t z;
                    351:   Q p,r;
                    352:
                    353:   if ( !n || UNIQ(n1) ) *nr = ONE;
                    354:   else if ( !n1 ) *nr = 0;
                    355:   else if ( !n->z || !n1->z )
                    356:     error("pwrz : invalid argument");
                    357:   else if ( MUNIQ(n1) ) {
                    358:     if ( mpz_even_p(BDY((Z)n)) ) *nr = ONE;
                    359:     else *nr = n1;
                    360:   } else if ( !smallz(n) ) {
                    361:     error("exponent too big."); *nr = 0;
                    362:   } else if ( n1->z && mpz_sgn(BDY((Z)n))>0 ) {
1.5       noro      363:     mpz_init(z); mpz_pow_ui(z,BDY(n1),ZTOS(n)); MPZTOZ(z,*nr);
1.1       noro      364:   } else {
                    365:     MPZTOMPQ(BDY(n1),q); MPQTOQ(q,r);
                    366:     pwrq(r,(Q)n,&p); *nr = (Z)p;
                    367:   }
                    368: }
                    369:
                    370: int cmpz(Z q1,Z q2)
                    371: {
                    372:   int sgn;
                    373:
                    374:   if ( !q1 ) {
                    375:     if ( !q2 )
                    376:       return 0;
                    377:     else
                    378:       return -mpz_sgn(BDY(q2));
                    379:   } else if ( !q2 )
                    380:     return mpz_sgn(BDY(q1));
                    381:   else if ( !q1->z || !q2->z )
                    382:     error("mpqz : invalid argument");
                    383:   else if ( (sgn = mpz_sgn(BDY(q1))) != mpz_sgn(BDY(q2)) )
                    384:       return sgn;
                    385:   else {
                    386:     sgn = mpz_cmp(BDY(q1),BDY(q2));
                    387:     if ( sgn > 0 ) return 1;
                    388:     else if ( sgn < 0 ) return -1;
                    389:     else return 0;
                    390:   }
                    391: }
                    392:
                    393: void gcdz(Z n1,Z n2,Z *nq)
                    394: {
                    395:   mpz_t t;
                    396:
                    397:   if ( !n1 ) *nq = n2;
                    398:   else if ( !n2 ) *nq = n1;
                    399:   else if ( !n1->z || !n2->z )
                    400:     error("gcdz : invalid argument");
                    401:   else {
                    402:     mpz_init(t); mpz_gcd(t,BDY(n1),BDY(n2));
                    403:     MPZTOZ(t,*nq);
                    404:   }
                    405: }
                    406:
                    407: void invz(Z n1,Z n2,Z *nq)
                    408: {
                    409:   mpz_t t;
                    410:
                    411:   if ( !n1 || !n2 || !n1->z || !n2->z )
                    412:     error("invz : invalid argument");
                    413:   mpz_init(t); mpz_invert(t,BDY(n1),BDY(n2));
                    414:   MPZTOZ(t,*nq);
                    415: }
                    416:
                    417: void lcmz(Z n1,Z n2,Z *nq)
                    418: {
                    419:   Z g,t;
                    420:
                    421:   if ( !n1 || !n2 ) *nq = 0;
                    422:   else if ( !n1->z || !n2->z )
                    423:     error("lcmz : invalid argument");
                    424:   else {
                    425:     gcdz(n1,n2,&g); divsz(n1,g,&t);
                    426:     mulz(n2,t,nq);
                    427:   }
                    428: }
                    429:
                    430: void gcdvz(VECT v,Z *q)
                    431: {
                    432:   int n,i;
                    433:   Z *b;
                    434:   Z g,g1;
                    435:
                    436:   n = v->len;
                    437:   b = (Z *)v->body;
                    438:   g = b[0];
                    439:   for ( i = 1; i < n; i++ ) {
                    440:     gcdz(g,b[i],&g1); g = g1;
                    441:   }
                    442:   *q = g;
                    443: }
                    444:
                    445: void gcdvz_estimate(VECT v,Z *q)
                    446: {
                    447:   int n,m,i;
                    448:   Z s,t,u;
                    449:   Z *b;
                    450:
                    451:   n = v->len;
                    452:   b = (Z *)v->body;
                    453:   if ( n == 1 ) {
                    454:     if ( mpz_sgn(BDY(b[0]))<0 ) chsgnz(b[0],q);
                    455:     else *q = b[0];
                    456:   }
                    457:   m = n/2;
                    458:   for ( i = 0, s = 0; i < m; i++ ) {
                    459:     if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subz(s,b[i],&u);
                    460:     else addz(s,b[i],&u);
                    461:     s = u;
                    462:   }
1.4       noro      463:   for ( t = 0; i < n; i++ ) {
1.1       noro      464:     if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subz(t,b[i],&u);
                    465:     else addz(t,b[i],&u);
                    466:     t = u;
                    467:   }
                    468:   gcdz(s,t,q);
                    469: }
                    470:
1.4       noro      471: void gcdv_mpz_estimate(mpz_t g,mpz_t *b,int n)
                    472: {
                    473:   int m,m2,i,j;
                    474:   mpz_t s,t;
                    475:
                    476:   mpz_init(g);
                    477:   for ( i = 0, m = 0; i < n; i++ )
                    478:     if ( mpz_sgn(b[i]) ) m++;
                    479:   if ( !m ) {
                    480:     mpz_set_ui(g,0);
                    481:     return;
                    482:   }
                    483:   if ( m == 1 ) {
                    484:     for ( i = 0, m = 0; i < n; i++ )
                    485:       if ( mpz_sgn(b[i]) ) break;
                    486:     if ( mpz_sgn(b[i])<0 ) mpz_neg(g,b[i]);
                    487:     else mpz_set(g,b[i]);
                    488:     return ;
                    489:   }
                    490:   m2 = m/2;
                    491:   mpz_init_set_ui(s,0);
                    492:   for ( i = j = 0; j < m2; i++ ) {
                    493:     if ( mpz_sgn(b[i]) ) {
                    494:       if ( mpz_sgn(b[i])<0 )
                    495:         mpz_sub(s,s,b[i]);
                    496:       else
                    497:         mpz_add(s,s,b[i]);
                    498:       j++;
                    499:     }
                    500:   }
                    501:   mpz_init_set_ui(t,0);
                    502:   for ( ; i < n; i++ ) {
                    503:     if ( mpz_sgn(b[i]) ) {
                    504:       if ( mpz_sgn(b[i])<0 )
                    505:         mpz_sub(t,t,b[i]);
                    506:       else
                    507:         mpz_add(t,t,b[i]);
                    508:     }
                    509:   }
                    510:   mpz_gcd(g,s,t);
                    511: }
                    512:
                    513:
1.1       noro      514: void factorialz(unsigned int n,Z *nr)
                    515: {
                    516:   mpz_t a;
                    517:   mpz_init(a);
                    518:   mpz_fac_ui(a,n);
                    519:   MPZTOZ(a,*nr);
                    520: }
                    521:
                    522: void randomz(int blen,Z *nr)
                    523: {
                    524:   mpz_t z;
                    525:
                    526:   mpz_init(z);
                    527:   mpz_urandomb(z,GMP_RAND,blen);
                    528:   MPZTOZ(z,*nr);
                    529: }
                    530:
                    531: int tstbitz(Z n,int k)
                    532: {
                    533:    if ( !n || !n->z )
                    534:     error("tstbitz : invalid argument");
                    535:    return !n ? 0 : mpz_tstbit(BDY(n),k);
                    536: }
                    537:
                    538: void addq(Q n1,Q n2,Q *nr)
                    539: {
                    540:   mpq_t q1,q2,t;
                    541:
                    542:   if ( !n1 ) *nr = n2;
                    543:   else if ( !n2 ) *nr = n1;
                    544:   else if ( n1->z && n2->z )
                    545:     addz((Z)n1,(Z)n2,(Z *)nr);
                    546:   else {
                    547:     if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
                    548:     else q1[0] = BDY(n1)[0];
                    549:     if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
                    550:     else q2[0] = BDY(n2)[0];
                    551:     mpq_init(t); mpq_add(t,q1,q2); *nr = mpqtozq(t);
                    552:   }
                    553: }
                    554:
                    555: void subq(Q n1,Q n2,Q *nr)
                    556: {
                    557:   mpq_t q1,q2,t;
                    558:
                    559:   if ( !n1 ) {
                    560:     if ( !n2 ) *nr = 0;
                    561:     else if ( n1->z ) chsgnz((Z)n1,(Z *)nr);
                    562:     else {
                    563:         mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOQ(t,*nr);
                    564:       }
                    565:   } else if ( !n2 ) *nr = n1;
                    566:   else if ( n1 == n2 ) *nr = 0;
                    567:   else if ( n1->z && n2->z )
                    568:     subz((Z)n1,(Z)n2,(Z *)nr);
                    569:   else {
                    570:     if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
                    571:     else q1[0] = BDY(n1)[0];
                    572:     if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
                    573:     else q2[0] = BDY(n2)[0];
                    574:     mpq_init(t); mpq_sub(t,q1,q2); *nr = mpqtozq(t);
                    575:   }
                    576: }
                    577:
                    578: void mulq(Q n1,Q n2,Q *nr)
                    579: {
                    580:   mpq_t t,q1,q2;
                    581:
                    582:   if ( !n1 || !n2 ) *nr = 0;
                    583:   else if ( n1->z && n2->z )
                    584:     mulz((Z)n1,(Z)n2,(Z *)nr);
                    585:   else {
                    586:     if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
                    587:     else q1[0] = BDY(n1)[0];
                    588:     if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
                    589:     else q2[0] = BDY(n2)[0];
                    590:     mpq_init(t); mpq_mul(t,q1,q2); *nr = mpqtozq(t);
                    591:   }
                    592: }
                    593:
                    594: void divq(Q n1,Q n2,Q *nq)
                    595: {
                    596:   mpq_t t,q1,q2;
                    597:
                    598:   if ( !n2 ) {
                    599:     error("division by 0");
                    600:     *nq = 0;
                    601:     return;
                    602:   } else if ( !n1 ) *nq = 0;
                    603:   else if ( n1 == n2 ) *nq = (Q)ONE;
                    604:   else {
                    605:     if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
                    606:     else q1[0] = BDY(n1)[0];
                    607:     if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
                    608:     else q2[0] = BDY(n2)[0];
                    609:     mpq_init(t); mpq_div(t,q1,q2); *nq = mpqtozq(t);
                    610:   }
                    611: }
                    612:
                    613: void invq(Q n,Q *nr)
                    614: {
                    615:   Z nm,dn;
                    616:
                    617:   if ( INT(n) )
                    618:     divq((Q)ONE,n,nr);
                    619:   else {
                    620:     nmq(n,&nm);
                    621:     dnq(n,&dn);
                    622:     divq((Q)dn,(Q)nm,nr);
                    623:   }
                    624: }
                    625:
                    626: void chsgnq(Q n,Q *nr)
                    627: {
                    628:   mpq_t t;
                    629:
                    630:   if ( !n ) *nr = 0;
                    631:   else if (n->z ) chsgnz((Z)n,(Z *)nr);
                    632:   else {
                    633:     mpq_init(t); mpq_neg(t,BDY(n)); MPQTOQ(t,*nr);
                    634:   }
                    635: }
                    636:
                    637: void absq(Q n,Q *nr)
                    638: {
                    639:   if ( !n ) *nr = 0;
                    640:   else if ( n->z ) absz((Z)n,(Z *)nr);
                    641:   else if ( sgnq(n) < 0 ) chsgnq(n,nr);
                    642:   else *nr = n;
                    643: }
                    644:
                    645: void pwrq(Q n1,Q n,Q *nr)
                    646: {
                    647:   int e;
                    648:   mpz_t nm,dn;
                    649:   mpq_t t;
                    650:
                    651:   if ( !n || UNIQ((Z)n1) || UNIQ(n1) ) *nr = (Q)ONE;
                    652:   else if ( !n1 ) *nr = 0;
                    653:   else if ( !INT(n) ) {
                    654:     error("can't calculate fractional power."); *nr = 0;
                    655:   } else if ( !smallz((Z)n) ) {
                    656:     error("exponent too big."); *nr = 0;
                    657:   } else {
1.5       noro      658:     e = ZTOS(n);
1.1       noro      659:     if ( e < 0 ) {
                    660:       e = -e;
                    661:       if ( n1->z ) {
                    662:         nm[0] = ONEMPZ[0];
                    663:         dn[0] = BDY((Z)n1)[0];
                    664:       } else {
                    665:         nm[0] = mpq_denref(BDY(n1))[0];
                    666:         dn[0] = mpq_numref(BDY(n1))[0];
                    667:       }
                    668:     } else {
                    669:       if ( n1->z ) {
                    670:         nm[0] = BDY((Z)n1)[0];
                    671:         dn[0] = ONEMPZ[0];
                    672:       } else {
                    673:         nm[0] = mpq_numref(BDY(n1))[0];
                    674:         dn[0] = mpq_denref(BDY(n1))[0];
                    675:       }
                    676:     }
                    677:     mpq_init(t);
                    678:     mpz_pow_ui(mpq_numref(t),nm,e); mpz_pow_ui(mpq_denref(t),dn,e);
                    679:     *nr = mpqtozq(t);
                    680:   }
                    681: }
                    682:
                    683: int cmpq(Q n1,Q n2)
                    684: {
                    685:   mpq_t q1,q2;
                    686:   int sgn;
                    687:
                    688:   if ( !n1 ) {
                    689:     if ( !n2 ) return 0;
                    690:     else return (n2->z) ? -mpz_sgn(BDY((Z)n2)) : -mpq_sgn(BDY(n2));
                    691:   } if ( !n2 ) return (n1->z) ? mpz_sgn(BDY((Z)n1)) : mpq_sgn(BDY(n1));
                    692:   else if ( n1->z && n2->z )
                    693:     return cmpz((Z)n1,(Z)n2);
                    694:   else if ( (sgn = mpq_sgn(BDY(n1))) != mpq_sgn(BDY(n2)) ) return sgn;
                    695:   else {
                    696:     if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
                    697:     else q1[0] = BDY(n1)[0];
                    698:     if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
                    699:     else q2[0] = BDY(n2)[0];
                    700:     sgn = mpq_cmp(q1,q2);
                    701:     if ( sgn > 0 ) return 1;
                    702:     else if ( sgn < 0 ) return -1;
                    703:     else return 0;
                    704:   }
                    705: }
                    706:
                    707: /* t = [nC0 nC1 ... nCn] */
                    708:
                    709: void mkbc(int n,Z *t)
                    710: {
                    711:   int i;
                    712:   Z c,d,iq;
                    713:
                    714:   for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
1.5       noro      715:     STOZ(n-i+1,c); mulz(t[i-1],c,&d);
                    716:     STOZ(i,iq); divsz(d,iq,&t[i]);
1.1       noro      717:   }
                    718:   for ( ; i <= n; i++ )
                    719:     t[i] = t[n-i];
                    720: }
                    721:
                    722: /*
                    723:  *  Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
                    724:  *
                    725:  *  t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
                    726:  *  where W(k,l,i) = i! * kCi * lCi
                    727:  */
                    728:
                    729: /* mod m table */
                    730: /* XXX : should be optimized */
                    731:
                    732: void mkwcm(int k,int l,int m,int *t)
                    733: {
                    734:   int i,n;
                    735:   Z *s;
                    736:
                    737:   n = MIN(k,l);
                    738:   s = (Z *)ALLOCA((n+1)*sizeof(Q));
                    739:   mkwc(k,l,s);
                    740:   for ( i = 0; i <= n; i++ ) {
                    741:     t[i] = remqi((Q)s[i],m);
                    742:   }
                    743: }
                    744:
                    745: void mkwc(int k,int l,Z *t)
                    746: {
                    747:   mpz_t a,b,q,nm,z,u;
                    748:   int i,n;
                    749:
                    750:   n = MIN(k,l);
                    751:   mpz_init_set_ui(z,1);
                    752:   mpz_init(u); mpz_set(u,z); MPZTOZ(u,t[0]);
                    753:   mpz_init(a); mpz_init(b); mpz_init(nm);
                    754:   for ( i = 1; i <= n; i++ ) {
                    755:     mpz_set_ui(a,k-i+1); mpz_set_ui(b,l-i+1); mpz_mul(nm,a,b);
                    756:     mpz_mul(z,BDY(t[i-1]),nm); mpz_fdiv_q_ui(z,z,i);
                    757:     mpz_init(u); mpz_set(u,z); MPZTOZ(u,t[i]);
                    758:   }
                    759: }
                    760:
                    761: void lgp(P p,Z *g,Z *l);
                    762:
                    763: void ptozp(P p,int sgn,Q *c,P *pr)
                    764: {
                    765:   Z nm,dn;
                    766:
                    767:   if ( !p ) {
                    768:     *c = 0; *pr = 0;
                    769:   } else {
                    770:     lgp(p,&nm,&dn);
                    771:     divz(nm,dn,(Z *)c);
                    772:     divsp(CO,p,(P)*c,pr);
                    773:   }
                    774: }
                    775:
                    776: void lgp(P p,Z *g,Z *l)
                    777: {
                    778:   DCP dc;
                    779:   Z g1,g2,l1,l2,l3,l4;
                    780:
                    781:   if ( NUM(p) ) {
                    782:     if ( ((Q)p)->z ) {
                    783:       MPZTOZ(BDY((Z)p),*g);
                    784:       *l = ONE;
                    785:     } else {
                    786:       MPZTOZ(mpq_numref(BDY((Q)p)),*g);
                    787:       MPZTOZ(mpq_denref(BDY((Q)p)),*l);
                    788:     }
                    789:   } else {
                    790:     dc = DC(p); lgp(COEF(dc),g,l);
                    791:     for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
                    792:       lgp(COEF(dc),&g1,&l1); gcdz(*g,g1,&g2); *g = g2;
                    793:       gcdz(*l,l1,&l2); mulz(*l,l1,&l3); divz(l3,l2,l);
                    794:     }
                    795:   }
                    796: }
                    797:
                    798: void qltozl(Q *w,int n,Z *dvr)
                    799: {
                    800:   Z nm,dn;
                    801:   Z g,g1,l1,l2,l3;
                    802:   Q c;
                    803:   int i;
                    804:   struct oVECT v;
                    805:
                    806:   for ( i = 0; i < n; i++ )
                    807:     if ( w[i] && !w[i]->z )
                    808:       break;
                    809:   if ( i == n ) {
                    810:     v.id = O_VECT; v.len = n; v.body = (pointer *)w;
                    811:     gcdvz(&v,dvr); return;
                    812:   }
                    813:   for ( i = 0; !w[i]; i++ );
                    814:   c = w[i];
                    815:   if ( !c->z ) {
                    816:     MPZTOZ(mpq_numref(BDY(c)),nm); MPZTOZ(mpq_denref(BDY(c)),dn);
                    817:   } else {
                    818:     MPZTOZ(BDY((Z)c),nm); dn = ONE;
                    819:   }
                    820:   for ( i++; i < n; i++ ) {
                    821:     c = w[i];
                    822:     if ( !c ) continue;
                    823:     if ( !c->z ) {
                    824:       MPZTOZ(mpq_numref(BDY(c)),g1); MPZTOZ(mpq_denref(BDY(c)),l1);
                    825:     } else {
                    826:       MPZTOZ(BDY((Z)c),g1); l1 = ONE;
                    827:     }
                    828:     gcdz(nm,g1,&g); nm = g;
                    829:     gcdz(dn,l1,&l2); mulz(dn,l1,&l3); divz(l3,l2,&dn);
                    830:   }
                    831:   divz(nm,dn,dvr);
                    832: }
                    833:
                    834: int z_bits(Q q)
                    835: {
                    836:   if ( !q ) return 0;
                    837:   else if ( q->z ) return mpz_sizeinbase(BDY((Z)q),2);
                    838:   else
                    839:     return mpz_sizeinbase(mpq_numref(BDY(q)),2)
                    840:       + mpz_sizeinbase(mpq_denref(BDY(q)),2);
                    841: }
                    842:
                    843: int zp_mag(P p)
                    844: {
                    845:   int s;
                    846:   DCP dc;
                    847:
                    848:   if ( !p ) return 0;
                    849:   else if ( OID(p) == O_N ) return z_bits((Q)p);
                    850:   else {
                    851:     for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) s += zp_mag(COEF(dc));
                    852:     return s;
                    853:   }
                    854: }
                    855:
                    856: void makesubstz(VL v,NODE *s)
                    857: {
                    858:   NODE r,r0;
                    859:   Z q;
                    860:   unsigned int n;
                    861:
                    862:   for ( r0 = 0; v; v = NEXT(v) ) {
                    863:     NEXTNODE(r0,r); BDY(r) = (pointer)v->v;
                    864: #if defined(_PA_RISC1_1)
                    865:     n = mrand48()&BMASK; q = utoz(n);
                    866: #else
                    867:     n = random(); q = utoz(n);
                    868: #endif
                    869:     NEXTNODE(r0,r); BDY(r) = (pointer)q;
                    870:   }
                    871:   if ( r0 ) NEXT(r) = 0;
                    872:   *s = r0;
                    873: }
                    874:
                    875: unsigned int remqi(Q a,unsigned int mod)
                    876: {
                    877:   unsigned int c,nm,dn;
                    878:   mpz_t r;
                    879:
                    880:   if ( !a ) return 0;
                    881:   else if ( a->z ) {
                    882:     mpz_init(r);
                    883:     c = mpz_fdiv_r_ui(r,BDY((Z)a),mod);
                    884:   } else {
                    885:     mpz_init(r);
                    886:     nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod);
                    887:     dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod);
                    888:     dn = invm(dn,mod);
                    889:     DMAR(nm,dn,0,mod,c);
                    890:   }
                    891:   return c;
                    892: }
                    893:
                    894: int generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
                    895: {
                    896:   int **wmat;
                    897:   Z **bmat,**tmat,*bmi,*tmi;
                    898:   Z q,m1,m2,m3,s,u;
                    899:   int *wmi,*colstat,*wcolstat,*rind,*cind;
                    900:   int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
                    901:   MAT r,crmat;
                    902:   int ret;
                    903:
1.6       noro      904: #if SIZEOF_LONG == 8
                    905:   return generic_gauss_elim64(mat,nm,dn,rindp,cindp);
                    906: #endif
1.1       noro      907:   bmat = (Z **)mat->body;
                    908:   row = mat->row; col = mat->col;
                    909:   wmat = (int **)almat(row,col);
                    910:   colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                    911:   wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                    912:   for ( ind = 0; ; ind++ ) {
                    913:     if ( DP_Print ) {
                    914:       fprintf(asir_out,"."); fflush(asir_out);
                    915:     }
                    916:     md = get_lprime(ind);
                    917:     for ( i = 0; i < row; i++ )
                    918:       for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
                    919:         wmi[j] = remqi((Q)bmi[j],md);
                    920:     rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
                    921:     if ( !ind ) {
                    922: RESET:
                    923:       m1 = utoz(md);
                    924:       rank0 = rank;
                    925:       bcopy(wcolstat,colstat,col*sizeof(int));
                    926:       MKMAT(crmat,rank,col-rank);
                    927:       MKMAT(r,rank,col-rank); *nm = r;
                    928:       tmat = (Z **)crmat->body;
                    929:       for ( i = 0; i < rank; i++ )
                    930:         for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                    931:           if ( !colstat[j] ) tmi[k++] = utoz(wmi[j]);
                    932:     } else {
                    933:       if ( rank < rank0 ) {
                    934:         if ( DP_Print ) {
                    935:           fprintf(asir_out,"lower rank matrix; continuing...\n");
                    936:           fflush(asir_out);
                    937:         }
                    938:         continue;
                    939:       } else if ( rank > rank0 ) {
                    940:         if ( DP_Print ) {
                    941:           fprintf(asir_out,"higher rank matrix; resetting...\n");
                    942:           fflush(asir_out);
                    943:         }
                    944:         goto RESET;
                    945:       } else {
                    946:         for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
                    947:         if ( j < col ) {
                    948:           if ( DP_Print ) {
                    949:             fprintf(asir_out,"inconsitent colstat; resetting...\n");
                    950:             fflush(asir_out);
                    951:           }
                    952:           goto RESET;
                    953:         }
                    954:       }
                    955:
                    956:       inv = invm(remqi((Q)m1,md),md);
                    957:       m2 = utoz(md); mulz(m1,m2,&m3);
                    958:       for ( i = 0; i < rank; i++ )
                    959:         for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                    960:           if ( !colstat[j] ) {
                    961:             if ( tmi[k] ) {
                    962:             /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                    963:               t = remqi((Q)tmi[k],md);
                    964:               if ( wmi[j] >= t ) t = wmi[j]-t;
                    965:               else t = md-(t-wmi[j]);
                    966:               DMAR(t,inv,0,md,t1)
                    967:               u = utoz(t1); mulz(m1,u,&s);
                    968:               addz(tmi[k],s,&u); tmi[k] = u;
                    969:             } else if ( wmi[j] ) {
                    970:             /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                    971:               DMAR(wmi[j],inv,0,md,t)
                    972:               u = utoz(t); mulz(m1,u,&s); tmi[k] = s;
                    973:             }
                    974:             k++;
                    975:           }
                    976:       m1 = m3;
                    977:       if ( ind % F4_INTRAT_PERIOD )
                    978:         ret = 0;
                    979:       else
                    980:         ret = intmtoratm(crmat,m1,*nm,dn);
                    981:       if ( ret ) {
                    982:         *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                    983:         *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
                    984:         for ( j = k = l = 0; j < col; j++ )
                    985:           if ( colstat[j] ) rind[k++] = j;
                    986:           else cind[l++] = j;
                    987:         if ( gensolve_check(mat,*nm,*dn,rind,cind) )
                    988:           return rank;
                    989:       }
                    990:     }
                    991:   }
                    992: }
                    993:
                    994: int generic_gauss_elim2(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
                    995: {
                    996:
                    997:   MAT full;
                    998:   Z **bmat,**b;
                    999:   Z *bmi;
                   1000:   Z dn0;
                   1001:   int row,col,md,i,j,rank,ret;
                   1002:   int **wmat;
                   1003:   int *wmi;
                   1004:   int *colstat,*rowstat;
                   1005:
                   1006:   bmat = (Z **)mat->body;
                   1007:   row = mat->row; col = mat->col;
                   1008:   wmat = (int **)almat(row,col);
                   1009:   colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                   1010:   rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
                   1011:   /* XXX */
                   1012:   md = get_lprime(0);
                   1013:   for ( i = 0; i < row; i++ )
                   1014:     for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
                   1015:       wmi[j] = remqi((Q)bmi[j],md);
                   1016:   rank = generic_gauss_elim_mod2(wmat,row,col,md,colstat,rowstat);
                   1017:   b = (Z **)MALLOC(rank*sizeof(Z));
                   1018:   for ( i = 0; i < rank; i++ ) b[i] = bmat[rowstat[i]];
                   1019:   NEWMAT(full); full->row = rank; full->col = col; full->body = (pointer **)b;
                   1020:   ret = generic_gauss_elim_full(full,nm,dn,rindp,cindp);
                   1021:   if ( !ret ) {
                   1022:     rank = generic_gauss_elim(mat,nm,&dn0,rindp,cindp);
                   1023:     for ( i = 0; i < rank; i++ ) dn[i] = dn0;
                   1024:   }
                   1025:   return rank;
                   1026: }
                   1027:
                   1028: int generic_gauss_elim_full(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
                   1029: {
                   1030:   int **wmat;
                   1031:   int *stat;
                   1032:   Z **bmat,**tmat,*bmi,*tmi,*ri;
                   1033:   Z q,m1,m2,m3,s,u;
                   1034:   int *wmi,*colstat,*wcolstat,*rind,*cind;
                   1035:   int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv,h;
                   1036:   MAT r,crmat;
                   1037:   int ret,initialized,done;
                   1038:
                   1039:   initialized = 0;
                   1040:   bmat = (Z **)mat->body;
                   1041:   row = mat->row; col = mat->col;
                   1042:   wmat = (int **)almat(row,col);
                   1043:   stat = (int *)MALLOC_ATOMIC(row*sizeof(int));
                   1044:   for ( i = 0; i < row; i++ ) stat[i] = 0;
                   1045:   colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                   1046:   wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                   1047:   for ( ind = 0; ; ind++ ) {
                   1048:     if ( DP_Print ) {
                   1049:       fprintf(asir_out,"."); fflush(asir_out);
                   1050:     }
                   1051:     md = get_lprime(ind);
                   1052:     for ( i = 0; i < row; i++ )
                   1053:       for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
                   1054:         wmi[j] = remqi((Q)bmi[j],md);
                   1055:     rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
                   1056:     if ( rank < row ) continue;
                   1057:     if ( !initialized ) {
                   1058:       m1 = utoz(md);
                   1059:       bcopy(wcolstat,colstat,col*sizeof(int));
                   1060:       MKMAT(crmat,row,col-row);
                   1061:       MKMAT(r,row,col-row); *nm = r;
                   1062:       tmat = (Z **)crmat->body;
                   1063:       for ( i = 0; i < row; i++ )
                   1064:         for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                   1065:           if ( !colstat[j] ) tmi[k++] = utoz(wmi[j]);
                   1066:       initialized = 1;
                   1067:     } else {
                   1068:       for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
                   1069:       if ( j < col ) continue;
                   1070:
                   1071:       inv = invm(remqi((Q)m1,md),md);
                   1072:       m2 = utoz(md); mulz(m1,m2,&m3);
                   1073:       for ( i = 0; i < row; i++ )
                   1074:         switch ( stat[i] ) {
                   1075:           case 1:
                   1076:             /* consistency check */
                   1077:             ri = (Z *)BDY(r)[i]; wmi = wmat[i];
                   1078:             for ( j = 0; j < col; j++ ) if ( colstat[j] ) break;
                   1079:             h = md-remqi((Q)dn[i],md);
                   1080:             for ( j++, k = 0; j < col; j++ )
                   1081:               if ( !colstat[j] ) {
                   1082:                 t = remqi((Q)ri[k],md);
                   1083:                 DMAR(wmi[i],h,t,md,t1);
                   1084:                 if ( t1 ) break;
                   1085:               }
                   1086:             if ( j == col ) { stat[i]++; break; }
                   1087:             else {
                   1088:               /* fall to the case 0 */
                   1089:               stat[i] = 0;
                   1090:             }
                   1091:           case 0:
                   1092:             tmi = tmat[i]; wmi = wmat[i];
                   1093:             for ( j = k = 0; j < col; j++ )
                   1094:               if ( !colstat[j] ) {
                   1095:                 if ( tmi[k] ) {
                   1096:                 /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                   1097:                   t = remqi((Q)tmi[k],md);
                   1098:                   if ( wmi[j] >= t ) t = wmi[j]-t;
                   1099:                   else t = md-(t-wmi[j]);
                   1100:                   DMAR(t,inv,0,md,t1)
                   1101:                   u = utoz(t1); mulz(m1,u,&s);
                   1102:                   addz(tmi[k],s,&u); tmi[k] = u;
                   1103:                 } else if ( wmi[j] ) {
                   1104:                 /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                   1105:                   DMAR(wmi[j],inv,0,md,t)
                   1106:                   u = utoz(t); mulz(m1,u,&s); tmi[k] = s;
                   1107:                 }
                   1108:                 k++;
                   1109:               }
                   1110:             break;
                   1111:           case 2: default:
                   1112:             break;
                   1113:         }
                   1114:       m1 = m3;
                   1115:       if ( ind % 4 )
                   1116:         ret = 0;
                   1117:       else
                   1118:         ret = intmtoratm2(crmat,m1,*nm,dn,stat);
                   1119:       if ( ret ) {
                   1120:         *rindp = rind = (int *)MALLOC_ATOMIC(row*sizeof(int));
                   1121:         *cindp = cind = (int *)MALLOC_ATOMIC((col-row)*sizeof(int));
                   1122:         for ( j = k = l = 0; j < col; j++ )
                   1123:           if ( colstat[j] ) rind[k++] = j;
                   1124:           else cind[l++] = j;
                   1125:         return gensolve_check2(mat,*nm,dn,rind,cind);
                   1126:       }
                   1127:     }
                   1128:   }
                   1129: }
                   1130:
                   1131: int generic_gauss_elim_direct(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp){
                   1132:   Z **bmat,*s;
                   1133:   Z u,v,w,x,d,t,y;
                   1134:   int row,col,i,j,k,l,m,rank;
                   1135:   int *colstat,*colpos,*cind;
                   1136:   MAT r,in;
                   1137:
                   1138:   row = mat->row; col = mat->col;
                   1139:   MKMAT(in,row,col);
                   1140:   for ( i = 0; i < row; i++ )
                   1141:     for ( j = 0; j < col; j++ ) in->body[i][j] = mat->body[i][j];
                   1142:   bmat = (Z **)in->body;
                   1143:   colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                   1144:   *rindp = colpos = (int *)MALLOC_ATOMIC(row*sizeof(int));
                   1145:   for ( j = 0, rank = 0, d = ONE; j < col; j++ ) {
                   1146:     for ( i = rank; i < row && !bmat[i][j]; i++  );
                   1147:     if ( i == row ) { colstat[j] = 0; continue; }
                   1148:     else { colstat[j] = 1; colpos[rank] = j; }
                   1149:     if ( i != rank )
                   1150:       for ( k = j; k < col; k++ ) {
                   1151:         t = bmat[i][k]; bmat[i][k] = bmat[rank][k]; bmat[rank][k] = t;
                   1152:       }
                   1153:     for ( i = rank+1, v = bmat[rank][j]; i < row; i++ )
                   1154:       for ( k = j, u = bmat[i][j]; k < col; k++ ) {
                   1155:         mulz(bmat[i][k],v,&w); mulz(bmat[rank][k],u,&x);
                   1156:         subz(w,x,&y); divsz(y,d,&bmat[i][k]);
                   1157:       }
                   1158:     d = v; rank++;
                   1159:   }
                   1160:   *dn = d;
                   1161:   s = (Z *)MALLOC(col*sizeof(Z));
                   1162:   for ( i = rank-1; i >= 0; i-- ) {
                   1163:     for ( k = colpos[i]; k < col; k++ ) mulz(bmat[i][k],d,&s[k]);
                   1164:     for ( m = rank-1; m > i; m-- ) {
                   1165:       for ( k = colpos[m], u = bmat[i][k]; k < col; k++ ) {
                   1166:         mulz(bmat[m][k],u,&w); subz(s[k],w,&x); s[k] = x;
                   1167:       }
                   1168:     }
                   1169:     for ( k = colpos[i], u = bmat[i][k]; k < col; k++ )
                   1170:       divz(s[k],u,&bmat[i][k]);
                   1171:   }
                   1172:   *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
                   1173:   MKMAT(r,rank,col-rank); *nm = r;
                   1174:   for ( j = 0, k = 0; j < col; j++ )
                   1175:     if ( !colstat[j] ) {
                   1176:       cind[k] = j;
                   1177:       for ( i = 0; i < rank; i++ ) r->body[i][k] = bmat[i][j];
                   1178:       k++;
                   1179:     }
                   1180:   return rank;
                   1181: }
                   1182:
                   1183: int intmtoratm(MAT mat,Z md,MAT nm,Z *dn)
                   1184: {
                   1185:   Z t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy;
                   1186:   int i,j,k,l,row,col,sgn,ret;
                   1187:   Z **rmat,**tmat,*tmi,*nmk;
                   1188:
                   1189:   if ( UNIQ(md) )
                   1190:     return 0;
                   1191:   row = mat->row; col = mat->col;
                   1192:   bshiftz(md,1,&t);
                   1193:   isqrtz(t,&s);
                   1194:   bshiftz(s,64,&b);
                   1195:   if ( !b ) b = ONE;
                   1196:   dn0 = ONE;
                   1197:   tmat = (Z **)mat->body;
                   1198:   rmat = (Z **)nm->body;
                   1199:   for ( i = 0; i < row; i++ )
                   1200:     for ( j = 0, tmi = tmat[i]; j < col; j++ )
                   1201:       if ( tmi[j] ) {
                   1202:         mulz(tmi[j],dn0,&s);
                   1203:         divqrz(s,md,&dmy,&u);
                   1204:         ret = inttorat(u,md,b,&nm1,&dn1);
                   1205:         if ( !ret ) return 0;
                   1206:         else {
                   1207:           if ( !UNIQ(dn1) ) {
                   1208:             for ( k = 0; k < i; k++ )
                   1209:               for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
                   1210:                 mulz(nmk[l],dn1,&q); nmk[l] = q;
                   1211:               }
                   1212:             for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
                   1213:               mulz(nmk[l],dn1,&q); nmk[l] = q;
                   1214:             }
                   1215:           }
                   1216:           rmat[i][j] = nm1;
                   1217:           mulz(dn0,dn1,&q); dn0 = q;
                   1218:         }
                   1219:       }
                   1220:   *dn = dn0;
                   1221:   return 1;
                   1222: }
                   1223:
                   1224: int intmtoratm2(MAT mat,Z md,MAT nm,Z *dn,int *stat)
                   1225: {
                   1226:   int row,col,i,j,ret;
                   1227:   Z dn0,dn1,t,s,b;
                   1228:   Z *w,*tmi;
                   1229:   Z **tmat;
                   1230:
                   1231:   bshiftz(md,1,&t);
                   1232:   isqrtz(t,&s);
                   1233:   bshiftz(s,64,&b);
                   1234:   tmat = (Z **)mat->body;
                   1235:   if ( UNIQ(md) ) return 0;
                   1236:   row = mat->row; col = mat->col;
                   1237:   dn0 = ONE;
                   1238:   for ( i = 0; i < row; i++ )
                   1239:     if ( cmpz(dn[i],dn0) > 0 ) dn0 = dn[i];
                   1240:   w = (Z *)MALLOC(col*sizeof(Z));
                   1241:   for ( i = 0; i < row; i++ )
                   1242:     if ( stat[i] == 0 ) {
                   1243:       for ( j = 0, tmi = tmat[i]; j < col; j++ )
                   1244:           mulz(tmi[j],dn0,&w[j]);
                   1245:       ret = intvtoratv(w,col,md,b,(Z *)BDY(nm)[i],&dn[i]);
                   1246:       if ( ret ) {
                   1247:         stat[i] = 1;
                   1248:         mulz(dn0,dn[i],&t); dn[i] = t; dn0 = t;
                   1249:       }
                   1250:     }
                   1251:   for ( i = 0; i < row; i++ ) if ( !stat[i] ) break;
                   1252:   if ( i == row ) return 1;
                   1253:   else return 0;
                   1254: }
                   1255:
                   1256: int intvtoratv(Z *v,int n,Z md,Z b,Z *nm,Z *dn)
                   1257: {
                   1258:   Z dn0,dn1,q,s,u,nm1,unm,udn,dmy;
                   1259:   Z *nmk;
                   1260:   int j,l,col,ret,sgn;
                   1261:
                   1262:   for ( j = 0; j < n; j++ ) nm[j] = 0;
                   1263:   dn0 = ONE;
                   1264:   for ( j = 0; j < n; j++ ) {
                   1265:     if ( !v[j] ) continue;
                   1266:     mulz(v[j],dn0,&s);
                   1267:     divqrz(s,md,&dmy,&u);
                   1268:     ret = inttorat(u,md,b,&nm1,&dn1);
                   1269:     if ( !ret ) return 0;
                   1270:     if ( !UNIQ(dn1) )
                   1271:       for ( l = 0; l < j; l++ ) {
                   1272:         mulz(nm[l],dn1,&q); nm[l] = q;
                   1273:       }
                   1274:     nm[j] = nm1;
                   1275:     mulz(dn0,dn1,&q); dn0 = q;
                   1276:   }
                   1277:   *dn = dn0;
                   1278:   return 1;
                   1279: }
                   1280:
                   1281: /* assuming 0 < c < m */
                   1282:
                   1283: int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
                   1284: {
                   1285:   Z qq,t,u1,v1,r1;
                   1286:   Z q,u2,v2,r2;
                   1287:
                   1288:   u1 = 0; v1 = ONE; u2 = m; v2 = c;
                   1289:   while ( cmpz(v2,b) >= 0 ) {
                   1290:     divqrz(u2,v2,&q,&r2); u2 = v2; v2 = r2;
                   1291:     mulz(q,v1,&t); subz(u1,t,&r1); u1 = v1; v1 = r1;
                   1292:   }
                   1293:   if ( cmpz(v1,b) >= 0 ) return 0;
                   1294:   else {
                   1295:     if ( mpz_sgn(BDY(v1))<0  ) {
                   1296:       chsgnz(v1,dnp); chsgnz(v2,nmp);
                   1297:     } else {
                   1298:       *dnp = v1; *nmp = v2;
                   1299:     }
                   1300:     return 1;
                   1301:   }
                   1302: }
                   1303:
                   1304: extern int f4_nocheck;
                   1305:
                   1306: int gensolve_check(MAT mat,MAT nm,Z dn,int *rind,int *cind)
                   1307: {
                   1308:   int row,col,rank,clen,i,j,k,l;
                   1309:   Z s,t;
                   1310:   Z *w;
                   1311:   Z *mati,*nmk;
                   1312:
                   1313:   if ( f4_nocheck ) return 1;
                   1314:   row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
                   1315:   w = (Z *)MALLOC(clen*sizeof(Z));
                   1316:   for ( i = 0; i < row; i++ ) {
                   1317:     mati = (Z *)mat->body[i];
                   1318:     bzero(w,clen*sizeof(Z));
                   1319:     for ( k = 0; k < rank; k++ )
                   1320:       for ( l = 0, nmk = (Z *)nm->body[k]; l < clen; l++ ) {
                   1321:         mulz(mati[rind[k]],nmk[l],&t); addz(w[l],t,&s); w[l] = s;
                   1322:       }
                   1323:     for ( j = 0; j < clen; j++ ) {
                   1324:       mulz(dn,mati[cind[j]],&t);
                   1325:       if ( cmpz(w[j],t) ) break;
                   1326:     }
                   1327:     if ( j != clen ) break;
                   1328:   }
                   1329:   if ( i != row ) return 0;
                   1330:   else return 1;
                   1331: }
                   1332:
                   1333: int gensolve_check2(MAT mat,MAT nm,Z *dn,int *rind,int *cind)
                   1334: {
                   1335:   int row,col,rank,clen,i,j,k,l;
                   1336:   Z s,t,u,d;
                   1337:   Z *w,*m;
                   1338:   Z *mati,*nmk;
                   1339:
                   1340:   if ( f4_nocheck ) return 1;
                   1341:   row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
                   1342:   w = (Z *)MALLOC(clen*sizeof(Z));
                   1343:   m = (Z *)MALLOC(clen*sizeof(Z));
                   1344:   for ( d = dn[0], i = 1; i < rank; i++ ) {
                   1345:     lcmz(d,dn[i],&t); d = t;
                   1346:   }
                   1347:   for ( i = 0; i < rank; i++ ) divsz(d,dn[i],&m[i]);
                   1348:   for ( i = 0; i < row; i++ ) {
                   1349:     mati = (Z *)mat->body[i];
                   1350:     bzero(w,clen*sizeof(Z));
                   1351:     for ( k = 0; k < rank; k++ ) {
                   1352:       mulz(mati[rind[k]],m[k],&u);
                   1353:       for ( l = 0, nmk = (Z *)nm->body[k]; l < clen; l++ ) {
                   1354:         mulz(u,nmk[l],&t); addz(w[l],t,&s); w[l] = s;
                   1355:       }
                   1356:     }
                   1357:     for ( j = 0; j < clen; j++ ) {
                   1358:       mulz(d,mati[cind[j]],&t);
                   1359:       if ( cmpz(w[j],t) ) break;
                   1360:     }
                   1361:     if ( j != clen ) break;
                   1362:   }
                   1363:   if ( i != row ) return 0;
                   1364:   else return 1;
                   1365: }
                   1366:
                   1367: void isqrtz(Z a,Z *r)
                   1368: {
                   1369:   int k;
                   1370:   Z x,t,x2,xh,quo,rem;
                   1371:   Z two;
                   1372:
                   1373:   if ( !a ) *r = 0;
                   1374:   else if ( UNIQ(a) ) *r = ONE;
                   1375:   else {
                   1376:     k = z_bits((Q)a); /* a <= 2^k-1 */
                   1377:     bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */
1.5       noro     1378:     STOZ(2,two);
1.1       noro     1379:     while ( 1 ) {
                   1380:       pwrz(x,two,&t);
                   1381:       if ( cmpz(t,a) <= 0 ) {
                   1382:         *r = x; return;
                   1383:       } else {
                   1384:         if ( mpz_tstbit(BDY(x),0) ) addz(x,a,&t);
                   1385:         else t = a;
                   1386:         bshiftz(x,-1,&x2); divqrz(t,x2,&quo,&rem);
                   1387:         bshiftz(x,1,&xh); addz(quo,xh,&x);
                   1388:       }
                   1389:     }
                   1390:   }
                   1391: }
                   1392:
                   1393: void bshiftz(Z a,int n,Z *r)
                   1394: {
                   1395:   mpz_t t;
                   1396:
                   1397:   if ( !a ) *r = 0;
                   1398:   else if ( n == 0 ) *r = a;
                   1399:   else if ( n < 0 ) {
                   1400:     mpz_init(t); mpz_mul_2exp(t,BDY(a),-n); MPZTOZ(t,*r);
                   1401:   } else {
                   1402:     mpz_init(t); mpz_fdiv_q_2exp(t,BDY(a),n);
                   1403:     if ( !mpz_sgn(t) ) *r = 0;
                   1404:     else MPZTOZ(t,*r);
                   1405:   }
                   1406: }
                   1407:
                   1408: void addlf(Z a,Z b,Z *c)
                   1409: {
                   1410:   addz(a,b,c);
                   1411:   if ( !lf_lazy ) {
                   1412:     if ( cmpz(*c,current_mod_lf) >= 0 ) {
                   1413:       subz(*c,current_mod_lf,c);
                   1414:     }
                   1415:   }
                   1416: }
                   1417:
                   1418: void sublf(Z a,Z b,Z *c)
                   1419: {
                   1420:   subz(a,b,c);
                   1421:   if ( !lf_lazy ) {
                   1422:     remz(*c,current_mod_lf,c);
                   1423:   }
                   1424: }
                   1425:
                   1426: void mullf(Z a,Z b,Z *c)
                   1427: {
                   1428:   mulz(a,b,c);
                   1429:   if ( !lf_lazy ) {
                   1430:     remz(*c,current_mod_lf,c);
                   1431:   }
                   1432: }
                   1433:
                   1434: void divlf(Z a,Z b,Z *c)
                   1435: {
                   1436:   Z inv;
                   1437:
                   1438:   invz(b,current_mod_lf,&inv);
                   1439:   mulz(a,inv,c);
                   1440:   if ( !lf_lazy ) {
                   1441:     remz(*c,current_mod_lf,c);
                   1442:   }
                   1443: }
                   1444:
                   1445: void chsgnlf(Z a,Z *c)
                   1446: {
                   1447:   chsgnz(a,c);
                   1448:   if ( !lf_lazy ) {
                   1449:     remz(*c,current_mod_lf,c);
                   1450:   }
                   1451: }
                   1452:
                   1453: void lmtolf(LM a,Z *b)
                   1454: {
                   1455:   if ( !a ) *b = 0;
                   1456:   else {
                   1457:     MPZTOZ(BDY(a),*b);
                   1458:   }
                   1459: }
                   1460:
                   1461: void setmod_lf(Z p)
                   1462: {
                   1463:     current_mod_lf = p;
                   1464:     current_mod_lf_size = mpz_size(BDY(current_mod_lf))+1;
                   1465: }
                   1466:
                   1467: void simplf_force(Z a,Z *b)
                   1468: {
                   1469:     remz(a,current_mod_lf,b);
                   1470: }
                   1471:
                   1472: int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp)
                   1473: {
                   1474:   MAT bmat,xmat;
                   1475:   Z **a0,**a,**b,**x,**nm;
                   1476:   Z *ai,*bi,*xi;
                   1477:   int row,col;
                   1478:   int **w;
                   1479:   int *wi;
                   1480:   int **wc;
                   1481:   Z mdq,q,s,u;
                   1482:   Z tn;
                   1483:   int ind,md,i,j,k,l,li,ri,rank;
                   1484:   unsigned int t;
                   1485:   int *cinfo,*rinfo;
                   1486:   int *rind,*cind;
                   1487:   int count;
                   1488:   int ret;
1.3       noro     1489:   struct oEGT eg_mul1,eg_mul2,tmp0,tmp1,tmp2;
1.1       noro     1490:   int period;
                   1491:   int *wx,*ptr;
                   1492:   int wxsize,nsize;
                   1493:   Z wn;
                   1494:   Z wq;
                   1495:
1.3       noro     1496: init_eg(&eg_mul1); init_eg(&eg_mul2);
1.1       noro     1497:   a0 = (Z **)mat->body;
                   1498:   row = mat->row; col = mat->col;
                   1499:   w = (int **)almat(row,col);
                   1500:   for ( ind = 0; ; ind++ ) {
                   1501:     md = get_lprime(ind);
1.5       noro     1502:     STOZ(md,mdq);
1.1       noro     1503:     for ( i = 0; i < row; i++ )
                   1504:       for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
                   1505:         wi[j] = remqi((Q)ai[j],md);
                   1506:
                   1507:     if ( DP_Print > 3 ) {
                   1508:       fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
                   1509:     }
                   1510:     rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
                   1511:     if ( DP_Print > 3 ) {
                   1512:       fprintf(asir_out,"done.\n"); fflush(asir_out);
                   1513:     }
                   1514:     a = (Z **)almat_pointer(rank,rank); /* lhs mat */
                   1515:     MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */
                   1516:     for ( j = li = ri = 0; j < col; j++ )
                   1517:       if ( cinfo[j] ) {
                   1518:         /* the column is in lhs */
                   1519:         for ( i = 0; i < rank; i++ ) {
                   1520:           w[i][li] = w[i][j];
                   1521:           a[i][li] = a0[rinfo[i]][j];
                   1522:         }
                   1523:         li++;
                   1524:       } else {
                   1525:         /* the column is in rhs */
                   1526:         for ( i = 0; i < rank; i++ )
                   1527:           b[i][ri] = a0[rinfo[i]][j];
                   1528:         ri++;
                   1529:       }
                   1530:
                   1531:       /* solve Ax=B; A: rank x rank, B: rank x ri */
                   1532:       /* algorithm
                   1533:          c <- B
                   1534:          x <- 0
                   1535:          q <- 1
                   1536:          do
                   1537:            t <- A^(-1)c mod p
                   1538:            x <- x+qt
                   1539:            c <- (c-At)/p
                   1540:            q <- qp
                   1541:          end do
                   1542:          then Ax-B=0 mod q and b=(B-Ax)/q hold after "do".
                   1543:       */
                   1544:       MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;
                   1545:       MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;
                   1546:       wc = (int **)almat(rank,ri);
                   1547:       *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                   1548:       *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
                   1549:
                   1550:       period = F4_INTRAT_PERIOD;
                   1551:       for ( q = ONE, count = 0; ; ) {
1.3       noro     1552:         /* check Ax=B mod q */
1.1       noro     1553:         if ( DP_Print > 3 )
                   1554:           fprintf(stderr,"o");
                   1555:         /* wc = b mod md */
                   1556:         for ( i = 0; i < rank; i++ )
1.3       noro     1557:           for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
1.1       noro     1558:             wi[j] = remqi((Q)bi[j],md);
1.3       noro     1559:         /* wc = A^(-1)wc; wc is not normalized */
                   1560:         solve_by_lu_mod(w,rank,md,wc,ri,0);
1.1       noro     1561:         /* x += q*wc */
1.3       noro     1562: get_eg(&tmp0);
1.1       noro     1563:         for ( i = 0; i < rank; i++ )
                   1564:           for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);
1.3       noro     1565:         /* b =(b-A*wc)/md */
                   1566: get_eg(&tmp1); add_eg(&eg_mul1,&tmp0,&tmp1);
1.1       noro     1567:         for ( i = 0; i < rank; i++ )
                   1568:           for ( j = 0; j < ri; j++ ) {
1.3       noro     1569:             mpz_t uz;
                   1570:
                   1571:             if ( b[i][j] )
                   1572:               mpz_init_set(uz,BDY(b[i][j]));
                   1573:             else
                   1574:               mpz_init_set_ui(uz,0);
                   1575:             for ( k = 0; k < rank; k++ ) {
                   1576:               if ( a[i][k] && wc[k][j] ) {
                   1577:                 if ( wc[k][j] < 0 )
                   1578:                   mpz_addmul_ui(uz,BDY(a[i][k]),-wc[k][j]);
                   1579:                 else
                   1580:                   mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]);
                   1581:               }
                   1582:             }
                   1583:             MPZTOZ(uz,u);
1.1       noro     1584:             divsz(u,mdq,&b[i][j]);
                   1585:           }
1.3       noro     1586: get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2);
1.1       noro     1587:         count++;
                   1588:         /* q = q*md */
                   1589:         mulz(q,mdq,&u); q = u;
                   1590:         if ( count == period ) {
                   1591:           ret = intmtoratm(xmat,q,*nmmat,dn);
                   1592:           if ( ret ) {
1.3       noro     1593:             print_eg("MUL1",&eg_mul1);
                   1594:             print_eg("MUL2",&eg_mul2);
1.1       noro     1595:             for ( j = k = l = 0; j < col; j++ )
                   1596:               if ( cinfo[j] )
                   1597:                 rind[k++] = j;
                   1598:               else
                   1599:                 cind[l++] = j;
                   1600:             ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
                   1601:             if ( ret ) {
                   1602:               *rindp = rind;
                   1603:               *cindp = cind;
                   1604:               for ( j = k = 0; j < col; j++ )
                   1605:                 if ( !cinfo[j] )
                   1606:                   cind[k++] = j;
                   1607:               return rank;
                   1608:             }
                   1609:           } else {
                   1610:             period = period*3/2;
                   1611:             count = 0;
                   1612:           }
                   1613:         }
                   1614:       }
                   1615:   }
                   1616: }
                   1617:
                   1618: /* for inv_or_split_dalg */
                   1619:
                   1620: int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT *nmmat,Z *dn,int **rindp,int **cindp)
                   1621: {
                   1622:   MAT bmat,xmat;
                   1623:   Z **a0,**a,**b,**x,**nm;
                   1624:   Z *ai,*bi,*xi;
                   1625:   int row,col;
                   1626:   int **w;
                   1627:   int *wi;
                   1628:   int **wc;
                   1629:   Z mdq,q,s,u;
                   1630:   Z tn;
                   1631:   int ind,md,i,j,k,l,li,ri,rank;
                   1632:   unsigned int t;
                   1633:   int *cinfo,*rinfo;
                   1634:   int *rind,*cind;
                   1635:   int count;
                   1636:   int ret;
                   1637:   struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
                   1638:   int period;
                   1639:   int *wx,*ptr;
                   1640:   int wxsize,nsize;
                   1641:   Z wn;
                   1642:   Z wq;
                   1643:   DP m;
                   1644:
                   1645:   a0 = (Z **)mat->body;
                   1646:   row = mat->row; col = mat->col;
                   1647:   w = (int **)almat(row,col);
                   1648:   for ( ind = 0; ; ind++ ) {
                   1649:     md = get_lprime(ind);
1.5       noro     1650:     STOZ(md,mdq);
1.1       noro     1651:     for ( i = 0; i < row; i++ )
                   1652:       for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
                   1653:         wi[j] = remqi((Q)ai[j],md);
                   1654:
                   1655:     if ( DP_Print > 3 ) {
                   1656:       fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
                   1657:     }
                   1658:     rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
                   1659:     if ( DP_Print > 3 ) {
                   1660:       fprintf(asir_out,"done.\n"); fflush(asir_out);
                   1661:     }
                   1662:
                   1663:     /* this part is added for inv_or_split_dalg */
                   1664:     for ( i = 0; i < col-1; i++ ) {
                   1665:       if ( !cinfo[i] ) {
                   1666:         m = mb[i];
                   1667:         for ( j = i+1; j < col-1; j++ )
                   1668:           if ( dp_redble(mb[j],m) )
                   1669:             cinfo[j] = -1;
                   1670:       }
                   1671:     }
                   1672:
                   1673:     a = (Z **)almat_pointer(rank,rank); /* lhs mat */
                   1674:     MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */
                   1675:     for ( j = li = ri = 0; j < col; j++ )
1.4       noro     1676:       if ( cinfo[j] > 0 ) {
1.1       noro     1677:         /* the column is in lhs */
                   1678:         for ( i = 0; i < rank; i++ ) {
                   1679:           w[i][li] = w[i][j];
                   1680:           a[i][li] = a0[rinfo[i]][j];
                   1681:         }
                   1682:         li++;
1.4       noro     1683:       } else if ( !cinfo[j] ) {
1.1       noro     1684:         /* the column is in rhs */
                   1685:         for ( i = 0; i < rank; i++ )
                   1686:           b[i][ri] = a0[rinfo[i]][j];
                   1687:         ri++;
                   1688:       }
                   1689:
                   1690:       /* solve Ax=B; A: rank x rank, B: rank x ri */
                   1691:       /* algorithm
                   1692:          c <- B
                   1693:          x <- 0
                   1694:          q <- 1
                   1695:          do
                   1696:            t <- A^(-1)c mod p
                   1697:            x <- x+qt
                   1698:            c <- (c-At)/p
                   1699:            q <- qp
                   1700:          end do
                   1701:          then Ax-B=0 mod q and b=(B-Ax)/q hold after "do".
                   1702:       */
                   1703:       MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;
                   1704:       MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;
                   1705:       wc = (int **)almat(rank,ri);
                   1706:       *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                   1707:       *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
                   1708:
                   1709:       period = F4_INTRAT_PERIOD;
                   1710:       for ( q = ONE, count = 0; ; ) {
                   1711:         if ( DP_Print > 3 )
                   1712:           fprintf(stderr,"o");
                   1713:         /* wc = b mod md */
                   1714:         for ( i = 0; i < rank; i++ )
1.3       noro     1715:           for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
1.1       noro     1716:             wi[j] = remqi((Q)bi[j],md);
                   1717:         /* wc = A^(-1)wc; wc is normalized */
                   1718:         solve_by_lu_mod(w,rank,md,wc,ri,1);
                   1719:         /* x += q*wc */
                   1720:         for ( i = 0; i < rank; i++ )
                   1721:           for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);
1.3       noro     1722:         /* b =(b-A*wc)/md */
1.1       noro     1723:         for ( i = 0; i < rank; i++ )
                   1724:           for ( j = 0; j < ri; j++ ) {
1.3       noro     1725:             mpz_t uz;
                   1726:
                   1727:             if ( b[i][j] )
                   1728:               mpz_init_set(uz,BDY(b[i][j]));
                   1729:             else
                   1730:               mpz_init_set_ui(uz,0);
                   1731:             for ( k = 0; k < rank; k++ ) {
                   1732:               if ( a[i][k] && wc[k][j] ) {
                   1733:                 if ( wc[k][j] < 0 )
                   1734:                   mpz_addmul_ui(uz,BDY(a[i][k]),-wc[k][j]);
                   1735:                 else
                   1736:                   mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]);
                   1737:               }
                   1738:             }
                   1739:             MPZTOZ(uz,u);
1.1       noro     1740:             divsz(u,mdq,&b[i][j]);
                   1741:           }
                   1742:         count++;
                   1743:         /* q = q*md */
                   1744:         mulz(q,mdq,&u); q = u;
                   1745:         if ( count == period ) {
                   1746:           ret = intmtoratm(xmat,q,*nmmat,dn);
                   1747:           if ( ret ) {
                   1748:             for ( j = k = l = 0; j < col; j++ )
                   1749:               if ( cinfo[j] > 0 )
                   1750:                 rind[k++] = j;
                   1751:               else if ( !cinfo[j] )
                   1752:                 cind[l++] = j;
                   1753:             ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
                   1754:             if ( ret ) {
                   1755:               *rindp = rind;
                   1756:               *cindp = cind;
                   1757:               for ( j = k = 0; j < col; j++ )
                   1758:                 if ( !cinfo[j] )
                   1759:                   cind[k++] = j;
                   1760:               return rank;
                   1761:             }
                   1762:           } else {
                   1763:             period = period*3/2;
                   1764:             count = 0;
                   1765:           }
                   1766:         }
                   1767:       }
                   1768:   }
                   1769: }
1.6       noro     1770:
                   1771: #if SIZEOF_LONG == 8
                   1772: mp_limb_t remqi64(Q a,mp_limb_t mod)
                   1773: {
                   1774:   mp_limb_t c,nm,dn;
                   1775:   mpz_t r;
                   1776:
                   1777:   if ( !a ) return 0;
                   1778:   else if ( a->z ) {
                   1779:     mpz_init(r);
                   1780:     c = mpz_fdiv_r_ui(r,BDY((Z)a),mod);
                   1781:   } else {
                   1782:     mpz_init(r);
                   1783:     nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod);
                   1784:     dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod);
                   1785:     dn = invmod64(dn,mod);
                   1786:     c = mulmod64(nm,dn,mod);
                   1787:   }
                   1788:   return c;
                   1789: }
                   1790:
                   1791: int generic_gauss_elim_mod64(mp_limb_t **mat,int row,int col,mp_limb_t md,int *colstat);
                   1792: mp_limb_t get_lprime64(int ind);
                   1793:
                   1794: int generic_gauss_elim64(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
                   1795: {
                   1796:   mp_limb_t **wmat;
                   1797:   mp_limb_t *wmi;
                   1798:   mp_limb_t md,inv,t,t1;
                   1799:   Z **bmat,**tmat,*bmi,*tmi;
                   1800:   Z q,m1,m2,m3,s,u;
                   1801:   int *colstat,*wcolstat,*rind,*cind;
                   1802:   int row,col,ind,i,j,k,l,rank,rank0;
                   1803:   MAT r,crmat;
                   1804:   int ret;
                   1805:
                   1806:   bmat = (Z **)mat->body;
                   1807:   row = mat->row; col = mat->col;
                   1808:   wmat = (mp_limb_t **)almat64(row,col);
                   1809:   colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                   1810:   wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                   1811:   for ( ind = 0; ; ind++ ) {
                   1812:     if ( DP_Print ) {
                   1813:       fprintf(asir_out,"."); fflush(asir_out);
                   1814:     }
                   1815:     md = get_lprime64(ind);
                   1816:     for ( i = 0; i < row; i++ )
                   1817:       for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
                   1818:         wmi[j] = remqi64((Q)bmi[j],md);
                   1819:     rank = generic_gauss_elim_mod64(wmat,row,col,md,wcolstat);
                   1820:     if ( !ind ) {
                   1821: RESET:
                   1822:       UTOZ(md,m1);
                   1823:       rank0 = rank;
                   1824:       bcopy(wcolstat,colstat,col*sizeof(int));
                   1825:       MKMAT(crmat,rank,col-rank);
                   1826:       MKMAT(r,rank,col-rank); *nm = r;
                   1827:       tmat = (Z **)crmat->body;
                   1828:       for ( i = 0; i < rank; i++ )
                   1829:         for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                   1830:           if ( !colstat[j] ) { UTOZ(wmi[j],tmi[k]); k++; }
                   1831:     } else {
                   1832:       if ( rank < rank0 ) {
                   1833:         if ( DP_Print ) {
                   1834:           fprintf(asir_out,"lower rank matrix; continuing...\n");
                   1835:           fflush(asir_out);
                   1836:         }
                   1837:         continue;
                   1838:       } else if ( rank > rank0 ) {
                   1839:         if ( DP_Print ) {
                   1840:           fprintf(asir_out,"higher rank matrix; resetting...\n");
                   1841:           fflush(asir_out);
                   1842:         }
                   1843:         goto RESET;
                   1844:       } else {
                   1845:         for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
                   1846:         if ( j < col ) {
                   1847:           if ( DP_Print ) {
                   1848:             fprintf(asir_out,"inconsitent colstat; resetting...\n");
                   1849:             fflush(asir_out);
                   1850:           }
                   1851:           goto RESET;
                   1852:         }
                   1853:       }
                   1854:
                   1855:       inv = invmod64(remqi64((Q)m1,md),md);
                   1856:       UTOZ(md,m2); mulz(m1,m2,&m3);
                   1857:       for ( i = 0; i < rank; i++ )
                   1858:         for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                   1859:           if ( !colstat[j] ) {
                   1860:             if ( tmi[k] ) {
                   1861:             /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                   1862:               t = remqi64((Q)tmi[k],md);
                   1863:               if ( wmi[j] >= t ) t = wmi[j]-t;
                   1864:               else t = md-(t-wmi[j]);
                   1865:               t1 = mulmod64(t,inv,md);
                   1866:               UTOZ(t1,u); mulz(m1,u,&s);
                   1867:               addz(tmi[k],s,&u); tmi[k] = u;
                   1868:             } else if ( wmi[j] ) {
                   1869:             /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                   1870:               t = mulmod64(wmi[j],inv,md);
                   1871:               UTOZ(t,u); mulz(m1,u,&s); tmi[k] = s;
                   1872:             }
                   1873:             k++;
                   1874:           }
                   1875:       m1 = m3;
                   1876:       if ( ind % F4_INTRAT_PERIOD )
                   1877:         ret = 0;
                   1878:       else
                   1879:         ret = intmtoratm(crmat,m1,*nm,dn);
                   1880:       if ( ret ) {
                   1881:         *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                   1882:         *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
                   1883:         for ( j = k = l = 0; j < col; j++ )
                   1884:           if ( colstat[j] ) rind[k++] = j;
                   1885:           else cind[l++] = j;
                   1886:         if ( gensolve_check(mat,*nm,*dn,rind,cind) )
                   1887:           return rank;
                   1888:       }
                   1889:     }
                   1890:   }
                   1891: }
                   1892: #endif

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