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

Annotation of OpenXM_contrib2/asir2000/builtin/array.c, Revision 1.41

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

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