[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.21

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

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