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

Annotation of OpenXM_contrib2/asir2018/builtin/array.c, Revision 1.3

1.1       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
                     26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
                     27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.3     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2018/builtin/array.c,v 1.2 2018/09/28 08:20:27 noro Exp $
1.1       noro       49: */
                     50: #include "ca.h"
                     51: #include "base.h"
                     52: #include "parse.h"
                     53: #include "inline.h"
                     54:
                     55: #include <sys/types.h>
                     56: #include <sys/stat.h>
                     57: #if !defined(_MSC_VER)
                     58: #include <unistd.h>
                     59: #endif
                     60:
                     61: #define F4_INTRAT_PERIOD 8
                     62:
                     63: #if 0
                     64: #undef DMAR
                     65: #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d);
                     66: #endif
                     67:
                     68: extern int DP_Print; /* XXX */
                     69:
                     70:
                     71: void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(), Ptriangleq();
                     72: void Pinvmat();
                     73: void Pnewbytearray(),Pmemoryplot_to_coord();
                     74:
                     75: void Pgeneric_gauss_elim();
                     76: void Pgeneric_gauss_elim_mod();
                     77:
                     78: void Pindep_rows_mod();
                     79:
                     80: void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();
                     81: void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(), Pltov();
                     82: void Pgeninv_sf_swap();
                     83: void sepvect();
                     84: void Pmulmat_gf2n();
                     85: void Pbconvmat_gf2n();
                     86: void Pmul_vect_mat_gf2n();
                     87: void PNBmul_gf2n();
                     88: void Pmul_mat_vect_int();
                     89: void Psepmat_destructive();
                     90: void Px962_irredpoly_up2();
                     91: void Pirredpoly_up2();
                     92: void Pnbpoly_up2();
                     93: void Pqsort();
                     94: void Pexponent_vector();
                     95: void Pmat_swap_row_destructive();
                     96: void Pmat_swap_col_destructive();
                     97: void Pvect();
                     98: void Pmat();
                     99: void Pmatc();
                    100: void Pnd_det();
                    101: void Plu_mat();
                    102: void Pmat_col();
                    103: void Plusolve_prep();
                    104: void Plusolve_main();
                    105:
                    106: struct ftab array_tab[] = {
                    107: #if 0
                    108:   {"lu_mat",Plu_mat,1},
                    109: #endif
                    110:   {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},
                    111:   {"lu_gfmmat",Plu_gfmmat,2},
                    112:   {"mat_to_gfmmat",Pmat_to_gfmmat,2},
                    113:   {"generic_gauss_elim",Pgeneric_gauss_elim,1},
                    114:   {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},
                    115:   {"indep_rows_mod",Pindep_rows_mod,2},
                    116:   {"newvect",Pnewvect,-2},
                    117:   {"vect",Pvect,-99999999},
                    118:   {"vector",Pnewvect,-2},
                    119:   {"exponent_vector",Pexponent_vector,-99999999},
                    120:   {"newmat",Pnewmat,-3},
                    121:   {"matrix",Pnewmat,-3},
                    122:   {"mat",Pmat,-99999999},
                    123:   {"matr",Pmat,-99999999},
                    124:   {"matc",Pmatc,-99999999},
                    125:   {"newbytearray",Pnewbytearray,-2},
                    126:   {"memoryplot_to_coord",Pmemoryplot_to_coord,1},
                    127:   {"sepmat_destructive",Psepmat_destructive,2},
                    128:   {"sepvect",Psepvect,2},
                    129:   {"qsort",Pqsort,-2},
                    130:   {"vtol",Pvtol,1},
                    131:   {"ltov",Pltov,1},
                    132:   {"size",Psize,1},
                    133:   {"det",Pdet,-2},
                    134:   {"nd_det",Pnd_det,-2},
                    135:   {"invmat",Pinvmat,-2},
                    136:   {"leqm",Pleqm,2},
                    137:   {"leqm1",Pleqm1,2},
                    138:   {"geninvm",Pgeninvm,2},
                    139:   {"geninvm_swap",Pgeninvm_swap,2},
                    140:   {"geninv_sf_swap",Pgeninv_sf_swap,1},
                    141:   {"remainder",Premainder,2},
                    142:   {"sremainder",Psremainder,2},
                    143:   {"mulmat_gf2n",Pmulmat_gf2n,1},
                    144:   {"bconvmat_gf2n",Pbconvmat_gf2n,-4},
                    145:   {"mul_vect_mat_gf2n",Pmul_vect_mat_gf2n,2},
                    146:   {"mul_mat_vect_int",Pmul_mat_vect_int,2},
                    147:   {"nbmul_gf2n",PNBmul_gf2n,3},
                    148:   {"x962_irredpoly_up2",Px962_irredpoly_up2,2},
                    149:   {"irredpoly_up2",Pirredpoly_up2,2},
                    150:   {"nbpoly_up2",Pnbpoly_up2,2},
                    151:   {"mat_swap_row_destructive",Pmat_swap_row_destructive,3},
                    152:   {"mat_swap_col_destructive",Pmat_swap_col_destructive,3},
                    153:   {"mat_col",Pmat_col,2},
                    154:   {"lusolve_prep",Plusolve_prep,1},
                    155:   {"lusolve_main",Plusolve_main,1},
                    156:   {"triangleq",Ptriangleq,1},
                    157:   {0,0,0},
                    158: };
                    159:
                    160: typedef struct _ent { int j; unsigned int e; } ent;
                    161:
                    162: ent *get_row(FILE *,int *l);
                    163: void put_row(FILE *out,int l,ent *a);
                    164: void lu_elim(int *l,ent **a,int k,int i,int mul,int mod);
                    165: void lu_append(int *,ent **,int *,int,int,int);
                    166: void solve_l(int *,ent **,int,int *,int);
                    167: void solve_u(int *,ent **,int,int *,int);
                    168:
                    169:
                    170: static int *ul,*ll;
                    171: static ent **u,**l;
                    172: static int modulus;
                    173:
                    174: void Plusolve_prep(NODE arg,Q *rp)
                    175: {
                    176:   char *fname;
                    177:   FILE *in;
                    178:   int len,i,rank;
                    179:   int *rhs;
                    180:
                    181:   fname = BDY((STRING)ARG0(arg));
                    182:   in = fopen(fname,"r");
                    183:   modulus = getw(in);
                    184:   len = getw(in);
                    185:   ul = (int *)MALLOC_ATOMIC(len*sizeof(int));
                    186:   u = (ent **)MALLOC(len*sizeof(ent *));
                    187:   ll = (int *)MALLOC_ATOMIC(len*sizeof(int));
                    188:   l = (ent **)MALLOC(len*sizeof(ent *));
                    189:   for ( i = 0; i < len; i++ ) {
                    190:     u[i] = get_row(in,&ul[i]);
                    191:   }
                    192:   for ( i = 0; i < len; i++ ) {
                    193:     l[i] = get_row(in,&ll[i]);
                    194:   }
                    195:   fclose(in);
                    196:   *rp = (Q)ONE;
                    197: }
                    198:
                    199: void Plusolve_main(NODE arg,VECT *rp)
                    200: {
                    201:   Z *d,*p;
                    202:   VECT v,r;
                    203:   int len,i;
                    204:   int *rhs;
                    205:
                    206:   v = (VECT)ARG0(arg); len = v->len;
                    207:   d = (Z *)BDY(v);
                    208:   rhs = (int *)MALLOC_ATOMIC(len*sizeof(int));
1.2       noro      209:   for ( i = 0; i < len; i++ ) rhs[i] = ZTOS(d[i]);
1.1       noro      210:   solve_l(ll,l,len,rhs,modulus);
                    211:   solve_u(ul,u,len,rhs,modulus);
                    212:   NEWVECT(r); r->len = len;
                    213:   r->body = (pointer *)MALLOC(len*sizeof(pointer));
                    214:   p = (Z *)r->body;
                    215:   for ( i = 0; i < len; i++ )
1.2       noro      216:     STOZ(rhs[i],p[i]);
1.1       noro      217:   *rp = r;
                    218: }
                    219:
                    220: ent *get_row(FILE *in,int *l)
                    221: {
                    222:   int len,i;
                    223:   ent *a;
                    224:
                    225:   *l = len = getw(in);
                    226:   a = (ent *)MALLOC_ATOMIC(len*sizeof(ent));
                    227:   for ( i = 0; i < len; i++ ) {
                    228:     a[i].j = getw(in);
                    229:     a[i].e = getw(in);
                    230:   }
                    231:   return a;
                    232: }
                    233:
                    234: void lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int mod)
                    235: {
                    236:   int i,j,k,s,mul;
                    237:   unsigned int inv;
                    238:   int *ll2;
                    239:
                    240:   ll2 = (int *)MALLOC_ATOMIC(n*sizeof(int));
                    241:   for ( i = 0; i < n; i++ ) ll2[i] = 0;
                    242:   for ( i = 0; i < n; i++ ) {
                    243:     fprintf(stderr,"i=%d\n",i);
                    244:     inv = invm(u[i][0].e,mod);
                    245:     for ( k = i+1; k < n; k++ )
                    246:       if ( u[k][0].j == n-i ) {
                    247:         s = u[k][0].e;
                    248:         DMAR(s,inv,0,mod,mul);
                    249:         lu_elim(ul,u,k,i,mul,mod);
                    250:         lu_append(ll,l,ll2,k,i,mul);
                    251:       }
                    252:   }
                    253: }
                    254:
                    255: #define INITLEN 10
                    256:
                    257: void lu_append(int *l,ent **a,int *l2,int k,int i,int mul)
                    258: {
                    259:   int len;
                    260:   ent *p;
                    261:
                    262:   len = l[k];
                    263:   if ( !len ) {
                    264:     a[k] = p = (ent *)MALLOC_ATOMIC(INITLEN*sizeof(ent));
                    265:     p[0].j = i; p[0].e = mul;
                    266:     l[k] = 1; l2[k] = INITLEN;
                    267:   } else {
                    268:     if ( l2[k] == l[k] ) {
                    269:       l2[k] *= 2;
                    270:       a[k] = REALLOC(a[k],l2[k]*sizeof(ent));
                    271:     }
                    272:     p =a[k];
                    273:     p[l[k]].j = i; p[l[k]].e = mul;
                    274:     l[k]++;
                    275:   }
                    276: }
                    277:
                    278: /* a[k] = a[k]-mul*a[i] */
                    279:
                    280: void lu_elim(int *l,ent **a,int k,int i,int mul,int mod)
                    281: {
                    282:   ent *ak,*ai,*w;
                    283:   int lk,li,j,m,p,q,r,s,t,j0;
                    284:
                    285:   ak = a[k]; ai = a[i]; lk = l[k]; li = l[i];
                    286:   w = (ent *)alloca((lk+li)*sizeof(ent));
                    287:   p = 0; q = 0; j = 0;
                    288:   mul = mod-mul;
                    289:   while ( p < lk && q < li ) {
                    290:     if ( ak[p].j > ai[q].j ) {
                    291:       w[j] = ak[p]; j++; p++;
                    292:     } else if ( ak[p].j < ai[q].j ) {
                    293:       w[j].j = ai[q].j;
                    294:       t = ai[q].e;
                    295:       DMAR(t,mul,0,mod,r);
                    296:       w[j].e = r;
                    297:       j++; q++;
                    298:     } else {
                    299:       t = ai[q].e; s = ak[p].e;
                    300:       DMAR(t,mul,s,mod,r);
                    301:       if ( r ) {
                    302:         w[j].j = ai[q].j; w[j].e = r; j++;
                    303:       }
                    304:       p++; q++;
                    305:     }
                    306:   }
                    307:   if ( q == li )
                    308:     while ( p < lk ) {
                    309:       w[j] = ak[p]; j++; p++;
                    310:     }
                    311:   else if ( p == lk )
                    312:     while ( q < li ) {
                    313:       w[j].j = ai[q].j;
                    314:       t = ai[q].e;
                    315:       DMAR(t,mul,0,mod,r);
                    316:       w[j].e = r;
                    317:       j++; q++;
                    318:     }
                    319:   if ( j <= lk ) {
                    320:     for ( m = 0; m < j; m++ ) ak[m] = w[m];
                    321:   } else {
                    322:     a[k] = ak = (ent *)MALLOC_ATOMIC(j*sizeof(ent));
                    323:     for ( m = 0; m < j; m++ ) ak[m] = w[m];
                    324:   }
                    325:   l[k] = j;
                    326: }
                    327:
                    328: void solve_l(int *ll,ent **l,int n,int *rhs,int mod)
                    329: {
                    330:   int j,k,s,len;
                    331:   ent *p;
                    332:
                    333:   for ( j = 0; j < n; j++ ) {
                    334:     len = ll[j]; p = l[j];
                    335:     for ( k = 0, s = 0; k < len; k++ )
                    336:       s = dmar(p[k].e,rhs[p[k].j],s,mod);
                    337:     rhs[j] -=  s;
                    338:     if ( rhs[j] < 0 ) rhs[j] += mod;
                    339:   }
                    340: }
                    341:
                    342: void solve_u(int *ul,ent **u,int n,int *rhs,int mod)
                    343: {
                    344:   int j,k,s,len,inv;
                    345:   ent *p;
                    346:
                    347:   for ( j = n-1; j >= 0; j-- ) {
                    348:     len = ul[j]; p = u[j];
                    349:     for ( k = 1, s = 0; k < len; k++ )
                    350:       s = dmar(p[k].e,rhs[p[k].j],s,mod);
                    351:     rhs[j] -=  s;
                    352:     if ( rhs[j] < 0 ) rhs[j] += mod;
                    353:     inv = invm((unsigned int)p[0].e,mod);
                    354:     rhs[j] = dmar(rhs[j],inv,0,mod);
                    355:   }
                    356: }
                    357:
                    358: int comp_obj(Obj *a,Obj *b)
                    359: {
                    360:   return arf_comp(CO,*a,*b);
                    361: }
                    362:
                    363: static FUNC generic_comp_obj_func;
                    364: static NODE generic_comp_obj_arg;
                    365: static NODE generic_comp_obj_option;
                    366:
                    367: int generic_comp_obj(Obj *a,Obj *b)
                    368: {
                    369:   Q r;
                    370:
                    371:   BDY(generic_comp_obj_arg)=(pointer)(*a);
                    372:   BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b);
                    373:   r = (Q)bevalf_with_opts(generic_comp_obj_func,generic_comp_obj_arg,generic_comp_obj_option);
                    374:   if ( !r )
                    375:     return 0;
                    376:   else
                    377:     return sgnq(r)>0?1:-1;
                    378: }
                    379:
                    380:
                    381: void Pqsort(NODE arg,LIST *rp)
                    382: {
                    383:   VECT vect;
                    384:   NODE n,n1;
                    385:   P p;
                    386:   V v;
                    387:   FUNC func;
                    388:   int len,i;
                    389:   pointer *a;
                    390:   Obj t;
                    391:
                    392:   t = ARG0(arg);
                    393:     if (OID(t) == O_LIST) {
                    394:         n = (NODE)BDY((LIST)t);
                    395:         len = length(n);
                    396:         MKVECT(vect,len);
                    397:         for ( i = 0; i < len; i++, n = NEXT(n) ) {
                    398:             BDY(vect)[i] = BDY(n);
                    399:         }
                    400:
                    401:     }else if (OID(t) != O_VECT) {
                    402:         error("qsort : invalid argument");
                    403:     }else {
                    404:         vect = (VECT)t;
                    405:     }
                    406:   if ( argc(arg) == 1 )
                    407:     qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj);
                    408:   else {
                    409:     p = (P)ARG1(arg);
                    410:     if ( !p || OID(p)!=2 )
                    411:       error("qsort : invalid argument");
                    412:     v = VR(p);
                    413:     gen_searchf(NAME(v),&func);
                    414:     if ( !func ) {
                    415:       if ( (long)v->attr != V_SR )
                    416:         error("qsort : no such function");
                    417:       func = (FUNC)v->priv;
                    418:     }
                    419:     generic_comp_obj_func = func;
                    420:     MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);
                    421:     generic_comp_obj_option = current_option;
                    422:     qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);
                    423:   }
                    424:     if (OID(t) == O_LIST) {
                    425:         a = BDY(vect);
                    426:         for ( i = len - 1, n = 0; i >= 0; i-- ) {
                    427:             MKNODE(n1,a[i],n); n = n1;
                    428:         }
                    429:         MKLIST(*rp,n);
                    430:     }else {
                    431:         *rp = (LIST)vect;
                    432:     }
                    433: }
                    434:
                    435: void PNBmul_gf2n(NODE arg,GF2N *rp)
                    436: {
                    437:   GF2N a,b;
                    438:   GF2MAT mat;
                    439:   int n,w;
                    440:   unsigned int *ab,*bb;
                    441:   UP2 r;
                    442:
                    443:   a = (GF2N)ARG0(arg);
                    444:   b = (GF2N)ARG1(arg);
                    445:   mat = (GF2MAT)ARG2(arg);
                    446:   if ( !a || !b )
                    447:     *rp = 0;
                    448:   else {
                    449:     n = mat->row;
                    450:     w = (n+BSH-1)/BSH;
                    451:
                    452:     ab = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
                    453:     bzero((char *)ab,w*sizeof(unsigned int));
                    454:     bcopy(a->body->b,ab,(a->body->w)*sizeof(unsigned int));
                    455:
                    456:     bb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
                    457:     bzero((char *)bb,w*sizeof(unsigned int));
                    458:     bcopy(b->body->b,bb,(b->body->w)*sizeof(unsigned int));
                    459:
                    460:     NEWUP2(r,w);
                    461:     bzero((char *)r->b,w*sizeof(unsigned int));
                    462:     mul_nb(mat,ab,bb,r->b);
                    463:     r->w = w;
                    464:     _adjup2(r);
                    465:     if ( !r->w )
                    466:       *rp = 0;
                    467:     else
                    468:       MKGF2N(r,*rp);
                    469:   }
                    470: }
                    471:
                    472: void Pmul_vect_mat_gf2n(NODE arg,GF2N *rp)
                    473: {
                    474:   GF2N a;
                    475:   GF2MAT mat;
                    476:   int n,w;
                    477:   unsigned int *b;
                    478:   UP2 r;
                    479:
                    480:   a = (GF2N)ARG0(arg);
                    481:   mat = (GF2MAT)ARG1(arg);
                    482:   if ( !a )
                    483:     *rp = 0;
                    484:   else {
                    485:     n = mat->row;
                    486:     w = (n+BSH-1)/BSH;
                    487:     b = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
                    488:     bzero((char *)b,w*sizeof(unsigned int));
                    489:     bcopy(a->body->b,b,(a->body->w)*sizeof(unsigned int));
                    490:     NEWUP2(r,w);
                    491:     bzero((char *)r->b,w*sizeof(unsigned int));
                    492:     mulgf2vectmat(mat->row,b,mat->body,r->b);
                    493:     r->w = w;
                    494:     _adjup2(r);
                    495:     if ( !r->w )
                    496:       *rp = 0;
                    497:     else {
                    498:       MKGF2N(r,*rp);
                    499:     }
                    500:   }
                    501: }
                    502:
                    503: void Pbconvmat_gf2n(NODE arg,LIST *rp)
                    504: {
                    505:   P p0,p1;
                    506:   int to;
                    507:   GF2MAT p01,p10;
                    508:   GF2N root;
                    509:   NODE n0,n1;
                    510:
                    511:   p0 = (P)ARG0(arg);
                    512:   p1 = (P)ARG1(arg);
                    513:   to = ARG2(arg)?1:0;
                    514:   if ( argc(arg) == 4 ) {
                    515:     root = (GF2N)ARG3(arg);
                    516:     compute_change_of_basis_matrix_with_root(p0,p1,to,root,&p01,&p10);
                    517:   } else
                    518:     compute_change_of_basis_matrix(p0,p1,to,&p01,&p10);
                    519:   MKNODE(n1,p10,0); MKNODE(n0,p01,n1);
                    520:   MKLIST(*rp,n0);
                    521: }
                    522:
                    523: void Pmulmat_gf2n(NODE arg,GF2MAT *rp)
                    524: {
                    525:   GF2MAT m;
                    526:
                    527:   if ( !compute_multiplication_matrix((P)ARG0(arg),&m) )
                    528:     error("mulmat_gf2n : input is not a normal polynomial");
                    529:   *rp = m;
                    530: }
                    531:
                    532: void Psepmat_destructive(NODE arg,LIST *rp)
                    533: {
                    534:   MAT mat,mat1;
                    535:   int i,j,row,col;
                    536:   Z **a,**a1;
                    537:   Z ent;
                    538:   Z nm,mod,rem,quo;
                    539:   int sgn;
                    540:   NODE n0,n1;
                    541:
                    542:   mat = (MAT)ARG0(arg); mod = (Z)ARG1(arg);
                    543:   row = mat->row; col = mat->col;
                    544:   MKMAT(mat1,row,col);
                    545:   a = (Z **)mat->body; a1 = (Z **)mat1->body;
                    546:   for ( i = 0; i < row; i++ )
                    547:     for ( j = 0; j < col; j++ ) {
                    548:       ent = a[i][j];
                    549:       if ( !ent )
                    550:         continue;
                    551:       sgn = sgnz(ent);
                    552:       absz(ent,&nm);
                    553:       divqrz(nm,mod,&quo,&rem);
                    554:       if ( sgn > 0 ) {
                    555:         a[i][j] = rem;
                    556:         a1[i][j] = quo;
                    557:       } else {
                    558:         chsgnz(rem,&a[i][j]);
                    559:         chsgnz(quo,&a1[i][j]);
                    560:       }
                    561:     }
                    562:   MKNODE(n1,mat1,0); MKNODE(n0,mat,n1);
                    563:   MKLIST(*rp,n0);
                    564: }
                    565:
                    566: void Psepvect(NODE arg,VECT *rp)
                    567: {
1.2       noro      568:   sepvect((VECT)ARG0(arg),ZTOS((Q)ARG1(arg)),rp);
1.1       noro      569: }
                    570:
                    571: void sepvect(VECT v,int d,VECT *rp)
                    572: {
                    573:   int i,j,k,n,q,q1,r;
                    574:   pointer *pv,*pw,*pu;
                    575:   VECT w,u;
                    576:
                    577:   n = v->len;
                    578:   if ( d > n )
                    579:     d = n;
                    580:   q = n/d; r = n%d; q1 = q+1;
                    581:   MKVECT(w,d); *rp = w;
                    582:   pv = BDY(v); pw = BDY(w); k = 0;
                    583:   for ( i = 0; i < r; i++ ) {
                    584:     MKVECT(u,q1); pw[i] = (pointer)u;
                    585:     for ( pu = BDY(u), j = 0; j < q1; j++, k++ )
                    586:       pu[j] = pv[k];
                    587:   }
                    588:   for ( ; i < d; i++ ) {
                    589:     MKVECT(u,q); pw[i] = (pointer)u;
                    590:     for ( pu = BDY(u), j = 0; j < q; j++, k++ )
                    591:       pu[j] = pv[k];
                    592:   }
                    593: }
                    594:
                    595: void Pnewvect(NODE arg,VECT *rp)
                    596: {
                    597:   int len,i,r;
                    598:   VECT vect;
                    599:   pointer *vb;
                    600:   LIST list;
                    601:   NODE tn;
                    602:
                    603:   asir_assert(ARG0(arg),O_N,"newvect");
1.2       noro      604:   len = ZTOS((Q)ARG0(arg));
1.1       noro      605:   if ( len < 0 )
                    606:     error("newvect : invalid size");
                    607:   MKVECT(vect,len);
                    608:   if ( argc(arg) == 2 ) {
                    609:     list = (LIST)ARG1(arg);
                    610:     asir_assert(list,O_LIST,"newvect");
                    611: #if 0
                    612:     for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
                    613:     if ( r > len ) {
                    614:       *rp = vect;
                    615:       return;
                    616:     }
                    617: #endif
                    618:     for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )
                    619:       vb[i] = (pointer)BDY(tn);
                    620:   }
                    621:   *rp = vect;
                    622: }
                    623:
                    624: void Pvect(NODE arg,VECT *rp) {
                    625:   int len,i;
                    626:   VECT vect;
                    627:   pointer *vb;
                    628:   NODE tn;
                    629:
                    630:   if ( !arg ) {
                    631:     *rp =0;
                    632:     return;
                    633:   }
                    634:
                    635:   for (len = 0, tn = arg; tn; tn = NEXT(tn), len++);
                    636:   if ( len == 1 ) {
                    637:     if ( ARG0(arg) != 0 ) {
                    638:       switch ( OID(ARG0(arg)) ) {
                    639:         case O_VECT:
                    640:           *rp = ARG0(arg);
                    641:           return;
                    642:         case O_LIST:
                    643:           for ( len = 0, tn = ARG0(arg); tn; tn = NEXT(tn), len++ );
                    644:           MKVECT(vect,len-1);
                    645:           for ( i = 0, tn = BDY((LIST)ARG0(arg)), vb =BDY(vect);
                    646:               tn; i++, tn = NEXT(tn) )
                    647:             vb[i] = (pointer)BDY(tn);
                    648:           *rp=vect;
                    649:           return;
                    650:       }
                    651:     }
                    652:   }
                    653:   MKVECT(vect,len);
                    654:   for ( i = 0, tn = arg, vb = BDY(vect); tn; i++, tn = NEXT(tn) )
                    655:     vb[i] = (pointer)BDY(tn);
                    656:   *rp = vect;
                    657: }
                    658:
                    659: void Pexponent_vector(NODE arg,DP *rp)
                    660: {
                    661:   nodetod(arg,rp);
                    662: }
                    663:
                    664: void Pnewbytearray(NODE arg,BYTEARRAY *rp)
                    665: {
                    666:   int len,i,r;
                    667:   BYTEARRAY array;
                    668:   unsigned char *vb;
                    669:   char *str;
                    670:   LIST list;
                    671:   NODE tn;
                    672:   int ac;
                    673:   struct stat sbuf;
                    674:   char *fname;
                    675:   FILE *fp;
                    676:
                    677:   ac = argc(arg);
                    678:   if ( ac == 1 ) {
                    679:     if ( !OID((Obj)ARG0(arg)) ) error("newbytearray : invalid argument");
                    680:     switch ( OID((Obj)ARG0(arg)) ) {
                    681:       case O_STR:
                    682:         fname = BDY((STRING)ARG0(arg));
                    683:         fp = fopen(fname,"rb");
                    684:         if ( !fp ) error("newbytearray : fopen failed");
                    685:         if ( stat(fname,&sbuf) < 0 )
                    686:           error("newbytearray : stat failed");
                    687:         len = sbuf.st_size;
                    688:         MKBYTEARRAY(array,len);
                    689:         fread(BDY(array),len,sizeof(char),fp);
                    690:         break;
                    691:       case O_N:
                    692:         if ( !RATN(ARG0(arg)) )
                    693:           error("newbytearray : invalid argument");
1.2       noro      694:         len = ZTOS((Q)ARG0(arg));
1.1       noro      695:         if ( len < 0 )
                    696:           error("newbytearray : invalid size");
                    697:         MKBYTEARRAY(array,len);
                    698:         break;
                    699:       default:
                    700:         error("newbytearray : invalid argument");
                    701:     }
                    702:   } else if ( ac == 2 ) {
                    703:     asir_assert(ARG0(arg),O_N,"newbytearray");
1.2       noro      704:     len = ZTOS((Q)ARG0(arg));
1.1       noro      705:     if ( len < 0 )
                    706:       error("newbytearray : invalid size");
                    707:     MKBYTEARRAY(array,len);
                    708:     if ( !ARG1(arg) )
                    709:       error("newbytearray : invalid initialization");
                    710:     switch ( OID((Obj)ARG1(arg)) ) {
                    711:       case O_LIST:
                    712:         list = (LIST)ARG1(arg);
                    713:         asir_assert(list,O_LIST,"newbytearray");
                    714:         for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
                    715:         if ( r <= len ) {
                    716:           for ( i = 0, tn = BDY(list), vb = BDY(array); tn;
                    717:             i++, tn = NEXT(tn) )
1.2       noro      718:             vb[i] = (unsigned char)ZTOS((Q)BDY(tn));
1.1       noro      719:         }
                    720:         break;
                    721:       case O_STR:
                    722:         str = BDY((STRING)ARG1(arg));
                    723:         r = strlen(str);
                    724:         if ( r <= len )
                    725:           bcopy(str,BDY(array),r);
                    726:         break;
                    727:       default:
                    728:         if ( !ARG1(arg) )
                    729:           error("newbytearray : invalid initialization");
                    730:     }
                    731:   } else
                    732:     error("newbytearray : invalid argument");
                    733:   *rp = array;
                    734: }
                    735:
                    736: #define MEMORY_GETPOINT(a,len,x,y) (((a)[(len)*(y)+((x)>>3)])&(1<<((x)&7)))
                    737:
                    738: void Pmemoryplot_to_coord(NODE arg,LIST *rp)
                    739: {
                    740:   int len,blen,y,i,j;
                    741:   unsigned char *a;
                    742:   NODE r0,r,n;
                    743:   LIST l;
                    744:   BYTEARRAY ba;
                    745:   Z iq,jq;
                    746:
                    747:   asir_assert(ARG0(arg),O_LIST,"memoryplot_to_coord");
                    748:   arg = BDY((LIST)ARG0(arg));
1.2       noro      749:   len = ZTOS((Q)ARG0(arg));
1.1       noro      750:   blen = (len+7)/8;
1.2       noro      751:   y = ZTOS((Q)ARG1(arg));
1.1       noro      752:   ba = (BYTEARRAY)ARG2(arg); a = ba->body;
                    753:   r0 = 0;
                    754:   for ( j = 0; j < y; j++ )
                    755:     for ( i = 0; i < len; i++ )
                    756:       if ( MEMORY_GETPOINT(a,blen,i,j) ) {
                    757:         NEXTNODE(r0,r);
1.2       noro      758:         STOZ(i,iq); STOZ(j,jq);
1.1       noro      759:         n = mknode(2,iq,jq);
                    760:         MKLIST(l,n);
                    761:         BDY(r) = l;
                    762:       }
                    763:   if ( r0 ) NEXT(r) = 0;
                    764:   MKLIST(*rp,r0);
                    765: }
                    766:
                    767: void Pnewmat(NODE arg,MAT *rp)
                    768: {
                    769:   int row,col;
                    770:   int i,j,r,c;
                    771:   NODE tn,sn;
                    772:   MAT m;
                    773:   pointer **mb;
                    774:   LIST list;
                    775:
                    776:   asir_assert(ARG0(arg),O_N,"newmat");
                    777:   asir_assert(ARG1(arg),O_N,"newmat");
1.2       noro      778:   row = ZTOS((Q)ARG0(arg)); col = ZTOS((Q)ARG1(arg));
1.1       noro      779:   if ( row < 0 || col < 0 )
                    780:     error("newmat : invalid size");
                    781:   MKMAT(m,row,col);
                    782:   if ( argc(arg) == 3 ) {
                    783:     list = (LIST)ARG2(arg);
                    784:     asir_assert(list,O_LIST,"newmat");
                    785:     for ( r = 0, c = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ) {
                    786:       for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) );
                    787:       c = MAX(c,j);
                    788:     }
                    789:     if ( (r > row) || (c > col) ) {
                    790:       *rp = m;
                    791:       return;
                    792:     }
                    793:     for ( i = 0, tn = BDY(list), mb = BDY(m); tn; i++, tn = NEXT(tn) ) {
                    794:       asir_assert(BDY(tn),O_LIST,"newmat");
                    795:       for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) )
                    796:         mb[i][j] = (pointer)BDY(sn);
                    797:     }
                    798:   }
                    799:   *rp = m;
                    800: }
                    801:
                    802: void Pmat(NODE arg, MAT *rp)
                    803: {
                    804:   int row,col;
                    805:   int i;
                    806:   MAT m;
                    807:   pointer **mb;
                    808:   pointer *ent;
                    809:   NODE tn, sn;
                    810:   VECT v;
                    811:
                    812:   if ( !arg ) {
                    813:     *rp =0;
                    814:     return;
                    815:   }
                    816:
                    817:   for (row = 0, tn = arg; tn; tn = NEXT(tn), row++);
                    818:   if ( row == 1 ) {
                    819:     if ( OID(ARG0(arg)) == O_MAT ) {
                    820:       *rp=ARG0(arg);
                    821:       return;
                    822:     } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {
                    823:       error("mat : invalid argument");
                    824:     }
                    825:   }
                    826:   if ( OID(ARG0(arg)) == O_VECT ) {
                    827:     v = ARG0(arg);
                    828:     col = v->len;
                    829:   } else if ( OID(ARG0(arg)) == O_LIST ) {
                    830:     for (col = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), col++);
                    831:   } else {
                    832:     error("mat : invalid argument");
                    833:   }
                    834:
                    835:   MKMAT(m,row,col);
                    836:   for (row = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), row++) {
                    837:     if ( BDY(tn) == 0 ) {
                    838:       error("mat : invalid argument");
                    839:     } else if ( OID(BDY(tn)) == O_VECT ) {
                    840:       v = tn->body;
                    841:       ent = BDY(v);
                    842:       for (i = 0; i < v->len; i++ ) mb[row][i] = (Obj)ent[i];
                    843:     } else if ( OID(BDY(tn)) == O_LIST ) {
                    844:       for (col = 0, sn = BDY((LIST)BDY(tn)); sn; col++, sn = NEXT(sn) )
                    845:         mb[row][col] = (pointer)BDY(sn);
                    846:     } else {
                    847:       error("mat : invalid argument");
                    848:     }
                    849:   }
                    850:   *rp = m;
                    851: }
                    852:
                    853: void Pmatc(NODE arg, MAT *rp)
                    854: {
                    855:   int row,col;
                    856:   int i;
                    857:   MAT m;
                    858:   pointer **mb;
                    859:   pointer *ent;
                    860:   NODE tn, sn;
                    861:   VECT v;
                    862:
                    863:   if ( !arg ) {
                    864:     *rp =0;
                    865:     return;
                    866:   }
                    867:
                    868:   for (col = 0, tn = arg; tn; tn = NEXT(tn), col++);
                    869:   if ( col == 1 ) {
                    870:     if ( OID(ARG0(arg)) == O_MAT ) {
                    871:       *rp=ARG0(arg);
                    872:       return;
                    873:     } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {
                    874:       error("matc : invalid argument");
                    875:     }
                    876:   }
                    877:   if ( OID(ARG0(arg)) == O_VECT ) {
                    878:     v = ARG0(arg);
                    879:     row = v->len;
                    880:   } else if ( OID(ARG0(arg)) == O_LIST ) {
                    881:     for (row = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), row++);
                    882:   } else {
                    883:     error("matc : invalid argument");
                    884:   }
                    885:
                    886:   MKMAT(m,row,col);
                    887:   for (col = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), col++) {
                    888:     if ( BDY(tn) == 0 ) {
                    889:       error("matc : invalid argument");
                    890:     } else if ( OID(BDY(tn)) == O_VECT ) {
                    891:       v = tn->body;
                    892:       ent = BDY(v);
                    893:       for (i = 0; i < v->len; i++ ) mb[i][col] = (Obj)ent[i];
                    894:     } else if ( OID(BDY(tn)) == O_LIST ) {
                    895:       for (row = 0, sn = BDY((LIST)BDY(tn)); sn; row++, sn = NEXT(sn) )
                    896:         mb[row][col] = (pointer)BDY(sn);
                    897:     } else {
                    898:       error("matc : invalid argument");
                    899:     }
                    900:   }
                    901:   *rp = m;
                    902: }
                    903:
                    904: void Pvtol(NODE arg,LIST *rp)
                    905: {
                    906:   NODE n,n1;
                    907:   VECT v;
                    908:   pointer *a;
                    909:   int len,i;
                    910:
                    911:   if ( OID(ARG0(arg)) == O_LIST ) {
                    912:     *rp = ARG0(arg);
                    913:     return;
                    914:   }
                    915:   asir_assert(ARG0(arg),O_VECT,"vtol");
                    916:   v = (VECT)ARG0(arg); len = v->len; a = BDY(v);
                    917:   for ( i = len - 1, n = 0; i >= 0; i-- ) {
                    918:     MKNODE(n1,a[i],n); n = n1;
                    919:   }
                    920:   MKLIST(*rp,n);
                    921: }
                    922:
                    923: void Pltov(NODE arg,VECT *rp)
                    924: {
                    925:   NODE n;
                    926:   VECT v,v0;
                    927:   int len,i;
                    928:
                    929:   if ( OID(ARG0(arg)) == O_VECT ) {
                    930:     v0 = (VECT)ARG0(arg); len = v0->len;
                    931:     MKVECT(v,len);
                    932:     for ( i = 0; i < len; i++ ) {
                    933:       BDY(v)[i] = BDY(v0)[i];
                    934:     }
                    935:     *rp = v;
                    936:     return;
                    937:   }
                    938:   asir_assert(ARG0(arg),O_LIST,"ltov");
                    939:   n = (NODE)BDY((LIST)ARG0(arg));
                    940:   len = length(n);
                    941:   MKVECT(v,len);
                    942:   for ( i = 0; i < len; i++, n = NEXT(n) )
                    943:     BDY(v)[i] = BDY(n);
                    944:   *rp = v;
                    945: }
                    946:
                    947: /* XXX it seems that remainder is not used */
                    948:
                    949: void Premainder(NODE arg,Obj *rp)
                    950: {
                    951:   Obj a;
                    952:   VECT v,w;
                    953:   MAT m,l;
                    954:   pointer *vb,*wb;
                    955:   pointer **mb,**lb;
                    956:   int id,i,j,n,row,col,t,smd,sgn;
                    957:   Z md,q,r;
                    958:
                    959:   a = (Obj)ARG0(arg); md = (Z)ARG1(arg);
                    960:   if ( !a )
                    961:     *rp = 0;
                    962:   else {
                    963:     id = OID(a);
                    964:     switch ( id ) {
                    965:       case O_N:
                    966:       case O_P:
                    967:         cmp(md,(P)a,(P *)rp); break;
                    968:       case O_VECT:
                    969:         /* strage spec */
1.2       noro      970:         smd = ZTOS(md);
1.1       noro      971:         v = (VECT)a; n = v->len; vb = v->body;
                    972:         MKVECT(w,n); wb = w->body;
                    973:         for ( i = 0; i < n; i++ ) {
                    974:           if ( sgnz(vb[i]) > 0 )
                    975:             remz(vb[i],md,(Z *)&wb[i]);
                    976:           else {
                    977:             absz(vb[i],&q);
                    978:             remz(q,md,&r);
                    979:             chsgnz(r,(Z *)&wb[i]);
                    980:           }
                    981:         }
                    982:         *rp = (Obj)w;
                    983:         break;
                    984:       case O_MAT:
                    985:         m = (MAT)a; row = m->row; col = m->col; mb = m->body;
                    986:         MKMAT(l,row,col); lb = l->body;
                    987:         for ( i = 0; i < row; i++ )
                    988:           for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
                    989:             cmp(md,(P)vb[j],(P *)&wb[j]);
                    990:         *rp = (Obj)l;
                    991:         break;
                    992:       default:
                    993:         error("remainder : invalid argument");
                    994:     }
                    995:   }
                    996: }
                    997:
                    998: void Psremainder(NODE arg,Obj *rp)
                    999: {
                   1000:   Obj a;
                   1001:   VECT v,w;
                   1002:   MAT m,l;
                   1003:   pointer *vb,*wb;
                   1004:   pointer **mb,**lb;
                   1005:   int id,i,j,n,row,col,t,smd,sgn;
                   1006:   Z md,q,r;
                   1007:
                   1008:   a = (Obj)ARG0(arg); md = (Z)ARG1(arg);
                   1009:   if ( !a )
                   1010:     *rp = 0;
                   1011:   else {
                   1012:     id = OID(a);
                   1013:     switch ( id ) {
                   1014:       case O_N:
                   1015:       case O_P:
                   1016:         cmp(md,(P)a,(P *)rp); break;
                   1017:       case O_VECT:
1.2       noro     1018:         smd = ZTOS(md);
1.1       noro     1019:         v = (VECT)a; n = v->len; vb = v->body;
                   1020:         MKVECT(w,n); wb = w->body;
                   1021:         for ( i = 0; i < n; i++ )
                   1022:           remz(vb[i],md,(Z *)&wb[i]);
                   1023:         *rp = (Obj)w;
                   1024:         break;
                   1025:       case O_MAT:
                   1026:         m = (MAT)a; row = m->row; col = m->col; mb = m->body;
                   1027:         MKMAT(l,row,col); lb = l->body;
                   1028:         for ( i = 0; i < row; i++ )
                   1029:           for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
                   1030:             cmp(md,(P)vb[j],(P *)&wb[j]);
                   1031:         *rp = (Obj)l;
                   1032:         break;
                   1033:       default:
                   1034:         error("remainder : invalid argument");
                   1035:     }
                   1036:   }
                   1037: }
                   1038:
                   1039: void Psize(NODE arg,LIST *rp)
                   1040: {
                   1041:
                   1042:   int n,m;
                   1043:   Z q;
                   1044:   NODE t,s;
                   1045:
                   1046:   if ( !ARG0(arg) )
                   1047:      t = 0;
                   1048:   else {
                   1049:     switch (OID(ARG0(arg))) {
                   1050:       case O_VECT:
                   1051:         n = ((VECT)ARG0(arg))->len;
1.2       noro     1052:         STOZ(n,q); MKNODE(t,q,0);
1.1       noro     1053:         break;
                   1054:       case O_MAT:
                   1055:         n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;
1.2       noro     1056:         STOZ(m,q); MKNODE(s,q,0); STOZ(n,q); MKNODE(t,q,s);
1.1       noro     1057:         break;
                   1058:       case O_IMAT:
                   1059:         n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col;
1.2       noro     1060:         STOZ(m,q); MKNODE(s,q,0); STOZ(n,q); MKNODE(t,q,s);
1.1       noro     1061:         break;
                   1062:       default:
                   1063:         error("size : invalid argument"); break;
                   1064:     }
                   1065:   }
                   1066:   MKLIST(*rp,t);
                   1067: }
                   1068:
                   1069: void Pdet(NODE arg,P *rp)
                   1070: {
                   1071:   MAT m;
                   1072:   int n,i,j,mod;
                   1073:   P d;
                   1074:   P **mat,**w;
                   1075:
                   1076:   m = (MAT)ARG0(arg);
                   1077:   asir_assert(m,O_MAT,"det");
                   1078:   if ( m->row != m->col )
                   1079:     error("det : non-square matrix");
                   1080:   else if ( argc(arg) == 1 )
                   1081:     detp(CO,(P **)BDY(m),m->row,rp);
                   1082:   else {
1.2       noro     1083:     n = m->row; mod = ZTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
1.1       noro     1084:     w = (P **)almat_pointer(n,n);
                   1085:     for ( i = 0; i < n; i++ )
                   1086:       for ( j = 0; j < n; j++ )
                   1087:         ptomp(mod,mat[i][j],&w[i][j]);
                   1088:     detmp(CO,mod,w,n,&d);
                   1089:     mptop(d,rp);
                   1090:   }
                   1091: }
                   1092:
                   1093: void Pinvmat(NODE arg,LIST *rp)
                   1094: {
                   1095:   MAT m,r;
                   1096:   int n,i,j,mod;
                   1097:   P dn;
                   1098:   P **mat,**imat,**w;
                   1099:   NODE nd;
                   1100:
                   1101:   m = (MAT)ARG0(arg);
                   1102:   asir_assert(m,O_MAT,"invmat");
                   1103:   if ( m->row != m->col )
                   1104:     error("invmat : non-square matrix");
                   1105:   else if ( argc(arg) == 1 ) {
                   1106:     n = m->row;
                   1107:     invmatp(CO,(P **)BDY(m),n,&imat,&dn);
                   1108:     NEWMAT(r); r->row = n; r->col = n; r->body = (pointer **)imat;
                   1109:     nd = mknode(2,r,dn);
                   1110:     MKLIST(*rp,nd);
                   1111:   } else {
1.2       noro     1112:     n = m->row; mod = ZTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
1.1       noro     1113:     w = (P **)almat_pointer(n,n);
                   1114:     for ( i = 0; i < n; i++ )
                   1115:       for ( j = 0; j < n; j++ )
                   1116:         ptomp(mod,mat[i][j],&w[i][j]);
                   1117: #if 0
                   1118:     detmp(CO,mod,w,n,&d);
                   1119:     mptop(d,rp);
                   1120: #else
                   1121:     error("not implemented yet");
                   1122: #endif
                   1123:   }
                   1124: }
                   1125:
                   1126: /*
                   1127:   input : a row x col matrix A
                   1128:     A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
                   1129:
                   1130:   output : [B,D,R,C]
                   1131:     B : a rank(A) x col-rank(A) matrix
                   1132:     D : the denominator
                   1133:     R : a vector of length rank(A)
                   1134:     C : a vector of length col-rank(A)
                   1135:     B[I] <-> D*x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
                   1136: */
                   1137:
                   1138: void Pgeneric_gauss_elim(NODE arg,LIST *rp)
                   1139: {
                   1140:   NODE n0,opt,p;
                   1141:   MAT m,nm;
                   1142:   int *ri,*ci;
                   1143:   VECT rind,cind;
                   1144:   Z dn,q;
                   1145:   int i,row,col,t,rank;
                   1146:   int is_hensel = 0;
                   1147:   char *key;
                   1148:   Obj value;
                   1149:
                   1150:   if ( current_option ) {
                   1151:     for ( opt = current_option; opt; opt = NEXT(opt) ) {
                   1152:       p = BDY((LIST)BDY(opt));
                   1153:       key = BDY((STRING)BDY(p));
                   1154:       value = (Obj)BDY(NEXT(p));
                   1155:       if ( !strcmp(key,"hensel") && value ) {
                   1156:         is_hensel = value ? 1 : 0;
                   1157:         break;
                   1158:       }
                   1159:     }
                   1160:   }
                   1161:   asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim");
                   1162:   m = (MAT)ARG0(arg);
                   1163:   row = m->row; col = m->col;
                   1164:   if ( is_hensel )
                   1165:     rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci);
                   1166:   else
1.3     ! noro     1167:     rank = generic_gauss_elim64(m,&nm,&dn,&ri,&ci);
1.1       noro     1168:   t = col-rank;
                   1169:   MKVECT(rind,rank);
                   1170:   MKVECT(cind,t);
                   1171:   for ( i = 0; i < rank; i++ ) {
1.2       noro     1172:     STOZ(ri[i],q);
1.1       noro     1173:     BDY(rind)[i] = (pointer)q;
                   1174:   }
                   1175:   for ( i = 0; i < t; i++ ) {
1.2       noro     1176:     STOZ(ci[i],q);
1.1       noro     1177:     BDY(cind)[i] = (pointer)q;
                   1178:   }
                   1179:   n0 = mknode(4,nm,dn,rind,cind);
                   1180:   MKLIST(*rp,n0);
                   1181: }
                   1182:
                   1183: int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat);
                   1184:
                   1185: void Pindep_rows_mod(NODE arg,VECT *rp)
                   1186: {
                   1187:   MAT m,mat;
                   1188:   VECT rind;
                   1189:   Z **tmat;
                   1190:   int **wmat,**row0;
                   1191:   Z *rib;
                   1192:   int *rowstat,*p;
                   1193:   Z q;
                   1194:   int md,i,j,k,l,row,col,t,rank;
                   1195:
                   1196:   asir_assert(ARG0(arg),O_MAT,"indep_rows_mod");
                   1197:   asir_assert(ARG1(arg),O_N,"indep_rows_mod");
1.2       noro     1198:   m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1       noro     1199:   row = m->row; col = m->col; tmat = (Z **)m->body;
                   1200:   wmat = (int **)almat(row,col);
                   1201:
                   1202:   row0 = (int **)ALLOCA(row*sizeof(int *));
                   1203:   for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
                   1204:
                   1205:   rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
                   1206:   for ( i = 0; i < row; i++ )
                   1207:     for ( j = 0; j < col; j++ )
                   1208:       wmat[i][j] = remqi((Q)tmat[i][j],md);
                   1209:   rank = indep_rows_mod(wmat,row,col,md,rowstat);
                   1210:
                   1211:   MKVECT(rind,rank);
                   1212:   rib = (Z *)rind->body;
                   1213:   for ( j = 0; j < rank; j++ ) {
1.2       noro     1214:     STOZ(rowstat[j],rib[j]);
1.1       noro     1215:   }
                   1216:     *rp = rind;
                   1217: }
                   1218:
                   1219: /*
                   1220:   input : a row x col matrix A
                   1221:     A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
                   1222:
                   1223:   output : [B,R,C]
                   1224:     B : a rank(A) x col-rank(A) matrix
                   1225:     R : a vector of length rank(A)
                   1226:     C : a vector of length col-rank(A)
                   1227:     RN : a vector of length rank(A) indicating useful rows
                   1228:
                   1229:     B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
                   1230: */
                   1231:
1.3     ! noro     1232: #if SIZEOF_LONG==8
        !          1233: void Pgeneric_gauss_elim_mod64(NODE arg,LIST *rp)
        !          1234: {
        !          1235:   NODE n0;
        !          1236:   MAT m,mat;
        !          1237:   VECT rind,cind,rnum;
        !          1238:   Z **tmat;
        !          1239:   mp_limb_t **wmat,**row0;
        !          1240:   Z *rib,*cib,*rnb;
        !          1241:   int *colstat;
        !          1242:   mp_limb_t *p;
        !          1243:   Z q;
        !          1244:   mp_limb_t md;
        !          1245:   int i,j,k,l,row,col,t,rank;
        !          1246:
        !          1247:   asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod64");
        !          1248:   asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod64");
        !          1249:   m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
        !          1250:   row = m->row; col = m->col; tmat = (Z **)m->body;
        !          1251:   wmat = (mp_limb_t **)almat64(row,col);
        !          1252:
        !          1253:   row0 = (mp_limb_t **)ALLOCA(row*sizeof(mp_limb_t *));
        !          1254:   for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
        !          1255:
        !          1256:   colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
        !          1257:   for ( i = 0; i < row; i++ )
        !          1258:     for ( j = 0; j < col; j++ )
        !          1259:       wmat[i][j] = remqi64((Q)tmat[i][j],md);
        !          1260:   rank = generic_gauss_elim_mod64(wmat,row,col,md,colstat);
        !          1261:
        !          1262:   MKVECT(rnum,rank);
        !          1263:   rnb = (Z *)rnum->body;
        !          1264:   for ( i = 0; i < rank; i++ )
        !          1265:     for ( j = 0, p = wmat[i]; j < row; j++ )
        !          1266:       if ( p == row0[j] )
        !          1267:         STOZ(j,rnb[i]);
        !          1268:
        !          1269:   MKMAT(mat,rank,col-rank);
        !          1270:   tmat = (Z **)mat->body;
        !          1271:   for ( i = 0; i < rank; i++ )
        !          1272:     for ( j = k = 0; j < col; j++ )
        !          1273:       if ( !colstat[j] ) {
        !          1274:         UTOZ(wmat[i][j],tmat[i][k]); k++;
        !          1275:       }
        !          1276:
        !          1277:   MKVECT(rind,rank);
        !          1278:   MKVECT(cind,col-rank);
        !          1279:   rib = (Z *)rind->body; cib = (Z *)cind->body;
        !          1280:   for ( j = k = l = 0; j < col; j++ )
        !          1281:     if ( colstat[j] ) {
        !          1282:       STOZ(j,rib[k]); k++;
        !          1283:     } else {
        !          1284:       STOZ(j,cib[l]); l++;
        !          1285:     }
        !          1286:   n0 = mknode(4,mat,rind,cind,rnum);
        !          1287:   MKLIST(*rp,n0);
        !          1288: }
        !          1289: #endif
        !          1290:
1.1       noro     1291: void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
                   1292: {
                   1293:   NODE n0;
                   1294:   MAT m,mat;
                   1295:   VECT rind,cind,rnum;
                   1296:   Z **tmat;
                   1297:   int **wmat,**row0;
                   1298:   Z *rib,*cib,*rnb;
                   1299:   int *colstat,*p;
                   1300:   Z q;
1.3     ! noro     1301:   long mdl;
1.1       noro     1302:   int md,i,j,k,l,row,col,t,rank;
                   1303:
                   1304:   asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
                   1305:   asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
1.3     ! noro     1306: #if SIZEOF_LONG==8
        !          1307:   mdl = ZTOS((Z)ARG1(arg));
        !          1308:   if ( mdl >= ((mp_limb_t)1)<<32 ) {
        !          1309:     Pgeneric_gauss_elim_mod64(arg,rp);
        !          1310:     return;
        !          1311:   }
        !          1312: #endif
        !          1313:   m = (MAT)ARG0(arg);
        !          1314:   md = ZTOS((Z)ARG1(arg));
1.1       noro     1315:   row = m->row; col = m->col; tmat = (Z **)m->body;
                   1316:   wmat = (int **)almat(row,col);
                   1317:
                   1318:   row0 = (int **)ALLOCA(row*sizeof(int *));
                   1319:   for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
                   1320:
                   1321:   colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
                   1322:   for ( i = 0; i < row; i++ )
                   1323:     for ( j = 0; j < col; j++ )
                   1324:       wmat[i][j] = remqi((Q)tmat[i][j],md);
                   1325:   rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
                   1326:
                   1327:   MKVECT(rnum,rank);
                   1328:   rnb = (Z *)rnum->body;
                   1329:   for ( i = 0; i < rank; i++ )
                   1330:     for ( j = 0, p = wmat[i]; j < row; j++ )
                   1331:       if ( p == row0[j] )
1.2       noro     1332:         STOZ(j,rnb[i]);
1.1       noro     1333:
                   1334:   MKMAT(mat,rank,col-rank);
                   1335:   tmat = (Z **)mat->body;
                   1336:   for ( i = 0; i < rank; i++ )
                   1337:     for ( j = k = 0; j < col; j++ )
                   1338:       if ( !colstat[j] ) {
1.2       noro     1339:         UTOZ(wmat[i][j],tmat[i][k]); k++;
1.1       noro     1340:       }
                   1341:
                   1342:   MKVECT(rind,rank);
                   1343:   MKVECT(cind,col-rank);
                   1344:   rib = (Z *)rind->body; cib = (Z *)cind->body;
                   1345:   for ( j = k = l = 0; j < col; j++ )
                   1346:     if ( colstat[j] ) {
1.2       noro     1347:       STOZ(j,rib[k]); k++;
1.1       noro     1348:     } else {
1.2       noro     1349:       STOZ(j,cib[l]); l++;
1.1       noro     1350:     }
                   1351:   n0 = mknode(4,mat,rind,cind,rnum);
                   1352:   MKLIST(*rp,n0);
                   1353: }
                   1354:
1.3     ! noro     1355:
1.1       noro     1356: void Pleqm(NODE arg,VECT *rp)
                   1357: {
                   1358:   MAT m;
                   1359:   VECT vect;
                   1360:   pointer **mat;
                   1361:   Z *v;
                   1362:   Z q;
                   1363:   int **wmat;
                   1364:   int md,i,j,row,col,t,n,status;
                   1365:
                   1366:   asir_assert(ARG0(arg),O_MAT,"leqm");
                   1367:   asir_assert(ARG1(arg),O_N,"leqm");
1.2       noro     1368:   m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1       noro     1369:   row = m->row; col = m->col; mat = m->body;
                   1370:   wmat = (int **)almat(row,col);
                   1371:   for ( i = 0; i < row; i++ )
                   1372:     for ( j = 0; j < col; j++ )
                   1373:       wmat[i][j] = remqi((Q)mat[i][j],md);
                   1374:   status = gauss_elim_mod(wmat,row,col,md);
                   1375:   if ( status < 0 )
                   1376:     *rp = 0;
                   1377:   else if ( status > 0 )
                   1378:     *rp = (VECT)ONE;
                   1379:   else {
                   1380:     n = col - 1;
                   1381:     MKVECT(vect,n);
                   1382:     for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {
1.2       noro     1383:       t = (md-wmat[i][n])%md; STOZ(t,v[i]);
1.1       noro     1384:     }
                   1385:     *rp = vect;
                   1386:   }
                   1387: }
                   1388:
                   1389: int gauss_elim_mod(int **mat,int row,int col,int md)
                   1390: {
                   1391:   int i,j,k,inv,a,n;
                   1392:   int *t,*pivot;
                   1393:
                   1394:   n = col - 1;
                   1395:   for ( j = 0; j < n; j++ ) {
                   1396:     for ( i = j; i < row && !mat[i][j]; i++ );
                   1397:     if ( i == row )
                   1398:       return 1;
                   1399:     if ( i != j ) {
                   1400:       t = mat[i]; mat[i] = mat[j]; mat[j] = t;
                   1401:     }
                   1402:     pivot = mat[j];
                   1403:     inv = invm(pivot[j],md);
                   1404:     for ( k = j; k <= n; k++ ) {
                   1405: /*      pivot[k] = dmar(pivot[k],inv,0,md); */
                   1406:       DMAR(pivot[k],inv,0,md,pivot[k])
                   1407:     }
                   1408:     for ( i = 0; i < row; i++ ) {
                   1409:       t = mat[i];
                   1410:       if ( i != j && (a = t[j]) )
                   1411:         for ( k = j, a = md - a; k <= n; k++ ) {
                   1412:           unsigned int tk;
                   1413: /*          t[k] = dmar(pivot[k],a,t[k],md); */
                   1414:           DMAR(pivot[k],a,t[k],md,tk)
                   1415:           t[k] = tk;
                   1416:         }
                   1417:     }
                   1418:   }
                   1419:   for ( i = n; i < row && !mat[i][n]; i++ );
                   1420:   if ( i == row )
                   1421:     return 0;
                   1422:   else
                   1423:     return -1;
                   1424: }
                   1425:
                   1426: struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;
                   1427: struct oEGT eg_conv;
                   1428:
                   1429: #if 0
                   1430: void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm);
                   1431:
                   1432: /* XXX broken */
                   1433: void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
                   1434: {
                   1435:   Q **a0,**b;
                   1436:   Q *aiq;
                   1437:   N **a;
                   1438:   N *ai;
                   1439:   Q q,q1,dn2,a1,q0,bik;
                   1440:   MAT m;
                   1441:   unsigned int md;
                   1442:   int n,ind,i,j,rank,t,inv,t1,ret,min,k;
                   1443:   int **w;
                   1444:   int *wi,*rinfo0,*rinfo;
                   1445:   N m1,m2,m3,u,s;
                   1446:
                   1447:   a0 = (Q **)mat->body;
                   1448:   n = mat->row;
                   1449:   if ( n != mat->col )
                   1450:     error("lu_dec_cr : non-square matrix");
                   1451:   w = (int **)almat(n,n);
                   1452:   MKMAT(m,n,n);
                   1453:   a = (N **)m->body;
                   1454:   UTON(1,m1);
                   1455:   rinfo0 = 0;
                   1456:   ind = 0;
                   1457:   while ( 1 ) {
                   1458:     md = get_lprime(ind);
                   1459:     /* mat mod md */
                   1460:     for ( i = 0; i < n; i++ )
                   1461:       for ( j = 0, aiq = a0[i], wi = w[i]; j < n; j++ )
                   1462:         if ( q = aiq[j] ) {
                   1463:           t = rem(NM(q),md);
                   1464:           if ( t && SGN(q) < 0 )
                   1465:             t = (md - t) % md;
                   1466:           wi[j] = t;
                   1467:         } else
                   1468:           wi[j] = 0;
                   1469:
                   1470:     if ( !lu_mod((unsigned int **)w,n,md,&rinfo) ) continue;
                   1471:     printf("."); fflush(stdout);
                   1472:     if ( !rinfo0 )
                   1473:       *perm = rinfo0 = rinfo;
                   1474:     else {
                   1475:       for ( i = 0; i < n; i++ )
                   1476:         if ( rinfo[i] != rinfo0[i] ) break;
                   1477:       if ( i < n ) continue;
                   1478:     }
                   1479:     if ( UNIN(m1) ) {
                   1480:       for ( i = 0; i < n; i++ )
                   1481:         for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) {
                   1482:           UTON(wi[j],u); ai[j] = u;
                   1483:         }
                   1484:       UTON(md,m1);
                   1485:     } else {
                   1486:       inv = invm(rem(m1,md),md);
                   1487:       UTON(md,m2); muln(m1,m2,&m3);
                   1488:       for ( i = 0; i < n; i++ )
                   1489:         for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ )
                   1490:           if ( ai[i] ) {
                   1491:           /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                   1492:             t = rem(ai[j],md);
                   1493:             if ( wi[j] >= t )
                   1494:               t = wi[j]-t;
                   1495:             else
                   1496:               t = md-(t-wi[j]);
                   1497:             DMAR(t,inv,0,md,t1)
                   1498:             UTON(t1,u);
                   1499:             muln(m1,u,&s);
                   1500:             addn(ai[j],s,&u); ai[j] = u;
                   1501:           } else if ( wi[j] ) {
                   1502:             /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                   1503:             DMAR(wi[j],inv,0,md,t)
                   1504:             UTON(t,u);
                   1505:             muln(m1,u,&s); ai[j] = s;
                   1506:           }
                   1507:       m1 = m3;
                   1508:     }
                   1509:     if ( (++ind%8) == 0 ) {
                   1510:       ret = intmtoratm(m,m1,lu,dn);
                   1511:       if ( ret ) {
                   1512:         b = (Q **)lu->body;
                   1513:         mulq(*dn,*dn,&dn2);
                   1514:         for ( i = 0; i < n; i++ ) {
                   1515:           for ( j = 0; j < n; j++ ) {
                   1516:             q = 0;
                   1517:             min = MIN(i,j);
                   1518:             for ( k = 0; k <= min; k++ ) {
                   1519:               bik = k==i ? *dn : b[i][k];
                   1520:               mulq(bik,b[k][j],&q0);
                   1521:               addq(q,q0,&q1); q = q1;
                   1522:             }
                   1523:             mulq(a0[rinfo0[i]][j],dn2,&q1);
                   1524:             if ( cmpq(q,q1) ) break;
                   1525:           }
                   1526:           if ( j < n ) break;
                   1527:         }
                   1528:         if ( i == n )
                   1529:           return;
                   1530:       }
                   1531:     }
                   1532:   }
                   1533: }
                   1534: #endif
                   1535:
                   1536:
                   1537: int f4_nocheck;
                   1538:
                   1539: #define ONE_STEP1  if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
                   1540:
                   1541: void reduce_reducers_mod(int **mat,int row,int col,int md)
                   1542: {
                   1543:   int i,j,k,l,hc,zzz;
                   1544:   int *t,*s,*tj,*ind;
                   1545:
                   1546:   /* reduce the reducers */
                   1547:   ind = (int *)ALLOCA(row*sizeof(int));
                   1548:   for ( i = 0; i < row; i++ ) {
                   1549:     t = mat[i];
                   1550:     for ( j = 0; j < col && !t[j]; j++ );
                   1551:     /* register the position of the head term */
                   1552:     ind[i] = j;
                   1553:     for ( l = i-1; l >= 0; l-- ) {
                   1554:       /* reduce mat[i] by mat[l] */
                   1555:       if ( hc = t[ind[l]] ) {
                   1556:         /* mat[i] = mat[i]-hc*mat[l] */
                   1557:         j = ind[l];
                   1558:         s = mat[l]+j;
                   1559:         tj = t+j;
                   1560:         hc = md-hc;
                   1561:         k = col-j;
                   1562:         for ( ; k >= 64; k -= 64 ) {
                   1563:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1564:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1565:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1566:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1567:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1568:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1569:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1570:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1571:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1572:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1573:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1574:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1575:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1576:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1577:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1578:           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                   1579:         }
                   1580:         for ( ; k > 0; k-- ) {
                   1581:           if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
                   1582:         }
                   1583:       }
                   1584:     }
                   1585:   }
                   1586: }
                   1587:
                   1588: /*
                   1589:   mat[i] : reducers (i=0,...,nred-1)
                   1590:            spolys (i=nred,...,row-1)
                   1591:   mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
                   1592:   1. reduce the reducers
                   1593:   2. reduce spolys by the reduced reducers
                   1594: */
                   1595:
                   1596: void pre_reduce_mod(int **mat,int row,int col,int nred,int md)
                   1597: {
                   1598:   int i,j,k,l,hc,inv;
                   1599:   int *t,*s,*tk,*ind;
                   1600:
                   1601: #if 1
                   1602:   /* reduce the reducers */
                   1603:   ind = (int *)ALLOCA(row*sizeof(int));
                   1604:   for ( i = 0; i < nred; i++ ) {
                   1605:     /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */
                   1606:     t = mat[i];
                   1607:     for ( j = 0; j < col && !t[j]; j++ );
                   1608:     /* register the position of the head term */
                   1609:     ind[i] = j;
                   1610:     inv = invm(t[j],md);
                   1611:     for ( k = j; k < col; k++ )
                   1612:       if ( t[k] )
                   1613:         DMAR(t[k],inv,0,md,t[k])
                   1614:     for ( l = i-1; l >= 0; l-- ) {
                   1615:       /* reduce mat[i] by mat[l] */
                   1616:       if ( hc = t[ind[l]] ) {
                   1617:         /* mat[i] = mat[i]-hc*mat[l] */
                   1618:         for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
                   1619:           k < col; k++, tk++, s++ )
                   1620:           if ( *s )
                   1621:             DMAR(*s,hc,*tk,md,*tk)
                   1622:       }
                   1623:     }
                   1624:   }
                   1625:   /* reduce the spolys */
                   1626:   for ( i = nred; i < row; i++ ) {
                   1627:     t = mat[i];
                   1628:     for ( l = nred-1; l >= 0; l-- ) {
                   1629:       /* reduce mat[i] by mat[l] */
                   1630:       if ( hc = t[ind[l]] ) {
                   1631:         /* mat[i] = mat[i]-hc*mat[l] */
                   1632:         for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
                   1633:           k < col; k++, tk++, s++ )
                   1634:           if ( *s )
                   1635:             DMAR(*s,hc,*tk,md,*tk)
                   1636:       }
                   1637:     }
                   1638:   }
                   1639: #endif
                   1640: }
                   1641: /*
                   1642:   mat[i] : reducers (i=0,...,nred-1)
                   1643:   mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
                   1644: */
                   1645:
                   1646: void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md)
                   1647: {
                   1648:   int i,j,k,hc,zzz;
                   1649:   int *s,*tj;
                   1650:
                   1651:   /* reduce the spolys by redmat */
                   1652:   for ( i = nred-1; i >= 0; i-- ) {
                   1653:     /* reduce sp by redmat[i] */
                   1654:     if ( hc = sp[ind[i]] ) {
                   1655:       /* sp = sp-hc*redmat[i] */
                   1656:       j = ind[i];
                   1657:       hc = md-hc;
                   1658:       s = redmat[i]+j;
                   1659:       tj = sp+j;
                   1660:       for ( k = col-j; k > 0; k-- ) {
                   1661:         if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
                   1662:       }
                   1663:     }
                   1664:   }
                   1665: }
                   1666:
                   1667: /*
                   1668:   mat[i] : compressed reducers (i=0,...,nred-1)
                   1669:   mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
                   1670: */
                   1671:
                   1672: void red_by_compress(int m,unsigned int *p,unsigned int *r,
                   1673:   unsigned int *ri,unsigned int hc,int len)
                   1674: {
                   1675:   unsigned int up,lo;
                   1676:   unsigned int dmy;
                   1677:   unsigned int *pj;
                   1678:
                   1679:   p[*ri] = 0; r++; ri++;
                   1680:   for ( len--; len; len--, r++, ri++ ) {
                   1681:     pj = p+ *ri;
                   1682:     DMA(*r,hc,*pj,up,lo);
                   1683:     if ( up ) {
                   1684:       DSAB(m,up,lo,dmy,*pj);
                   1685:     } else
                   1686:       *pj = lo;
                   1687:   }
                   1688: }
                   1689:
                   1690: /* p -= hc*r */
                   1691:
                   1692: void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
                   1693: {
                   1694:   unsigned int up,lo,dmy;
                   1695:
                   1696:   *p++ = 0; r++; len--;
                   1697:   for ( ; len; len--, r++, p++ )
                   1698:     if ( *r ) {
                   1699:       DMA(*r,hc,*p,up,lo);
                   1700:       if ( up ) {
                   1701:         DSAB(m,up,lo,dmy,*p);
                   1702:       } else
                   1703:         *p = lo;
                   1704:     }
                   1705: }
                   1706:
1.3     ! noro     1707: #if SIZEOF_LONG==8
        !          1708: /* mp_limb_t vector += U32 vector(normalized)*U32 */
1.1       noro     1709:
1.3     ! noro     1710: void red_by_vect64(int m, mp_limb_t *p,unsigned int *c,mp_limb_t *r,unsigned int hc,int len)
1.1       noro     1711: {
1.3     ! noro     1712:   mp_limb_t t;
1.1       noro     1713:
                   1714:   /* (p[0],c[0]) is normalized */
                   1715:   *p++ = 0; *c++ = 0; r++; len--;
                   1716:   for ( ; len; len--, r++, p++, c++ )
                   1717:     if ( *r ) {
                   1718:       t = (*p)+(*r)*hc;
                   1719:       if ( t < *p ) (*c)++;
                   1720:       *p = t;
                   1721:     }
                   1722: }
1.3     ! noro     1723:
        !          1724: /* mp_limb_t vector = (mp_limb_t vector+mp_limb_t vector*mp_limb_t)%mp_limb_t */
        !          1725:
        !          1726: void red_by_vect64mod(mp_limb_t m, mp_limb_t *p,mp_limb_t *r,mp_limb_t hc,int len)
        !          1727: {
        !          1728:   *p++ = 0; r++; len--;
        !          1729:   for ( ; len; len--, r++, p++ )
        !          1730:     if ( *r )
        !          1731:       *p = muladdmod64(*r,hc,*p,m);
        !          1732: }
        !          1733:
        !          1734: int generic_gauss_elim_mod64(mp_limb_t **mat,int row,int col,mp_limb_t md,int *colstat)
        !          1735: {
        !          1736:   int i,j,k,l,rank;
        !          1737:   mp_limb_t inv,a;
        !          1738:   mp_limb_t *t,*pivot,*pk;
        !          1739:
        !          1740:   for ( rank = 0, j = 0; j < col; j++ ) {
        !          1741:     for ( i = rank; i < row; i++ )
        !          1742:       if ( mat[i][j] )
        !          1743:         break;
        !          1744:     if ( i == row ) {
        !          1745:       colstat[j] = 0;
        !          1746:       continue;
        !          1747:     } else
        !          1748:       colstat[j] = 1;
        !          1749:     if ( i != rank ) {
        !          1750:       t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
        !          1751:     }
        !          1752:     pivot = mat[rank];
        !          1753:     inv = invmod64(pivot[j],md);
        !          1754:     for ( k = j, pk = pivot+k; k < col; k++, pk++ )
        !          1755:       if ( *pk )
        !          1756:         *pk = mulmod64(*pk,inv,md);
        !          1757:     for ( i = rank+1; i < row; i++ ) {
        !          1758:       t = mat[i];
        !          1759:       if ( a = t[j] )
        !          1760:         red_by_vect64mod(md,t+j,pivot+j,md-a,col-j);
        !          1761:     }
        !          1762:     rank++;
        !          1763:   }
        !          1764:   for ( j = col-1, l = rank-1; j >= 0; j-- )
        !          1765:     if ( colstat[j] ) {
        !          1766:       pivot = mat[l];
        !          1767:       for ( i = 0; i < l; i++ ) {
        !          1768:         t = mat[i];
        !          1769:         if ( a = t[j] )
        !          1770:           red_by_vect64mod(md,t+j,pivot+j,md-a,col-j);
        !          1771:       }
        !          1772:       l--;
        !          1773:     }
        !          1774:   return rank;
        !          1775: }
        !          1776:
1.1       noro     1777: #endif
                   1778:
                   1779: void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
                   1780: {
                   1781:   *p++ = 0; r++; len--;
                   1782:   for ( ; len; len--, r++, p++ )
                   1783:     if ( *r )
                   1784:       *p = _addsf(_mulsf(*r,hc),*p);
                   1785: }
                   1786:
                   1787: extern Z current_mod_lf;
                   1788: extern int current_mod_lf_size;
                   1789:
                   1790: void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len)
                   1791: {
                   1792:   mpz_set_ui(*p++,0); r++; len--;
                   1793:   for ( ; len; len--, r++, p++ ) {
                   1794:        mpz_addmul(*p,*r,hc);
                   1795: #if 0
                   1796:        if ( mpz_size(*p) > current_mod_lf_size )
                   1797:          mpz_mod(*p,*p,BDY(current_mod_lf));
                   1798: #endif
                   1799:     }
                   1800: }
                   1801:
                   1802:
                   1803: extern unsigned int **psca;
                   1804:
                   1805: void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,
                   1806:   int nred,int col,int md)
                   1807: {
                   1808:   int i,len;
                   1809:   CDP ri;
                   1810:   unsigned int hc;
                   1811:   unsigned int *usp;
                   1812:
                   1813:   usp = (unsigned int *)sp;
                   1814:   /* reduce the spolys by redmat */
                   1815:   for ( i = nred-1; i >= 0; i-- ) {
                   1816:     /* reduce sp by redmat[i] */
                   1817:     usp[ind[i]] %= md;
                   1818:     if ( hc = usp[ind[i]] ) {
                   1819:       /* sp = sp-hc*redmat[i] */
                   1820:       hc = md-hc;
                   1821:       ri = redmat[i];
                   1822:       len = ri->len;
                   1823:       red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len);
                   1824:     }
                   1825:   }
                   1826:   for ( i = 0; i < col; i++ )
                   1827:     if ( usp[i] >= (unsigned int)md )
                   1828:       usp[i] %= md;
                   1829: }
                   1830:
                   1831: #define ONE_STEP2  if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
                   1832:
                   1833: int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat)
                   1834: {
                   1835:   int i,j,k,l,inv,a,rank;
                   1836:   unsigned int *t,*pivot,*pk;
                   1837:   unsigned int **mat;
                   1838:
                   1839:   mat = (unsigned int **)mat0;
                   1840:   for ( rank = 0, j = 0; j < col; j++ ) {
                   1841:     for ( i = rank; i < row; i++ )
                   1842:       mat[i][j] %= md;
                   1843:     for ( i = rank; i < row; i++ )
                   1844:       if ( mat[i][j] )
                   1845:         break;
                   1846:     if ( i == row ) {
                   1847:       colstat[j] = 0;
                   1848:       continue;
                   1849:     } else
                   1850:       colstat[j] = 1;
                   1851:     if ( i != rank ) {
                   1852:       t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                   1853:     }
                   1854:     pivot = mat[rank];
                   1855:     inv = invm(pivot[j],md);
                   1856:     for ( k = j, pk = pivot+k; k < col; k++, pk++ )
                   1857:       if ( *pk ) {
                   1858:         if ( *pk >= (unsigned int)md )
                   1859:           *pk %= md;
                   1860:         DMAR(*pk,inv,0,md,*pk)
                   1861:       }
                   1862:     for ( i = rank+1; i < row; i++ ) {
                   1863:       t = mat[i];
                   1864:       if ( a = t[j] )
                   1865:         red_by_vect(md,t+j,pivot+j,md-a,col-j);
                   1866:     }
                   1867:     rank++;
                   1868:   }
                   1869:   for ( j = col-1, l = rank-1; j >= 0; j-- )
                   1870:     if ( colstat[j] ) {
                   1871:       pivot = mat[l];
                   1872:       for ( i = 0; i < l; i++ ) {
                   1873:         t = mat[i];
                   1874:         t[j] %= md;
                   1875:         if ( a = t[j] )
                   1876:           red_by_vect(md,t+j,pivot+j,md-a,col-j);
                   1877:       }
                   1878:       l--;
                   1879:     }
                   1880:   for ( j = 0, l = 0; l < rank; j++ )
                   1881:     if ( colstat[j] ) {
                   1882:       t = mat[l];
                   1883:       for ( k = j; k < col; k++ )
                   1884:         if ( t[k] >= (unsigned int)md )
                   1885:           t[k] %= md;
                   1886:       l++;
                   1887:     }
                   1888:   return rank;
                   1889: }
                   1890:
                   1891: int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,int *rowstat)
                   1892: {
                   1893:   int i,j,k,l,inv,a,rank;
                   1894:   unsigned int *t,*pivot,*pk;
                   1895:   unsigned int **mat;
                   1896:
                   1897:   for ( i = 0; i < row; i++ ) rowstat[i] = i;
                   1898:   mat = (unsigned int **)mat0;
                   1899:   for ( rank = 0, j = 0; j < col; j++ ) {
                   1900:     for ( i = rank; i < row; i++ )
                   1901:       mat[i][j] %= md;
                   1902:     for ( i = rank; i < row; i++ )
                   1903:       if ( mat[i][j] )
                   1904:         break;
                   1905:     if ( i == row ) {
                   1906:       colstat[j] = 0;
                   1907:       continue;
                   1908:     } else
                   1909:       colstat[j] = 1;
                   1910:     if ( i != rank ) {
                   1911:       t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                   1912:       k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
                   1913:     }
                   1914:     pivot = mat[rank];
                   1915:     inv = invm(pivot[j],md);
                   1916:     for ( k = j, pk = pivot+k; k < col; k++, pk++ )
                   1917:       if ( *pk ) {
                   1918:         if ( *pk >= (unsigned int)md )
                   1919:           *pk %= md;
                   1920:         DMAR(*pk,inv,0,md,*pk)
                   1921:       }
                   1922:     for ( i = rank+1; i < row; i++ ) {
                   1923:       t = mat[i];
                   1924:       if ( a = t[j] )
                   1925:         red_by_vect(md,t+j,pivot+j,md-a,col-j);
                   1926:     }
                   1927:     rank++;
                   1928:   }
                   1929:   for ( j = col-1, l = rank-1; j >= 0; j-- )
                   1930:     if ( colstat[j] ) {
                   1931:       pivot = mat[l];
                   1932:       for ( i = 0; i < l; i++ ) {
                   1933:         t = mat[i];
                   1934:         t[j] %= md;
                   1935:         if ( a = t[j] )
                   1936:           red_by_vect(md,t+j,pivot+j,md-a,col-j);
                   1937:       }
                   1938:       l--;
                   1939:     }
                   1940:   for ( j = 0, l = 0; l < rank; j++ )
                   1941:     if ( colstat[j] ) {
                   1942:       t = mat[l];
                   1943:       for ( k = j; k < col; k++ )
                   1944:         if ( t[k] >= (unsigned int)md )
                   1945:           t[k] %= md;
                   1946:       l++;
                   1947:     }
                   1948:   return rank;
                   1949: }
                   1950:
                   1951: int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat)
                   1952: {
                   1953:   int i,j,k,l,inv,a,rank;
                   1954:   unsigned int *t,*pivot,*pk;
                   1955:   unsigned int **mat;
                   1956:
                   1957:   for ( i = 0; i < row; i++ ) rowstat[i] = i;
                   1958:   mat = (unsigned int **)mat0;
                   1959:   for ( rank = 0, j = 0; j < col; j++ ) {
                   1960:     for ( i = rank; i < row; i++ )
                   1961:       mat[i][j] %= md;
                   1962:     for ( i = rank; i < row; i++ )
                   1963:       if ( mat[i][j] )
                   1964:         break;
                   1965:     if ( i == row ) continue;
                   1966:     if ( i != rank ) {
                   1967:       t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                   1968:       k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
                   1969:     }
                   1970:     pivot = mat[rank];
                   1971:     inv = invm(pivot[j],md);
                   1972:     for ( k = j, pk = pivot+k; k < col; k++, pk++ )
                   1973:       if ( *pk ) {
                   1974:         if ( *pk >= (unsigned int)md )
                   1975:           *pk %= md;
                   1976:         DMAR(*pk,inv,0,md,*pk)
                   1977:       }
                   1978:     for ( i = rank+1; i < row; i++ ) {
                   1979:       t = mat[i];
                   1980:       if ( a = t[j] )
                   1981:         red_by_vect(md,t+j,pivot+j,md-a,col-j);
                   1982:     }
                   1983:     rank++;
                   1984:   }
                   1985:   return rank;
                   1986: }
                   1987:
                   1988: int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)
                   1989: {
                   1990:   int i,j,k,l,inv,a,rank;
                   1991:   unsigned int *t,*pivot,*pk;
                   1992:   unsigned int **mat;
                   1993:
                   1994:   mat = (unsigned int **)mat0;
                   1995:   for ( rank = 0, j = 0; j < col; j++ ) {
                   1996:     for ( i = rank; i < row; i++ )
                   1997:       if ( mat[i][j] )
                   1998:         break;
                   1999:     if ( i == row ) {
                   2000:       colstat[j] = 0;
                   2001:       continue;
                   2002:     } else
                   2003:       colstat[j] = 1;
                   2004:     if ( i != rank ) {
                   2005:       t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                   2006:     }
                   2007:     pivot = mat[rank];
                   2008:     inv = _invsf(pivot[j]);
                   2009:     for ( k = j, pk = pivot+k; k < col; k++, pk++ )
                   2010:       if ( *pk )
                   2011:         *pk = _mulsf(*pk,inv);
                   2012:     for ( i = rank+1; i < row; i++ ) {
                   2013:       t = mat[i];
                   2014:       if ( a = t[j] )
                   2015:         red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
                   2016:     }
                   2017:     rank++;
                   2018:   }
                   2019:   for ( j = col-1, l = rank-1; j >= 0; j-- )
                   2020:     if ( colstat[j] ) {
                   2021:       pivot = mat[l];
                   2022:       for ( i = 0; i < l; i++ ) {
                   2023:         t = mat[i];
                   2024:         if ( a = t[j] )
                   2025:           red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
                   2026:       }
                   2027:       l--;
                   2028:     }
                   2029:   return rank;
                   2030: }
                   2031:
                   2032: /* LU decomposition; a[i][i] = 1/U[i][i] */
                   2033:
                   2034: int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
                   2035: {
                   2036:   int row,col;
                   2037:   int i,j,k;
                   2038:   unsigned int *t,*pivot;
                   2039:   unsigned int **a;
                   2040:   unsigned int inv,m;
                   2041:
                   2042:   row = mat->row; col = mat->col;
                   2043:   a = mat->body;
                   2044:   bzero(perm,row*sizeof(int));
                   2045:
                   2046:   for ( i = 0; i < row; i++ )
                   2047:     perm[i] = i;
                   2048:   for ( k = 0; k < col; k++ ) {
                   2049:     for ( i = k; i < row && !a[i][k]; i++ );
                   2050:     if ( i == row )
                   2051:       return 0;
                   2052:     if ( i != k ) {
                   2053:       j = perm[i]; perm[i] = perm[k]; perm[k] = j;
                   2054:       t = a[i]; a[i] = a[k]; a[k] = t;
                   2055:     }
                   2056:     pivot = a[k];
                   2057:     pivot[k] = inv = invm(pivot[k],md);
                   2058:     for ( i = k+1; i < row; i++ ) {
                   2059:       t = a[i];
                   2060:       if ( m = t[k] ) {
                   2061:         DMAR(inv,m,0,md,t[k])
                   2062:         for ( j = k+1, m = md - t[k]; j < col; j++ )
                   2063:           if ( pivot[j] ) {
                   2064:             unsigned int tj;
                   2065:
                   2066:             DMAR(m,pivot[j],t[j],md,tj)
                   2067:             t[j] = tj;
                   2068:           }
                   2069:       }
                   2070:     }
                   2071:   }
                   2072:   return 1;
                   2073: }
                   2074:
                   2075: /*
                   2076:  Input
                   2077:   a: a row x col matrix
                   2078:   md : a modulus
                   2079:
                   2080:  Output:
                   2081:   return : d = the rank of mat
                   2082:   a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])
                   2083:   rinfo: array of length row
                   2084:   cinfo: array of length col
                   2085:     i-th row in new a <-> rinfo[i]-th row in old a
                   2086:   cinfo[j]=1 <=> j-th column is contained in the LU decomp.
                   2087: */
                   2088:
                   2089: int find_lhs_and_lu_mod(unsigned int **a,int row,int col,
                   2090:   unsigned int md,int **rinfo,int **cinfo)
                   2091: {
                   2092:   int i,j,k,d;
                   2093:   int *rp,*cp;
                   2094:   unsigned int *t,*pivot;
                   2095:   unsigned int inv,m;
                   2096:
                   2097:   *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
                   2098:   *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
                   2099:   for ( i = 0; i < row; i++ )
                   2100:     rp[i] = i;
                   2101:   for ( k = 0, d = 0; k < col; k++ ) {
                   2102:     for ( i = d; i < row && !a[i][k]; i++ );
                   2103:     if ( i == row ) {
                   2104:       cp[k] = 0;
                   2105:       continue;
                   2106:     } else
                   2107:       cp[k] = 1;
                   2108:     if ( i != d ) {
                   2109:       j = rp[i]; rp[i] = rp[d]; rp[d] = j;
                   2110:       t = a[i]; a[i] = a[d]; a[d] = t;
                   2111:     }
                   2112:     pivot = a[d];
                   2113:     pivot[k] = inv = invm(pivot[k],md);
                   2114:     for ( i = d+1; i < row; i++ ) {
                   2115:       t = a[i];
                   2116:       if ( m = t[k] ) {
                   2117:         DMAR(inv,m,0,md,t[k])
                   2118:         for ( j = k+1, m = md - t[k]; j < col; j++ )
                   2119:           if ( pivot[j] ) {
                   2120:             unsigned int tj;
                   2121:             DMAR(m,pivot[j],t[j],md,tj)
                   2122:             t[j] = tj;
                   2123:           }
                   2124:       }
                   2125:     }
                   2126:     d++;
                   2127:   }
                   2128:   return d;
                   2129: }
                   2130:
                   2131: int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo)
                   2132: {
                   2133:   int i,j,k;
                   2134:   int *rp;
                   2135:   unsigned int *t,*pivot;
                   2136:   unsigned int inv,m;
                   2137:
                   2138:   *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int));
                   2139:   for ( i = 0; i < n; i++ ) rp[i] = i;
                   2140:   for ( k = 0; k < n; k++ ) {
                   2141:     for ( i = k; i < n && !a[i][k]; i++ );
                   2142:     if ( i == n ) return 0;
                   2143:     if ( i != k ) {
                   2144:       j = rp[i]; rp[i] = rp[k]; rp[k] = j;
                   2145:       t = a[i]; a[i] = a[k]; a[k] = t;
                   2146:     }
                   2147:     pivot = a[k];
                   2148:     inv = invm(pivot[k],md);
                   2149:     for ( i = k+1; i < n; i++ ) {
                   2150:       t = a[i];
                   2151:       if ( m = t[k] ) {
                   2152:         DMAR(inv,m,0,md,t[k])
                   2153:         for ( j = k+1, m = md - t[k]; j < n; j++ )
                   2154:           if ( pivot[j] ) {
                   2155:             unsigned int tj;
                   2156:             DMAR(m,pivot[j],t[j],md,tj)
                   2157:             t[j] = tj;
                   2158:           }
                   2159:       }
                   2160:     }
                   2161:   }
                   2162:   return 1;
                   2163: }
                   2164:
                   2165: /*
                   2166:   Input
                   2167:   a : n x n matrix; a result of LU-decomposition
                   2168:   md : modulus
                   2169:   b : n x l matrix
                   2170:  Output
                   2171:   b = a^(-1)b
                   2172:  */
                   2173:
                   2174: void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize)
                   2175: {
                   2176:   unsigned int *y,*c;
                   2177:   int i,j,k;
                   2178:   unsigned int t,m,m2;
                   2179:
                   2180:   y = (int *)MALLOC_ATOMIC(n*sizeof(int));
                   2181:   c = (int *)MALLOC_ATOMIC(n*sizeof(int));
                   2182:   m2 = md>>1;
                   2183:   for ( k = 0; k < l; k++ ) {
                   2184:     /* copy b[.][k] to c */
                   2185:     for ( i = 0; i < n; i++ )
                   2186:       c[i] = (unsigned int)b[i][k];
                   2187:     /* solve Ly=c */
                   2188:     for ( i = 0; i < n; i++ ) {
                   2189:       for ( t = c[i], j = 0; j < i; j++ )
                   2190:         if ( a[i][j] ) {
                   2191:           m = md - a[i][j];
                   2192:           DMAR(m,y[j],t,md,t)
                   2193:         }
                   2194:       y[i] = t;
                   2195:     }
                   2196:     /* solve Uc=y */
                   2197:     for ( i = n-1; i >= 0; i-- ) {
                   2198:       for ( t = y[i], j =i+1; j < n; j++ )
                   2199:         if ( a[i][j] ) {
                   2200:           m = md - a[i][j];
                   2201:           DMAR(m,c[j],t,md,t)
                   2202:         }
                   2203:       /* a[i][i] = 1/U[i][i] */
                   2204:       DMAR(t,a[i][i],0,md,c[i])
                   2205:     }
                   2206:     /* copy c to b[.][k] with normalization */
                   2207:     if ( normalize )
                   2208:       for ( i = 0; i < n; i++ )
                   2209:         b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
                   2210:     else
                   2211:       for ( i = 0; i < n; i++ )
                   2212:         b[i][k] = c[i];
                   2213:   }
                   2214: }
                   2215:
                   2216: void Pleqm1(NODE arg,VECT *rp)
                   2217: {
                   2218:   MAT m;
                   2219:   VECT vect;
                   2220:   pointer **mat;
                   2221:   Z *v;
                   2222:   Z q;
                   2223:   int **wmat;
                   2224:   int md,i,j,row,col,t,n,status;
                   2225:
                   2226:   asir_assert(ARG0(arg),O_MAT,"leqm1");
                   2227:   asir_assert(ARG1(arg),O_N,"leqm1");
1.2       noro     2228:   m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1       noro     2229:   row = m->row; col = m->col; mat = m->body;
                   2230:   wmat = (int **)almat(row,col);
                   2231:   for ( i = 0; i < row; i++ )
                   2232:     for ( j = 0; j < col; j++ )
                   2233:       wmat[i][j] = remqi((Q)mat[i][j],md);
                   2234:   status = gauss_elim_mod1(wmat,row,col,md);
                   2235:   if ( status < 0 )
                   2236:     *rp = 0;
                   2237:   else if ( status > 0 )
                   2238:     *rp = (VECT)ONE;
                   2239:   else {
                   2240:     n = col - 1;
                   2241:     MKVECT(vect,n);
                   2242:     for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {
1.2       noro     2243:       t = (md-wmat[i][n])%md; STOZ(t,v[i]);
1.1       noro     2244:     }
                   2245:     *rp = vect;
                   2246:   }
                   2247: }
                   2248:
                   2249: int gauss_elim_mod1(int **mat,int row,int col,int md)
                   2250: {
                   2251:   int i,j,k,inv,a,n;
                   2252:   int *t,*pivot;
                   2253:
                   2254:   n = col - 1;
                   2255:   for ( j = 0; j < n; j++ ) {
                   2256:     for ( i = j; i < row && !mat[i][j]; i++ );
                   2257:     if ( i == row )
                   2258:       return 1;
                   2259:     if ( i != j ) {
                   2260:       t = mat[i]; mat[i] = mat[j]; mat[j] = t;
                   2261:     }
                   2262:     pivot = mat[j];
                   2263:     inv = invm(pivot[j],md);
                   2264:     for ( k = j; k <= n; k++ )
                   2265:       pivot[k] = dmar(pivot[k],inv,0,md);
                   2266:     for ( i = j+1; i < row; i++ ) {
                   2267:       t = mat[i];
                   2268:       if ( i != j && (a = t[j]) )
                   2269:         for ( k = j, a = md - a; k <= n; k++ )
                   2270:           t[k] = dmar(pivot[k],a,t[k],md);
                   2271:     }
                   2272:   }
                   2273:   for ( i = n; i < row && !mat[i][n]; i++ );
                   2274:   if ( i == row ) {
                   2275:     for ( j = n-1; j >= 0; j-- ) {
                   2276:       for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
                   2277:         mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
                   2278:         mat[i][j] = 0;
                   2279:       }
                   2280:     }
                   2281:     return 0;
                   2282:   } else
                   2283:     return -1;
                   2284: }
                   2285:
                   2286: void Pgeninvm(NODE arg,LIST *rp)
                   2287: {
                   2288:   MAT m;
                   2289:   pointer **mat;
                   2290:   Z **tmat;
                   2291:   Z q;
                   2292:   unsigned int **wmat;
                   2293:   int md,i,j,row,col,t,status;
                   2294:   MAT mat1,mat2;
                   2295:   NODE node1,node2;
                   2296:
                   2297:   asir_assert(ARG0(arg),O_MAT,"leqm1");
                   2298:   asir_assert(ARG1(arg),O_N,"leqm1");
1.2       noro     2299:   m = (MAT)ARG0(arg); md = ZTOS((Q)ARG1(arg));
1.1       noro     2300:   row = m->row; col = m->col; mat = m->body;
                   2301:   wmat = (unsigned int **)almat(row,col+row);
                   2302:   for ( i = 0; i < row; i++ ) {
                   2303:     bzero((char *)wmat[i],(col+row)*sizeof(int));
                   2304:     for ( j = 0; j < col; j++ )
                   2305:       wmat[i][j] = remqi((Q)mat[i][j],md);
                   2306:   }
                   2307:   status = gauss_elim_geninv_mod(wmat,row,col,md);
                   2308:   if ( status > 0 )
                   2309:     *rp = 0;
                   2310:   else {
                   2311:     MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
                   2312:     for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )
                   2313:       for ( j = 0; j < row; j++ )
1.2       noro     2314:         UTOZ(wmat[i][j+col],tmat[i][j]);
1.1       noro     2315:     for ( tmat = (Z **)mat2->body; i < row; i++ )
                   2316:       for ( j = 0; j < row; j++ )
1.2       noro     2317:         UTOZ(wmat[i][j+col],tmat[i-col][j]);
1.1       noro     2318:      MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
                   2319:   }
                   2320: }
                   2321:
                   2322: int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)
                   2323: {
                   2324:   int i,j,k,inv,a,n,m;
                   2325:   unsigned int *t,*pivot;
                   2326:
                   2327:   n = col; m = row+col;
                   2328:   for ( j = 0; j < n; j++ ) {
                   2329:     for ( i = j; i < row && !mat[i][j]; i++ );
                   2330:     if ( i == row )
                   2331:       return 1;
                   2332:     if ( i != j ) {
                   2333:       t = mat[i]; mat[i] = mat[j]; mat[j] = t;
                   2334:     }
                   2335:     pivot = mat[j];
                   2336:     inv = invm(pivot[j],md);
                   2337:     for ( k = j; k < m; k++ )
                   2338:       pivot[k] = dmar(pivot[k],inv,0,md);
                   2339:     for ( i = j+1; i < row; i++ ) {
                   2340:       t = mat[i];
                   2341:       if ( a = t[j] )
                   2342:         for ( k = j, a = md - a; k < m; k++ )
                   2343:           t[k] = dmar(pivot[k],a,t[k],md);
                   2344:     }
                   2345:   }
                   2346:   for ( j = n-1; j >= 0; j-- ) {
                   2347:     pivot = mat[j];
                   2348:     for ( i = j-1; i >= 0; i-- ) {
                   2349:       t = mat[i];
                   2350:       if ( a = t[j] )
                   2351:         for ( k = j, a = md - a; k < m; k++ )
                   2352:           t[k] = dmar(pivot[k],a,t[k],md);
                   2353:     }
                   2354:   }
                   2355:   return 0;
                   2356: }
                   2357:
                   2358: void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
                   2359: {
                   2360:   GFMMAT lu;
                   2361:   Z *perm,*rhs,*v;
                   2362:   int n,i;
                   2363:   unsigned int md;
                   2364:   unsigned int *b,*sol;
                   2365:   VECT r;
                   2366:
                   2367:   lu = (GFMMAT)ARG0(arg);
                   2368:   perm = (Z *)BDY((VECT)ARG1(arg));
                   2369:   rhs = (Z *)BDY((VECT)ARG2(arg));
1.2       noro     2370:   md = (unsigned int)ZTOS((Z)ARG3(arg));
1.1       noro     2371:   n = lu->col;
                   2372:   b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
                   2373:   sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
                   2374:   for ( i = 0; i < n; i++ )
1.2       noro     2375:     b[i] = ZTOS(rhs[ZTOS(perm[i])]);
1.1       noro     2376:   solve_by_lu_gfmmat(lu,md,b,sol);
                   2377:   MKVECT(r,n);
                   2378:   for ( i = 0, v = (Z *)r->body; i < n; i++ )
1.2       noro     2379:       UTOZ(sol[i],v[i]);
1.1       noro     2380:   *rp = r;
                   2381: }
                   2382:
                   2383: void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
                   2384:   unsigned int *b,unsigned int *x)
                   2385: {
                   2386:   int n;
                   2387:   unsigned int **a;
                   2388:   unsigned int *y;
                   2389:   int i,j;
                   2390:   unsigned int t,m;
                   2391:
                   2392:   n = lu->col;
                   2393:   a = lu->body;
                   2394:   y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
                   2395:   /* solve Ly=b */
                   2396:   for ( i = 0; i < n; i++ ) {
                   2397:     for ( t = b[i], j = 0; j < i; j++ )
                   2398:       if ( a[i][j] ) {
                   2399:         m = md - a[i][j];
                   2400:         DMAR(m,y[j],t,md,t)
                   2401:       }
                   2402:     y[i] = t;
                   2403:   }
                   2404:   /* solve Ux=y */
                   2405:   for ( i = n-1; i >= 0; i-- ) {
                   2406:     for ( t = y[i], j =i+1; j < n; j++ )
                   2407:       if ( a[i][j] ) {
                   2408:         m = md - a[i][j];
                   2409:         DMAR(m,x[j],t,md,t)
                   2410:       }
                   2411:     /* a[i][i] = 1/U[i][i] */
                   2412:     DMAR(t,a[i][i],0,md,x[i])
                   2413:   }
                   2414: }
                   2415:
                   2416: #if 0
                   2417: void Plu_mat(NODE arg,LIST *rp)
                   2418: {
                   2419:   MAT m,lu;
                   2420:   Q dn;
                   2421:   Q *v;
                   2422:   int n,i;
                   2423:   int *iperm;
                   2424:   VECT perm;
                   2425:   NODE n0;
                   2426:
                   2427:   asir_assert(ARG0(arg),O_MAT,"lu_mat");
                   2428:   m = (MAT)ARG0(arg);
                   2429:   n = m->row;
                   2430:   MKMAT(lu,n,n);
                   2431:   lu_dec_cr(m,lu,&dn,&iperm);
                   2432:   MKVECT(perm,n);
                   2433:   for ( i = 0, v = (Q *)perm->body; i < n; i++ )
1.2       noro     2434:     STOZ(iperm[i],v[i]);
1.1       noro     2435:   n0 = mknode(3,lu,dn,perm);
                   2436:   MKLIST(*rp,n0);
                   2437: }
                   2438: #endif
                   2439:
                   2440: void Plu_gfmmat(NODE arg,LIST *rp)
                   2441: {
                   2442:   MAT m;
                   2443:   GFMMAT mm;
                   2444:   unsigned int md;
                   2445:   int i,row,col,status;
                   2446:   int *iperm;
                   2447:   Z *v;
                   2448:   VECT perm;
                   2449:   NODE n0;
                   2450:
                   2451:   asir_assert(ARG0(arg),O_MAT,"lu_gfmmat");
                   2452:   asir_assert(ARG1(arg),O_N,"lu_gfmmat");
1.2       noro     2453:   m = (MAT)ARG0(arg); md = (unsigned int)ZTOS((Q)ARG1(arg));
1.1       noro     2454:   mat_to_gfmmat(m,md,&mm);
                   2455:   row = m->row;
                   2456:   col = m->col;
                   2457:   iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
                   2458:   status = lu_gfmmat(mm,md,iperm);
                   2459:   if ( !status )
                   2460:     n0 = 0;
                   2461:   else {
                   2462:     MKVECT(perm,row);
                   2463:     for ( i = 0, v = (Z *)perm->body; i < row; i++ )
1.2       noro     2464:       STOZ(iperm[i],v[i]);
1.1       noro     2465:     n0 = mknode(2,mm,perm);
                   2466:   }
                   2467:   MKLIST(*rp,n0);
                   2468: }
                   2469:
                   2470: void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
                   2471: {
                   2472:   MAT m;
                   2473:   unsigned int md;
                   2474:
                   2475:   asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
                   2476:   asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1.2       noro     2477:   m = (MAT)ARG0(arg); md = (unsigned int)ZTOS((Q)ARG1(arg));
1.1       noro     2478:   mat_to_gfmmat(m,md,rp);
                   2479: }
                   2480:
                   2481: void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
                   2482: {
                   2483:   unsigned int **wmat;
                   2484:   unsigned int t;
                   2485:   Z **mat;
                   2486:   Z q;
                   2487:   int i,j,row,col;
                   2488:
                   2489:   row = m->row; col = m->col; mat = (Z **)m->body;
                   2490:   wmat = (unsigned int **)almat(row,col);
                   2491:   for ( i = 0; i < row; i++ ) {
                   2492:     bzero((char *)wmat[i],col*sizeof(unsigned int));
                   2493:     for ( j = 0; j < col; j++ )
                   2494:       wmat[i][j] = remqi((Q)mat[i][j],md);
                   2495:   }
                   2496:   TOGFMMAT(row,col,wmat,*rp);
                   2497: }
                   2498:
                   2499: void Pgeninvm_swap(NODE arg,LIST *rp)
                   2500: {
                   2501:   MAT m;
                   2502:   pointer **mat;
                   2503:   Z **tmat;
                   2504:   Z *tvect;
                   2505:   Z q;
                   2506:   unsigned int **wmat,**invmat;
                   2507:   int *index;
                   2508:   unsigned int t,md;
                   2509:   int i,j,row,col,status;
                   2510:   MAT mat1;
                   2511:   VECT vect1;
                   2512:   NODE node1,node2;
                   2513:
                   2514:   asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
                   2515:   asir_assert(ARG1(arg),O_N,"geninvm_swap");
1.2       noro     2516:   m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1       noro     2517:   row = m->row; col = m->col; mat = m->body;
                   2518:   wmat = (unsigned int **)almat(row,col+row);
                   2519:   for ( i = 0; i < row; i++ ) {
                   2520:     bzero((char *)wmat[i],(col+row)*sizeof(int));
                   2521:     for ( j = 0; j < col; j++ )
                   2522:       wmat[i][j] = remqi((Q)mat[i][j],md);
                   2523:     wmat[i][col+i] = 1;
                   2524:   }
                   2525:   status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
                   2526:   if ( status > 0 )
                   2527:     *rp = 0;
                   2528:   else {
                   2529:     MKMAT(mat1,col,col);
                   2530:     for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )
                   2531:       for ( j = 0; j < col; j++ )
1.2       noro     2532:         UTOZ(invmat[i][j],tmat[i][j]);
1.1       noro     2533:     MKVECT(vect1,row);
                   2534:     for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )
1.2       noro     2535:       STOZ(index[i],tvect[i]);
1.1       noro     2536:      MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
                   2537:   }
                   2538: }
                   2539:
                   2540: int gauss_elim_geninv_mod_swap(unsigned int **mat,int row,int col,unsigned int md,
                   2541:     unsigned int ***invmatp,int **indexp)
                   2542: {
                   2543:   int i,j,k,inv,a,n,m;
                   2544:   unsigned int *t,*pivot,*s;
                   2545:   int *index;
                   2546:   unsigned int **invmat;
                   2547:
                   2548:   n = col; m = row+col;
                   2549:   *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
                   2550:   for ( i = 0; i < row; i++ )
                   2551:     index[i] = i;
                   2552:   for ( j = 0; j < n; j++ ) {
                   2553:     for ( i = j; i < row && !mat[i][j]; i++ );
                   2554:     if ( i == row ) {
                   2555:       *indexp = 0; *invmatp = 0; return 1;
                   2556:     }
                   2557:     if ( i != j ) {
                   2558:       t = mat[i]; mat[i] = mat[j]; mat[j] = t;
                   2559:       k = index[i]; index[i] = index[j]; index[j] = k;
                   2560:     }
                   2561:     pivot = mat[j];
                   2562:     inv = (unsigned int)invm(pivot[j],md);
                   2563:     for ( k = j; k < m; k++ )
                   2564:       if ( pivot[k] )
                   2565:         pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
                   2566:     for ( i = j+1; i < row; i++ ) {
                   2567:       t = mat[i];
                   2568:       if ( a = t[j] )
                   2569:         for ( k = j, a = md - a; k < m; k++ )
                   2570:           if ( pivot[k] )
                   2571:             t[k] = dmar(pivot[k],a,t[k],md);
                   2572:     }
                   2573:   }
                   2574:   for ( j = n-1; j >= 0; j-- ) {
                   2575:     pivot = mat[j];
                   2576:     for ( i = j-1; i >= 0; i-- ) {
                   2577:       t = mat[i];
                   2578:       if ( a = t[j] )
                   2579:         for ( k = j, a = md - a; k < m; k++ )
                   2580:           if ( pivot[k] )
                   2581:             t[k] = dmar(pivot[k],a,t[k],md);
                   2582:     }
                   2583:   }
                   2584:   *invmatp = invmat = (unsigned int **)almat(col,col);
                   2585:   for ( i = 0; i < col; i++ )
                   2586:     for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
                   2587:       s[j] = t[col+index[j]];
                   2588:   return 0;
                   2589: }
                   2590:
                   2591: int gauss_elim_geninv_sf_swap(int **mat,int row,int col, int ***invmatp,int **indexp);
                   2592:
                   2593: void Pgeninv_sf_swap(NODE arg,LIST *rp)
                   2594: {
                   2595:   MAT m;
                   2596:   GFS **mat,**tmat;
                   2597:   Z *tvect;
                   2598:   GFS q;
                   2599:   int **wmat,**invmat;
                   2600:   int *index;
                   2601:   unsigned int t;
                   2602:   int i,j,row,col,status;
                   2603:   MAT mat1;
                   2604:   VECT vect1;
                   2605:   NODE node1,node2;
                   2606:
                   2607:   asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap");
                   2608:   m = (MAT)ARG0(arg);
                   2609:   row = m->row; col = m->col; mat = (GFS **)m->body;
                   2610:   wmat = (int **)almat(row,col+row);
                   2611:   for ( i = 0; i < row; i++ ) {
                   2612:     bzero((char *)wmat[i],(col+row)*sizeof(int));
                   2613:     for ( j = 0; j < col; j++ )
                   2614:       if ( q = (GFS)mat[i][j] )
                   2615:         wmat[i][j] = FTOIF(CONT(q));
                   2616:     wmat[i][col+i] = _onesf();
                   2617:   }
                   2618:   status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index);
                   2619:   if ( status > 0 )
                   2620:     *rp = 0;
                   2621:   else {
                   2622:     MKMAT(mat1,col,col);
                   2623:     for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )
                   2624:       for ( j = 0; j < col; j++ )
                   2625:         if ( t = invmat[i][j] ) {
                   2626:           MKGFS(IFTOF(t),tmat[i][j]);
                   2627:         }
                   2628:     MKVECT(vect1,row);
                   2629:     for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )
1.2       noro     2630:       STOZ(index[i],tvect[i]);
1.1       noro     2631:     MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
                   2632:   }
                   2633: }
                   2634:
                   2635: int gauss_elim_geninv_sf_swap(int **mat,int row,int col, int ***invmatp,int **indexp)
                   2636: {
                   2637:   int i,j,k,inv,a,n,m,u;
                   2638:   int *t,*pivot,*s;
                   2639:   int *index;
                   2640:   int **invmat;
                   2641:
                   2642:   n = col; m = row+col;
                   2643:   *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
                   2644:   for ( i = 0; i < row; i++ )
                   2645:     index[i] = i;
                   2646:   for ( j = 0; j < n; j++ ) {
                   2647:     for ( i = j; i < row && !mat[i][j]; i++ );
                   2648:     if ( i == row ) {
                   2649:       *indexp = 0; *invmatp = 0; return 1;
                   2650:     }
                   2651:     if ( i != j ) {
                   2652:       t = mat[i]; mat[i] = mat[j]; mat[j] = t;
                   2653:       k = index[i]; index[i] = index[j]; index[j] = k;
                   2654:     }
                   2655:     pivot = mat[j];
                   2656:     inv = _invsf(pivot[j]);
                   2657:     for ( k = j; k < m; k++ )
                   2658:       if ( pivot[k] )
                   2659:         pivot[k] = _mulsf(pivot[k],inv);
                   2660:     for ( i = j+1; i < row; i++ ) {
                   2661:       t = mat[i];
                   2662:       if ( a = t[j] )
                   2663:         for ( k = j, a = _chsgnsf(a); k < m; k++ )
                   2664:           if ( pivot[k] ) {
                   2665:             u = _mulsf(pivot[k],a);
                   2666:             t[k] = _addsf(u,t[k]);
                   2667:           }
                   2668:     }
                   2669:   }
                   2670:   for ( j = n-1; j >= 0; j-- ) {
                   2671:     pivot = mat[j];
                   2672:     for ( i = j-1; i >= 0; i-- ) {
                   2673:       t = mat[i];
                   2674:       if ( a = t[j] )
                   2675:         for ( k = j, a = _chsgnsf(a); k < m; k++ )
                   2676:           if ( pivot[k] ) {
                   2677:             u = _mulsf(pivot[k],a);
                   2678:             t[k] = _addsf(u,t[k]);
                   2679:           }
                   2680:     }
                   2681:   }
                   2682:   *invmatp = invmat = (int **)almat(col,col);
                   2683:   for ( i = 0; i < col; i++ )
                   2684:     for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
                   2685:       s[j] = t[col+index[j]];
                   2686:   return 0;
                   2687: }
                   2688:
                   2689: void inner_product_int(Z *a,Z *b,int n,Z *r)
                   2690: {
                   2691:   int i;
                   2692:   Z t;
                   2693:
                   2694:   t = 0;
                   2695:   for ( i = 0; i < n; i++ )
                   2696:     muladdtoz(a[i],b[i],&t);
                   2697:   *r = t;
                   2698: }
                   2699:
                   2700: void Pmul_mat_vect_int(NODE arg,VECT *rp)
                   2701: {
                   2702:   MAT mat;
                   2703:   VECT vect,r;
                   2704:   int row,col,i;
                   2705:
                   2706:   mat = (MAT)ARG0(arg);
                   2707:   vect = (VECT)ARG1(arg);
                   2708:   row = mat->row;
                   2709:   col = mat->col;
                   2710:   MKVECT(r,row);
                   2711:   for ( i = 0; i < row; i++ ) {
                   2712:     inner_product_int((Z *)mat->body[i],(Z *)vect->body,col,(Z *)&r->body[i]);
                   2713:   }
                   2714:   *rp = r;
                   2715: }
                   2716:
                   2717: void Pnbpoly_up2(NODE arg,GF2N *rp)
                   2718: {
                   2719:   int m,type,ret;
                   2720:   UP2 r;
                   2721:
1.2       noro     2722:   m = ZTOS((Z)ARG0(arg));
                   2723:   type = ZTOS((Z)ARG1(arg));
1.1       noro     2724:   ret = generate_ONB_polynomial(&r,m,type);
                   2725:   if ( ret == 0 )
                   2726:     MKGF2N(r,*rp);
                   2727:   else
                   2728:     *rp = 0;
                   2729: }
                   2730:
                   2731: void Px962_irredpoly_up2(NODE arg,GF2N *rp)
                   2732: {
                   2733:   int m,ret,w;
                   2734:   GF2N prev;
                   2735:   UP2 r;
                   2736:
1.2       noro     2737:   m = ZTOS((Q)ARG0(arg));
1.1       noro     2738:   prev = (GF2N)ARG1(arg);
                   2739:   if ( !prev ) {
                   2740:     w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
                   2741:     bzero((char *)r->b,w*sizeof(unsigned int));
                   2742:   } else {
                   2743:     r = prev->body;
                   2744:     if ( degup2(r) != m ) {
                   2745:       w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
                   2746:       bzero((char *)r->b,w*sizeof(unsigned int));
                   2747:     }
                   2748:   }
                   2749:   ret = _generate_irreducible_polynomial(r,m);
                   2750:   if ( ret == 0 )
                   2751:     MKGF2N(r,*rp);
                   2752:   else
                   2753:     *rp = 0;
                   2754: }
                   2755:
                   2756: void Pirredpoly_up2(NODE arg,GF2N *rp)
                   2757: {
                   2758:   int m,ret,w;
                   2759:   GF2N prev;
                   2760:   UP2 r;
                   2761:
1.2       noro     2762:   m = ZTOS((Q)ARG0(arg));
1.1       noro     2763:   prev = (GF2N)ARG1(arg);
                   2764:   if ( !prev ) {
                   2765:     w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
                   2766:     bzero((char *)r->b,w*sizeof(unsigned int));
                   2767:   } else {
                   2768:     r = prev->body;
                   2769:     if ( degup2(r) != m ) {
                   2770:       w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
                   2771:       bzero((char *)r->b,w*sizeof(unsigned int));
                   2772:     }
                   2773:   }
                   2774:   ret = _generate_good_irreducible_polynomial(r,m);
                   2775:   if ( ret == 0 )
                   2776:     MKGF2N(r,*rp);
                   2777:   else
                   2778:     *rp = 0;
                   2779: }
                   2780:
                   2781: void Pmat_swap_row_destructive(NODE arg, MAT *m)
                   2782: {
                   2783:   int i1,i2;
                   2784:   pointer *t;
                   2785:   MAT mat;
                   2786:
                   2787:   asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive");
                   2788:   asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");
                   2789:   asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");
                   2790:   mat = (MAT)ARG0(arg);
1.2       noro     2791:   i1 = ZTOS((Q)ARG1(arg));
                   2792:   i2 = ZTOS((Q)ARG2(arg));
1.1       noro     2793:   if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )
                   2794:     error("mat_swap_row_destructive : Out of range");
                   2795:   t = mat->body[i1];
                   2796:   mat->body[i1] = mat->body[i2];
                   2797:   mat->body[i2] = t;
                   2798:   *m = mat;
                   2799: }
                   2800:
                   2801: void Pmat_swap_col_destructive(NODE arg, MAT *m)
                   2802: {
                   2803:   int j1,j2,i,n;
                   2804:   pointer *mi;
                   2805:   pointer t;
                   2806:   MAT mat;
                   2807:
                   2808:   asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive");
                   2809:   asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");
                   2810:   asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");
                   2811:   mat = (MAT)ARG0(arg);
1.2       noro     2812:   j1 = ZTOS((Q)ARG1(arg));
                   2813:   j2 = ZTOS((Q)ARG2(arg));
1.1       noro     2814:   if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )
                   2815:     error("mat_swap_col_destructive : Out of range");
                   2816:   n = mat->row;
                   2817:   for ( i = 0; i < n; i++ ) {
                   2818:     mi = mat->body[i];
                   2819:     t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t;
                   2820:   }
                   2821:   *m = mat;
                   2822: }
                   2823: /*
                   2824:  * f = type 'type' normal polynomial of degree m if exists
                   2825:  * IEEE P1363 A.7.2
                   2826:  *
                   2827:  * return value : 0  --- exists
                   2828:  *                1  --- does not exist
                   2829:  *                -1 --- failure (memory allocation error)
                   2830:  */
                   2831:
                   2832: int generate_ONB_polynomial(UP2 *rp,int m,int type)
                   2833: {
                   2834:   int i,r;
                   2835:   int w;
                   2836:   UP2 f,f0,f1,f2,t;
                   2837:
                   2838:   w = (m>>5)+1;
                   2839:   switch ( type ) {
                   2840:     case 1:
                   2841:       if ( !TypeT_NB_check(m,1) ) return 1;
                   2842:       NEWUP2(f,w); *rp = f; f->w = w;
                   2843:       /* set all the bits */
                   2844:       for ( i = 0; i < w; i++ )
                   2845:         f->b[i] = 0xffffffff;
                   2846:       /* mask the top word if necessary */
                   2847:       if ( r = (m+1)&31 )
                   2848:         f->b[w-1] &= (1<<r)-1;
                   2849:       return 0;
                   2850:       break;
                   2851:     case 2:
                   2852:       if ( !TypeT_NB_check(m,2) ) return 1;
                   2853:       NEWUP2(f,w); *rp = f;
                   2854:       W_NEWUP2(f0,w);
                   2855:       W_NEWUP2(f1,w);
                   2856:       W_NEWUP2(f2,w);
                   2857:
                   2858:       /* recursion for genrating Type II normal polynomial */
                   2859:
                   2860:       /* f0 = 1, f1 = t+1 */
                   2861:       f0->w = 1; f0->b[0] = 1;
                   2862:       f1->w = 1; f1->b[0] = 3;
                   2863:       for ( i = 2; i <= m; i++ ) {
                   2864:         /* f2 = t*f1+f0 */
                   2865:         _bshiftup2(f1,-1,f2);
                   2866:         _addup2_destructive(f2,f0);
                   2867:         /* cyclic change of the variables */
                   2868:         t = f0; f0 = f1; f1 = f2; f2 = t;
                   2869:       }
                   2870:       _copyup2(f1,f);
                   2871:       return 0;
                   2872:       break;
                   2873:     default:
                   2874:       return -1;
                   2875:       break;
                   2876:     }
                   2877: }
                   2878:
                   2879: /*
                   2880:  * f = an irreducible trinomial or pentanomial of degree d 'after' f
                   2881:  * return value : 0  --- exists
                   2882:  *                1  --- does not exist (exhaustion)
                   2883:  */
                   2884:
                   2885: int _generate_irreducible_polynomial(UP2 f,int d)
                   2886: {
                   2887:   int ret,i,j,k,nz,i0,j0,k0;
                   2888:   int w;
                   2889:   unsigned int *fd;
                   2890:
                   2891:   /*
                   2892:    * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
                   2893:    * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
                   2894:    * otherwise i0,j0,k0 is set to 0.
                   2895:    */
                   2896:
                   2897:   fd = f->b;
                   2898:   w = (d>>5)+1;
                   2899:   if ( f->w && (d==degup2(f)) ) {
                   2900:     for ( nz = 0, i = d; i >= 0; i-- )
                   2901:       if ( fd[i>>5]&(1<<(i&31)) ) nz++;
                   2902:     switch ( nz ) {
                   2903:       case 3:
                   2904:         for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
                   2905:         /* reset i0-th bit */
                   2906:         fd[i0>>5] &= ~(1<<(i0&31));
                   2907:         j0 = k0 = 0;
                   2908:         break;
                   2909:       case 5:
                   2910:         for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
                   2911:         /* reset i0-th bit */
                   2912:         fd[i0>>5] &= ~(1<<(i0&31));
                   2913:         for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
                   2914:         /* reset j0-th bit */
                   2915:         fd[j0>>5] &= ~(1<<(j0&31));
                   2916:         for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
                   2917:         /* reset k0-th bit */
                   2918:         fd[k0>>5] &= ~(1<<(k0&31));
                   2919:         break;
                   2920:       default:
                   2921:         f->w = 0; break;
                   2922:     }
                   2923:   } else
                   2924:     f->w = 0;
                   2925:
                   2926:   if ( !f->w ) {
                   2927:     fd = f->b;
                   2928:     f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
                   2929:     i0 = j0 = k0 = 0;
                   2930:   }
                   2931:   /* if j0 > 0 then f is already a pentanomial */
                   2932:   if ( j0 > 0 ) goto PENTA;
                   2933:
                   2934:   /* searching for an irreducible trinomial */
                   2935:
                   2936:   for ( i = 1; 2*i <= d; i++ ) {
                   2937:     /* skip the polynomials 'before' f */
                   2938:     if ( i < i0 ) continue;
                   2939:     if ( i == i0 ) { i0 = 0; continue; }
                   2940:     /* set i-th bit */
                   2941:     fd[i>>5] |= (1<<(i&31));
                   2942:     ret = irredcheck_dddup2(f);
                   2943:     if ( ret == 1 ) return 0;
                   2944:     /* reset i-th bit */
                   2945:     fd[i>>5] &= ~(1<<(i&31));
                   2946:   }
                   2947:
                   2948:   /* searching for an irreducible pentanomial */
                   2949: PENTA:
                   2950:   for ( i = 1; i < d; i++ ) {
                   2951:     /* skip the polynomials 'before' f */
                   2952:     if ( i < i0 ) continue;
                   2953:     if ( i == i0 ) i0 = 0;
                   2954:     /* set i-th bit */
                   2955:     fd[i>>5] |= (1<<(i&31));
                   2956:     for ( j = i+1; j < d; j++ ) {
                   2957:       /* skip the polynomials 'before' f */
                   2958:       if ( j < j0 ) continue;
                   2959:       if ( j == j0 ) j0 = 0;
                   2960:       /* set j-th bit */
                   2961:       fd[j>>5] |= (1<<(j&31));
                   2962:       for ( k = j+1; k < d; k++ ) {
                   2963:         /* skip the polynomials 'before' f */
                   2964:         if ( k < k0 ) continue;
                   2965:         else if ( k == k0 ) { k0 = 0; continue; }
                   2966:         /* set k-th bit */
                   2967:         fd[k>>5] |= (1<<(k&31));
                   2968:         ret = irredcheck_dddup2(f);
                   2969:         if ( ret == 1 ) return 0;
                   2970:         /* reset k-th bit */
                   2971:         fd[k>>5] &= ~(1<<(k&31));
                   2972:       }
                   2973:       /* reset j-th bit */
                   2974:       fd[j>>5] &= ~(1<<(j&31));
                   2975:     }
                   2976:     /* reset i-th bit */
                   2977:     fd[i>>5] &= ~(1<<(i&31));
                   2978:   }
                   2979:   /* exhausted */
                   2980:   return 1;
                   2981: }
                   2982:
                   2983: /*
                   2984:  * f = an irreducible trinomial or pentanomial of degree d 'after' f
                   2985:  *
                   2986:  * searching strategy:
                   2987:  *   trinomial x^d+x^i+1:
                   2988:  *         i is as small as possible.
                   2989:  *   trinomial x^d+x^i+x^j+x^k+1:
                   2990:  *         i is as small as possible.
                   2991:  *         For such i, j is as small as possible.
                   2992:  *         For such i and j, 'k' is as small as possible.
                   2993:  *
                   2994:  * return value : 0  --- exists
                   2995:  *                1  --- does not exist (exhaustion)
                   2996:  */
                   2997:
                   2998: int _generate_good_irreducible_polynomial(UP2 f,int d)
                   2999: {
                   3000:   int ret,i,j,k,nz,i0,j0,k0;
                   3001:   int w;
                   3002:   unsigned int *fd;
                   3003:
                   3004:   /*
                   3005:    * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
                   3006:    * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
                   3007:    * otherwise i0,j0,k0 is set to 0.
                   3008:    */
                   3009:
                   3010:   fd = f->b;
                   3011:   w = (d>>5)+1;
                   3012:   if ( f->w && (d==degup2(f)) ) {
                   3013:     for ( nz = 0, i = d; i >= 0; i-- )
                   3014:       if ( fd[i>>5]&(1<<(i&31)) ) nz++;
                   3015:     switch ( nz ) {
                   3016:       case 3:
                   3017:         for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
                   3018:         /* reset i0-th bit */
                   3019:         fd[i0>>5] &= ~(1<<(i0&31));
                   3020:         j0 = k0 = 0;
                   3021:         break;
                   3022:       case 5:
                   3023:         for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
                   3024:         /* reset i0-th bit */
                   3025:         fd[i0>>5] &= ~(1<<(i0&31));
                   3026:         for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
                   3027:         /* reset j0-th bit */
                   3028:         fd[j0>>5] &= ~(1<<(j0&31));
                   3029:         for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
                   3030:         /* reset k0-th bit */
                   3031:         fd[k0>>5] &= ~(1<<(k0&31));
                   3032:         break;
                   3033:       default:
                   3034:         f->w = 0; break;
                   3035:     }
                   3036:   } else
                   3037:     f->w = 0;
                   3038:
                   3039:   if ( !f->w ) {
                   3040:     fd = f->b;
                   3041:     f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
                   3042:     i0 = j0 = k0 = 0;
                   3043:   }
                   3044:   /* if j0 > 0 then f is already a pentanomial */
                   3045:   if ( j0 > 0 ) goto PENTA;
                   3046:
                   3047:   /* searching for an irreducible trinomial */
                   3048:
                   3049:   for ( i = 1; 2*i <= d; i++ ) {
                   3050:     /* skip the polynomials 'before' f */
                   3051:     if ( i < i0 ) continue;
                   3052:     if ( i == i0 ) { i0 = 0; continue; }
                   3053:     /* set i-th bit */
                   3054:     fd[i>>5] |= (1<<(i&31));
                   3055:     ret = irredcheck_dddup2(f);
                   3056:     if ( ret == 1 ) return 0;
                   3057:     /* reset i-th bit */
                   3058:     fd[i>>5] &= ~(1<<(i&31));
                   3059:   }
                   3060:
                   3061:   /* searching for an irreducible pentanomial */
                   3062: PENTA:
                   3063:   for ( i = 3; i < d; i++ ) {
                   3064:     /* skip the polynomials 'before' f */
                   3065:     if ( i < i0 ) continue;
                   3066:     if ( i == i0 ) i0 = 0;
                   3067:     /* set i-th bit */
                   3068:     fd[i>>5] |= (1<<(i&31));
                   3069:      for ( j = 2; j < i; j++ ) {
                   3070:       /* skip the polynomials 'before' f */
                   3071:       if ( j < j0 ) continue;
                   3072:       if ( j == j0 ) j0 = 0;
                   3073:       /* set j-th bit */
                   3074:       fd[j>>5] |= (1<<(j&31));
                   3075:        for ( k = 1; k < j; k++ ) {
                   3076:         /* skip the polynomials 'before' f */
                   3077:         if ( k < k0 ) continue;
                   3078:         else if ( k == k0 ) { k0 = 0; continue; }
                   3079:         /* set k-th bit */
                   3080:         fd[k>>5] |= (1<<(k&31));
                   3081:         ret = irredcheck_dddup2(f);
                   3082:         if ( ret == 1 ) return 0;
                   3083:         /* reset k-th bit */
                   3084:         fd[k>>5] &= ~(1<<(k&31));
                   3085:       }
                   3086:       /* reset j-th bit */
                   3087:       fd[j>>5] &= ~(1<<(j&31));
                   3088:     }
                   3089:     /* reset i-th bit */
                   3090:     fd[i>>5] &= ~(1<<(i&31));
                   3091:   }
                   3092:   /* exhausted */
                   3093:   return 1;
                   3094: }
                   3095:
                   3096: void printqmat(Q **mat,int row,int col)
                   3097: {
                   3098:   int i,j;
                   3099:
                   3100:   for ( i = 0; i < row; i++ ) {
                   3101:     for ( j = 0; j < col; j++ ) {
                   3102:       printnum((Num)mat[i][j]); printf(" ");
                   3103:     }
                   3104:     printf("\n");
                   3105:   }
                   3106: }
                   3107:
                   3108: void printimat(int **mat,int row,int col)
                   3109: {
                   3110:   int i,j;
                   3111:
                   3112:   for ( i = 0; i < row; i++ ) {
                   3113:     for ( j = 0; j < col; j++ ) {
                   3114:       printf("%d ",mat[i][j]);
                   3115:     }
                   3116:     printf("\n");
                   3117:   }
                   3118: }
                   3119:
                   3120: void Pnd_det(NODE arg,P *rp)
                   3121: {
                   3122:   if ( argc(arg) == 1 )
                   3123:     nd_det(0,ARG0(arg),rp);
                   3124:   else
1.2       noro     3125:     nd_det(ZTOS((Q)ARG1(arg)),ARG0(arg),rp);
1.1       noro     3126: }
                   3127:
                   3128: void Pmat_col(NODE arg,VECT *rp)
                   3129: {
                   3130:   int i,j,n;
                   3131:   MAT mat;
                   3132:   VECT vect;
                   3133:
                   3134:   asir_assert(ARG0(arg),O_MAT,"mat_col");
                   3135:   asir_assert(ARG1(arg),O_N,"mat_col");
                   3136:   mat = (MAT)ARG0(arg);
1.2       noro     3137:   j = ZTOS((Q)ARG1(arg));
1.1       noro     3138:   if ( j < 0 || j >= mat->col) {
                   3139:     error("mat_col : Out of range");
                   3140:   }
                   3141:   n = mat->row;
                   3142:   MKVECT(vect,n);
                   3143:   for(i=0; i<n; i++) {
                   3144:     BDY(vect)[i] = BDY(mat)[i][j];
                   3145:   }
                   3146:   *rp = vect;
                   3147: }
                   3148:
                   3149: NODE triangleq(NODE e)
                   3150: {
                   3151:   int n,i,k;
                   3152:   V v;
                   3153:   VL vl;
                   3154:   P *p;
                   3155:   NODE r,r1;
                   3156:
                   3157:   n = length(e);
                   3158:   p = (P *)MALLOC(n*sizeof(P));
                   3159:   for ( i = 0; i < n; i++, e = NEXT(e) ) p[i] = (P)BDY(e);
                   3160:   i = 0;
                   3161:   while ( 1 ) {
                   3162:     for ( ; i < n && !p[i]; i++ );
                   3163:     if ( i == n ) break;
                   3164:     if ( OID(p[i]) == O_N ) return 0;
                   3165:     v = p[i]->v;
                   3166:     for ( k = i+1; k < n; k++ )
                   3167:       if ( p[k] ) {
                   3168:         if ( OID(p[k]) == O_N ) return 0;
                   3169:         if ( p[k]->v == v ) p[k] = 0;
                   3170:       }
                   3171:     i++;
                   3172:   }
                   3173:   for ( r = 0, i = 0; i < n; i++ ) {
                   3174:     if ( p[i] ) {
                   3175:       MKNODE(r1,p[i],r); r = r1;
                   3176:     }
                   3177:   }
                   3178:   return r;
                   3179: }
                   3180:
                   3181: void Ptriangleq(NODE arg,LIST *rp)
                   3182: {
                   3183:   NODE ret;
                   3184:
                   3185:   asir_assert(ARG0(arg),O_LIST,"sparseleq");
                   3186:   ret = triangleq(BDY((LIST)ARG0(arg)));
                   3187:   MKLIST(*rp,ret);
                   3188: }

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