[BACK]Return to weight CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / lib

Annotation of OpenXM_contrib2/asir2000/lib/weight, Revision 1.32

1.4       kimura      1: load("solve")$
                      2: load("gr")$
                      3:
1.18      kimura      4: #define EPS 1E-6
                      5: #define TINY 1E-20
                      6: #define MAX_ITER 100
                      7:
                      8: def rotate(A,I,J,K,L,C,S){
                      9:
                     10:        X=A[I][J];
                     11:        Y=A[K][L];
                     12:        A[I][J]=X*C-Y*S;
                     13:        A[K][L]=X*S+Y*C;
                     14:
                     15:        return 1;
                     16: }
                     17:
                     18: def jacobi(N,A,W){
                     19:
                     20:        S=OFFDIAG=0.0;
                     21:
                     22:        for(J=0;J<N;J++){
                     23:
                     24:                for(K=0;K<N;K++)
                     25:                        W[J][K]=0.0;
                     26:
                     27:                W[J][J]=1.0;
                     28:                S+=A[J][J]*A[J][J];
                     29:
                     30:                for(K=J+1;K<N;K++)
                     31:                        OFFDIAG+=A[J][K]*A[J][K];
                     32:        }
                     33:
                     34:        TOLERANCE=EPS*EPS*(S/2+OFFDIAG);
                     35:
                     36:        for(ITER=1;ITER<=MAX_ITER;ITER++){
                     37:
                     38:                OFFDIAG=0.0;
                     39:                for(J=0;J<N-1;J++)
                     40:                        for(K=J+1;K<N;K++)
                     41:                                OFFDIAG+=A[J][K]*A[J][K];
                     42:
                     43:                if(OFFDIAG < TOLERANCE)
                     44:                        break;
                     45:
                     46:                for(J=0;J<N-1;J++){
                     47:                        for(K=J+1;K<N;K++){
                     48:
                     49:                                if(dabs(A[J][K])<TINY)
                     50:                                        continue;
                     51:
                     52:                                T=(A[K][K]-A[J][J])/(2.0*A[J][K]);
                     53:
                     54:                                if(T>=0.0)
                     55:                                        T=1.0/(T+dsqrt(T*T+1));
                     56:                                else
                     57:                                        T=1.0/(T-dsqrt(T*T+1));
                     58:
                     59:                                C=1.0/dsqrt(T*T+1);
                     60:
                     61:                                S=T*C;
                     62:
                     63:                                T*=A[J][K];
                     64:
                     65:                                A[J][J]-=T;
                     66:                                A[K][K]+=T;
                     67:                                A[J][K]=0.0;
                     68:
                     69:                                for(I=0;I<J;I++)
                     70:                                        rotate(A,I,J,I,K,C,S);
                     71:
                     72:                                for(I=J+1;I<K;I++)
                     73:                                        rotate(A,J,I,I,K,C,S);
                     74:
                     75:                                for(I=K+1;I<N;I++)
                     76:                                        rotate(A,J,I,K,I,C,S);
                     77:
                     78:                                for(I=0;I<N;I++)
                     79:                                        rotate(W,J,I,K,I,C,S);
                     80:
                     81:                        }
                     82:                }
                     83:        }
                     84:
                     85:        if (ITER > MAX_ITER)
                     86:                return 0;
                     87:
                     88:        for(I=0;I<N-1;I++){
                     89:
                     90:                K=I;
                     91:
                     92:                T=A[K][K];
                     93:
                     94:                for(J=I+1;J<N;J++)
                     95:                        if(A[J][J]>T){
                     96:                                K=J;
                     97:                                T=A[K][K];
                     98:                        }
                     99:
                    100:                A[K][K]=A[I][I];
                    101:
                    102:                A[I][I]=T;
                    103:
                    104:                V=W[K];
                    105:
                    106:                W[K]=W[I];
                    107:
                    108:                W[I]=V;
                    109:        }
                    110:
                    111:        return 1;
                    112: }
                    113:
1.24      kimura    114: def interval2value(A,Vars){
                    115:
                    116:        B=atl(A)$
                    117:
                    118:        if(length(B)>2){
                    119:                print("bug")$
                    120:                return []$
                    121:        }
1.26      kimura    122:        else if(length(B)==0){
                    123:                if(fop(A)==0)
1.24      kimura    124:                        return [Vars,1]$
                    125:                else
                    126:                        return []$
                    127:        }
                    128:        else if(length(B)==1){
                    129:
                    130:                C=fargs(B[0])$
                    131:                D=vars(C)$
                    132:                E=solve(C,D)$
                    133:
                    134:                if(fop(B[0])==15)
                    135:                        return [Vars,E[0][1]+1]$
                    136:                else if(fop(B[0])==11)
                    137:                        return [Vars,E[0][1]-1]$
1.25      kimura    138:                else if(fop(B[0])==8)
                    139:                        return [Vars,E[0][1]]$
1.24      kimura    140:                else
                    141:                        return []$
                    142:        }
                    143:        else{
                    144:
                    145:                C=fargs(B[0])$
                    146:                D=vars(C)$
                    147:                E=solve(C,D)$
                    148:
                    149:                C=fargs(B[1])$
                    150:                D=vars(C)$
                    151:                F=solve(C,D)$
                    152:
                    153:                return [Vars,(E[0][1]+F[0][1])/2]$
                    154:        }
                    155:
                    156: }
                    157:
                    158: def fixpointmain(F,Vars){
                    159:
                    160:        RET=[]$
                    161:        for(I=length(Vars)-1;I>=1;I--){
                    162:
                    163:                for(H=[],J=0;J<I;J++)
                    164:                        H=cons(Vars[J],H)$
                    165:
                    166:                G=interval2value(qe(ex(H,F)),Vars[I])$
                    167:
                    168:                if(G==[])
                    169:                        return RET$
                    170:                else
                    171:                        RET=cons(G,RET)$
                    172:
                    173:                F=subf(F,G[0],G[1])$
                    174:        }
                    175:
                    176:        G=interval2value(simpl(F),Vars[0])$
                    177:
                    178:        if(G==[])
                    179:                return RET$
                    180:        else
                    181:                RET=cons(G,RET)$
                    182:
                    183:        return RET$
                    184: }
                    185:
                    186:
                    187: def fixedpoint(A,FLAG){
                    188:
                    189:        Vars=vars(A)$
                    190:
                    191:        N=length(A)$
                    192:
                    193:        if (FLAG==0)
                    194:                for(F=@true,I=0;I < N; I++ ) { F = F @&& A[I] @> 0$ }
                    195:        else if (FLAG==1)
                    196:                for(F=@true,I=0;I < N; I++ ) { F = F @&& A[I] @< 0$ }
                    197:
                    198:        return fixpointmain(F,Vars)$
                    199: }
                    200:
1.8       kimura    201: def junban(A,B){
                    202:        return (A<B ? 1:(A>B ? -1:0))$
                    203: }
                    204:
1.10      kimura    205: def bsort(A){
1.8       kimura    206:
1.10      kimura    207:        K=size(A)[0]-1$
                    208:        while(K>=0){
                    209:                J=-1$
                    210:                for(I=1;I<=K;I++)
                    211:                        if(A[I-1][0]<A[I][0]){
                    212:                                J=I-1$
                    213:                                X=A[J]$
                    214:                                A[J]=A[I]$
                    215:                                A[I]=X$
                    216:                        }
                    217:                K=J$
                    218:        }
                    219:        return A$
                    220: }
                    221:
1.26      kimura    222: def wsort(A,B,C,ID){
1.10      kimura    223:
1.26      kimura    224:        D=newvect(length(B))$
                    225:        for(I=0;I<length(B);I++)
                    226:                D[I]=[A[I],B[I],C[I]]$
1.10      kimura    227:
1.26      kimura    228:        D=bsort(D)$
1.10      kimura    229:
1.26      kimura    230:        for(E=[],I=0;I<length(B);I++)
                    231:                E=cons(D[I][1],E)$
                    232:        E=reverse(E)$
1.10      kimura    233:
1.26      kimura    234:        for(F=[],I=0;I<length(B);I++)
                    235:                F=cons(D[I][2],F)$
                    236:        F=reverse(F)$
1.8       kimura    237:
1.26      kimura    238:        return [[ID,E,F]]$
1.8       kimura    239: }
                    240:
1.4       kimura    241: def nonposdegchk(Res){
                    242:
                    243:        for(I=0;I<length(Res);I++)
                    244:                if(Res[I][1]<=0)
                    245:                        return 0$
                    246:
                    247:        return 1$
                    248: }
                    249:
1.8       kimura    250: def getgcd(A,B){
                    251:
1.31      kimura    252:        Anum=length(A)$
1.8       kimura    253:
1.31      kimura    254:        TMP=[]$
                    255:        for(I=0;I<length(B);I++){
1.8       kimura    256:
1.31      kimura    257:                for(J=0;J<Anum;J++)
                    258:                        if(B[I]==A[J][0])
                    259:                                break;
1.8       kimura    260:
1.31      kimura    261:                if(J==Anum)
                    262:                        TMP=cons([B[I],B[I]],TMP)$
                    263:        }
                    264:        A=append(A,TMP)$
1.8       kimura    265:
1.31      kimura    266:        Anum=length(A)$
                    267:        A=map(ltov,A)$
1.8       kimura    268:
1.31      kimura    269:        for(D=0,I=0;I<Anum;I++)
                    270:                D=gcd(D,A[I][1])$
1.8       kimura    271:
                    272:        if(D!=0){
1.31      kimura    273:                for(I=0;I<Anum;I++)
                    274:                        A[I][1]=red(A[I][1]/D)$
1.8       kimura    275:        }
1.6       kimura    276:
1.31      kimura    277:        for(L=1,D=0,I=0;I<Anum;I++){
                    278:                if(type(TMP=dn(A[I][1]))==1)
1.8       kimura    279:                        L=ilcm(L,TMP)$
                    280:
1.31      kimura    281:                if(type(TMP=nm(A[I][1]))==1)
1.8       kimura    282:                        D=igcd(D,TMP)$
                    283:        }
                    284:
1.31      kimura    285:        for(I=0;I<Anum;I++)
                    286:                A[I][1]=A[I][1]*L$
                    287:
                    288:        if(D!=0){
1.6       kimura    289:
1.31      kimura    290:                for(I=0;I<Anum;I++)
                    291:                        A[I][1]=A[I][1]/D$
                    292:        }
1.6       kimura    293:
1.31      kimura    294:        return map(vtol,A)$
1.6       kimura    295: }
                    296:
1.10      kimura    297: def makeret(Res,Vars,FLAG){
1.4       kimura    298:
1.8       kimura    299:        ResNum=length(Res)$
1.4       kimura    300:        VarsNum=length(Vars)$
                    301:
1.17      kimura    302:        ResVec=newvect(ResNum)$
1.4       kimura    303:
1.29      kimura    304:        if(FLAG)
                    305:                M=0$
                    306:        else
                    307:                M=-1$
                    308:
                    309:         for(I=0;I<ResNum;I++){
1.17      kimura    310:                 if(member(Res[I][0],Vars)){
                    311:                         ResVec[I]=Res[I][1]$
1.16      kimura    312:
1.29      kimura    313:                         if(FLAG){
                    314:                                if(type(ResVec[I])==1){
                    315:                                        if(M==0)
                    316:                                                M=ResVec[I]$
                    317:                                        else
                    318:                                                if(ResVec[I]<M)
                    319:                                                        M=ResVec[I]$
                    320:                                }
                    321:                                else
                    322:                                        M=-1$
                    323:                        }
                    324:                }
1.17      kimura    325:         }
1.16      kimura    326:
1.29      kimura    327:
                    328:        if(M!=-1)
1.8       kimura    329:                ResVec=ResVec/M;
1.4       kimura    330:
1.17      kimura    331:        RET=newvect(VarsNum,Vars)$
1.4       kimura    332:
1.8       kimura    333:        for(I=0;I<ResNum;I++){
                    334:                for(J=0;J<VarsNum;J++)
                    335:                        if(Vars[J]==Res[I][0])
                    336:                                break$
1.4       kimura    337:
1.8       kimura    338:                if(J<VarsNum)
                    339:                        RET[J]=ResVec[I]$
                    340:        }
                    341:
1.4       kimura    342:        for(I=0;I<VarsNum;I++)
1.8       kimura    343:                if(type(RET[I])!=1)
                    344:                        return [1,RET]$
1.4       kimura    345:
1.8       kimura    346:        return [0,RET]$
1.4       kimura    347: }
                    348:
                    349: def roundret(V){
                    350:
1.8       kimura    351:        VN=size(V)[0]$
1.4       kimura    352:
1.32    ! kimura    353:        K=1$
1.8       kimura    354:        RET0=V$
1.32    ! kimura    355:        RET1=map(drint,RET0)$
        !           356:        S=0$
        !           357:        for(J=0;J<VN;J++)
        !           358:                S+=(RET0[J]-RET1[J])^2$
        !           359:
        !           360:        for(I=2;I<10;I++){
        !           361:                RET0=I*V$
        !           362:                RET1=map(drint,RET0)$
        !           363:
        !           364:                T=0$
        !           365:                for(J=0;J<VN;J++)
        !           366:                        T+=(RET0[J]-RET1[J])^2$
        !           367:
        !           368:                if(T<S){
        !           369:                        K=I$
        !           370:                        S=T$
        !           371:                }
1.4       kimura    372:        }
                    373:
1.32    ! kimura    374:        return map(drint,K*V)$
1.4       kimura    375: }
                    376:
                    377: def chkou(L,ExpMat,CHAGORD){
                    378:
1.8       kimura    379:        for(P=1,I=0;I<L;I++){
1.4       kimura    380:                Q=ExpMat[L][CHAGORD[I]]$
                    381:                for(J=0;J<size(ExpMat[0])[0];J++){
                    382:                        ExpMat[L][CHAGORD[J]]=red((ExpMat[I][CHAGORD[I]]
                    383:                                *ExpMat[L][CHAGORD[J]]-
                    384:                                        Q*ExpMat[I][CHAGORD[J]])/P)$
                    385:                }
                    386:
                    387:                P=ExpMat[I][CHAGORD[I]]$
                    388:        }
                    389:
                    390:        for(J=0;J<size(ExpMat[0])[0];J++)
                    391:                if(ExpMat[L][CHAGORD[J]]!=0)
                    392:                        break$
                    393:
                    394:        if(J==size(ExpMat[0])[0])
                    395:                return L$
                    396:        else{
                    397:                TMP=CHAGORD[L]$
                    398:                CHAGORD[L]=CHAGORD[J]$
                    399:                CHAGORD[J]=TMP$
                    400:                return (L+1)$
                    401:        }
                    402: }
                    403:
1.8       kimura    404: def qcheckmain(PolyList,Vars){
1.4       kimura    405:
                    406:        RET=[]$
                    407:        PolyListNum=length(PolyList)$
                    408:        VarsNum=length(Vars)$
                    409:
                    410:        ExpMat=newvect(VarsNum)$
                    411:        CHAGORD=newvect(VarsNum)$
                    412:        for(I=0;I<VarsNum;I++)
                    413:                CHAGORD[I]=I$
                    414:
                    415:        L=0$
                    416:        for(I=0;I<PolyListNum;I++){
                    417:                Poly=dp_ptod(PolyList[I],Vars)$
                    418:                BASE0=dp_etov(dp_ht(Poly))$
                    419:                Poly=dp_rest(Poly)$
                    420:                for(;Poly!=0;Poly=dp_rest(Poly)){
                    421:                        ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
                    422:                        L=chkou(L,ExpMat,CHAGORD)$
1.8       kimura    423:                        if(L==VarsNum-1)
                    424:                                return [L,CHAGORD,ExpMat]$
1.4       kimura    425:                }
                    426:        }
                    427:
1.8       kimura    428:        return [L,CHAGORD,ExpMat]$
1.4       kimura    429: }
                    430:
                    431: def inner(A,B){
                    432:
                    433:        SUM=0$
                    434:        for(I=0;I<size(A)[0];I++)
                    435:                SUM+=A[I]*B[I]$
                    436:
                    437:        return SUM$
                    438: }
                    439:
                    440: def checktd(PolyList,Vars,ResVars){
                    441:
                    442:        PolyListNum=length(PolyList)$
                    443:        VarsNum=length(Vars)$
                    444:
                    445:        L=0$
                    446:        for(I=0;I<PolyListNum;I++){
                    447:                Poly=dp_ptod(PolyList[I],Vars)$
                    448:                J0=inner(dp_etov(dp_ht(Poly)),ResVars)$
                    449:                Poly=dp_rest(Poly)$
                    450:                for(;Poly!=0;Poly=dp_rest(Poly))
                    451:                        if(J0!=inner(dp_etov(dp_ht(Poly)),ResVars))
                    452:                                return 0$
                    453:        }
                    454:
                    455:        return 1$
                    456: }
                    457:
1.24      kimura    458: def value2(Vars,Ans){
                    459:
                    460:        N=length(Vars)$
                    461:        Res=newvect(N)$
                    462:        for(I=0;I<N;I++){
                    463:                Res[I]=newvect(2)$
                    464:                Res[I][0]=Vars[I]$
                    465:                Res[I][1]=Ans[I]$
                    466:        }
1.31      kimura    467:        Res=map(vtol,Res)$
                    468:        Res=vtol(Res)$
1.24      kimura    469:
                    470:        Res=getgcd(Res,Vars)$
                    471:
                    472:        if(nonposdegchk(Res)){
                    473:                TMP1=makeret(Res,Vars,1)$
                    474:                return vtol(TMP1[1])$
                    475:        }
                    476:        else
                    477:                return []$
                    478: }
                    479:
1.10      kimura    480: def qcheck(PolyList,Vars,FLAG){
1.4       kimura    481:
1.10      kimura    482:        RET=[]$
1.8       kimura    483:        Res=qcheckmain(PolyList,Vars)$
1.4       kimura    484:        VarsNum=length(Vars)$
                    485:
                    486:        IndNum=Res[0]$
                    487:        CHAGORD=Res[1]$
                    488:        ExpMat=Res[2]$
                    489:
                    490:        SolveList=[]$
                    491:        for(I=0;I<IndNum;I++){
                    492:                TMP=0$
                    493:                for(J=0;J<VarsNum;J++)
                    494:                        TMP+=ExpMat[I][CHAGORD[J]]*Vars[CHAGORD[J]]$
                    495:
                    496:                SolveList=cons(TMP,SolveList)$
                    497:        }
                    498:
1.8       kimura    499:        Rea=vars(SolveList)$
                    500:
1.4       kimura    501:        VarsList=[]$
                    502:        for(I=0;I<VarsNum;I++)
1.31      kimura    503:                if(member(TMP0=Vars[CHAGORD[I]],Rea))
1.8       kimura    504:                        VarsList=cons(Vars[CHAGORD[I]],VarsList)$
1.4       kimura    505:
                    506:        Res=solve(reverse(SolveList),reverse(VarsList))$
1.8       kimura    507:        Res=getgcd(Res,Rea)$
1.4       kimura    508:
                    509:        if(nonposdegchk(Res)){
                    510:
1.26      kimura    511:                TMP1=makeret(Res,Vars,0)$
1.4       kimura    512:
1.26      kimura    513:                if(checktd(PolyList,Vars,TMP1[1])==1){
1.18      kimura    514:
1.26      kimura    515:                        if(FLAG==0){
1.24      kimura    516:
1.26      kimura    517:                                if(TMP1[0]==0)
                    518:                                        RET=append(RET,wsort(TMP1[1],Vars,TMP1[1],0))$
                    519:                                else{
1.24      kimura    520:
1.26      kimura    521:                                        TMP=vtol(TMP1[1])$
1.27      kimura    522:                                        RET0=[]$
1.26      kimura    523:                                        if((TMP0=fixedpoint(TMP,0))!=[]){
1.24      kimura    524:
1.26      kimura    525:                                                for(I=0;I<length(TMP0);I++)
                    526:                                                        TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.27      kimura    527:                                                RET0=value2(Vars,TMP)$
                    528:                                                if(RET0!=[])
1.32    ! kimura    529:                                                        RET0=wsort(RET0,Vars,RET0,-1)$
1.27      kimura    530:                                        }
                    531:
                    532:                                        TMP=vtol(TMP1[1])$
                    533:                                        if(RET0==[] && (TMP0=fixedpoint(TMP,1))!=[]){
1.24      kimura    534:
1.26      kimura    535:                                                for(I=0;I<length(TMP0);I++)
                    536:                                                        TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.27      kimura    537:                                                RET0=value2(Vars,TMP)$
1.26      kimura    538:
1.27      kimura    539:                                                if(RET0!=[])
1.32    ! kimura    540:                                                        RET0=wsort(RET0,Vars,RET0,-1)$
1.26      kimura    541:                                        }
1.27      kimura    542:                                        RET=append(RET,RET0)$
1.26      kimura    543:                                }
1.10      kimura    544:                        }
1.26      kimura    545:                        else if(FLAG==1)
                    546:                                RET=append(RET,[[0,Vars,vtol(TMP1[1])]])$
1.4       kimura    547:                }
                    548:        }
                    549:
1.26      kimura    550:        return RET$
1.4       kimura    551: }
                    552:
1.32    ! kimura    553: def unitweight2(NormMat0,ExpMat,Vars,FLAG,ID){
1.8       kimura    554:
                    555:        RET=[]$
1.4       kimura    556:
                    557:        ExpMatRowNum=size(ExpMat)[0]$
                    558:        ExpMatColNum=size(ExpMat[0])[0]$
1.32    ! kimura    559:        ExtMatColNum=ExpMatColNum+1$
        !           560:
        !           561:        ExtVars=append(Vars,[uc()])$
1.4       kimura    562:
1.8       kimura    563:        if(NormMat==0){
1.32    ! kimura    564:
        !           565:                NormMat0=newvect(ExtMatColNum)$
        !           566:                for(I=0;I<ExtMatColNum;I++)
        !           567:                        NormMat0[I]=newvect(ExtMatColNum)$
1.8       kimura    568:
                    569:                for(I=0;I<ExpMatColNum;I++)
                    570:                        for(J=I;J<ExpMatColNum;J++)
                    571:                                for(K=0;K<ExpMatRowNum;K++)
1.32    ! kimura    572:                                        NormMat0[I][J]+=
        !           573:                                                ExpMat[K][I]*
        !           574:                                                ExpMat[K][J]$
1.8       kimura    575:        }
1.4       kimura    576:
                    577:        for(I=0;I<ExpMatColNum;I++)
1.32    ! kimura    578:                for(K=0;K<ExpMatRowNum;K++)
        !           579:                        NormMat0[I][ExpMatColNum]-=ExpMat[K][I]$
1.4       kimura    580:
1.32    ! kimura    581:        NormMat0[ExpMatColNum][ExpMatColNum]=ExpMatRowNum$
1.8       kimura    582:
1.32    ! kimura    583:        WorkMat=newvect(ExtMatColNum)$
        !           584:        for(I=0;I<ExtMatColNum;I++)
        !           585:                WorkMat[I]=newvect(ExtMatColNum)$
1.4       kimura    586:
1.32    ! kimura    587:        if(jacobi(ExtMatColNum,NormMat0,WorkMat)){
1.4       kimura    588:
1.32    ! kimura    589:                Res=newvect(ExtMatColNum)$
        !           590:                for(I=0;I<ExtMatColNum;I++){
        !           591:                        Res[I]=newvect(2)$
        !           592:                        Res[I][0]=ExtVars[I]$
        !           593:                        Res[I][1]=WorkMat[ExtMatColNum-1][I]$
        !           594:                }
1.8       kimura    595:
1.32    ! kimura    596:                if(nonposdegchk(Res)){
1.8       kimura    597:
1.32    ! kimura    598:                        TMP1=makeret(Res,Vars,1)$
1.8       kimura    599:
1.32    ! kimura    600:                        if(FLAG==0){
1.26      kimura    601:                                TMP=roundret(TMP1[1])$
1.18      kimura    602:
1.26      kimura    603:                                RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]*1.0),ID))$
1.8       kimura    604:
1.26      kimura    605:                                if(TMP!=[])
                    606:                                        RET=append(RET,wsort(TMP1[1],Vars,TMP,ID+1))$
                    607:                        }
1.32    ! kimura    608:                        else if(FLAG==1)
        !           609:                                RET=append(RET,[[ID,Vars,vtol(TMP1[1])]])$
        !           610:                }
        !           611:        }
1.27      kimura    612:
1.26      kimura    613:        return RET$
1.8       kimura    614: }
                    615:
1.32    ! kimura    616: def unitweight1(ExpMat,Vars,PolyListNum,OneMat,FLAG){
1.8       kimura    617:
                    618:        RET=[]$
                    619:
                    620:        ExpMatRowNum=size(ExpMat)[0]$
                    621:        ExpMatColNum=size(ExpMat[0])[0]$
                    622:        ExtMatColNum=ExpMatColNum+PolyListNum$
                    623:
                    624:        ExtVars=reverse(Vars)$
                    625:        for(I=0;I<PolyListNum;I++)
                    626:                ExtVars=cons(uc(),ExtVars)$
                    627:
                    628:        ExtVars=reverse(ExtVars)$
                    629:
1.32    ! kimura    630:        NormMat0=newvect(ExpMatColNum+1)$
1.18      kimura    631:        for(I=0;I<ExpMatColNum;I++)
1.32    ! kimura    632:                NormMat0[I]=newvect(ExpMatColNum+1)$
1.8       kimura    633:
                    634:        for(I=0;I<ExpMatColNum;I++)
                    635:                for(J=I;J<ExpMatColNum;J++)
                    636:                        for(K=0;K<ExpMatRowNum;K++)
1.18      kimura    637:                                NormMat0[I][J]+=
1.8       kimura    638:                                        ExpMat[K][I]*
                    639:                                        ExpMat[K][J]$
                    640:
1.18      kimura    641:        NormMat1=newvect(ExtMatColNum)$
                    642:        for(I=0;I<ExtMatColNum;I++)
                    643:                NormMat1[I]=newvect(ExtMatColNum)$
                    644:
                    645:        WorkMat=newvect(ExtMatColNum)$
                    646:        for(I=0;I<ExtMatColNum;I++)
                    647:                WorkMat[I]=newvect(ExtMatColNum)$
                    648:
                    649:        for(I=0;I<ExpMatColNum;I++)
                    650:                for(J=I;J<ExpMatColNum;J++)
                    651:                        NormMat1[I][J]=NormMat0[I][J]$
                    652:
1.8       kimura    653:        for(I=0;I<ExpMatColNum;I++)
                    654:                for(J=0;J<PolyListNum;J++)
                    655:                        for(K=OneMat[J];K<OneMat[J+1];K++)
1.26      kimura    656:                                NormMat1[I][J+ExpMatColNum]-=ExpMat[K][I]$
1.8       kimura    657:
                    658:        for(I=0;I<PolyListNum;I++)
1.18      kimura    659:                NormMat1[I+ExpMatColNum][I+ExpMatColNum]=OneMat[I+1]-OneMat[I]$
                    660:
                    661:        if(jacobi(ExtMatColNum,NormMat1,WorkMat)){
1.8       kimura    662:
1.30      kimura    663:                Res=newvect(ExtMatColNum)$
                    664:                for(I=0;I<ExtMatColNum;I++){
1.18      kimura    665:                        Res[I]=newvect(2)$
1.30      kimura    666:                        Res[I][0]=ExtVars[I]$
1.18      kimura    667:                        Res[I][1]=WorkMat[ExtMatColNum-1][I]$
                    668:                }
1.8       kimura    669:
                    670:                if(nonposdegchk(Res)){
                    671:
1.10      kimura    672:                        TMP1=makeret(Res,Vars,1)$
1.16      kimura    673:
1.26      kimura    674:                        if(FLAG==0){
                    675:                                TMP=roundret(TMP1[1])$
1.8       kimura    676:
1.26      kimura    677:                                RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]*1.0),1))$
1.18      kimura    678:
1.26      kimura    679:                                if(TMP!=[])
                    680:                                        RET=append(RET,wsort(TMP1[1],Vars,TMP,2))$
                    681:                        }
1.27      kimura    682:                        else if(FLAG==1)
1.26      kimura    683:                                RET=append(RET,[[1,Vars,vtol(TMP1[1])]])$
1.4       kimura    684:                }
1.18      kimura    685:        }
                    686:
                    687:        return [NormMat0,RET]$
1.6       kimura    688: }
                    689:
1.18      kimura    690: def weight(PolyList,Vars,FLAG){
1.6       kimura    691:
                    692:        Vars0=vars(PolyList)$
                    693:        Vars1=[]$
                    694:        for(I=0;I<length(Vars);I++)
                    695:                if(member(Vars[I],Vars0))
                    696:                        Vars1=cons(Vars[I],Vars1)$
                    697:
                    698:        Vars=reverse(Vars1)$
                    699:
                    700:        RET=[]$
                    701:
1.18      kimura    702:        TMP=qcheck(PolyList,Vars,FLAG)$
1.6       kimura    703:
1.8       kimura    704:        if(TMP!=[]){
                    705:                RET=append(RET,TMP)$
1.18      kimura    706:                return RET$
1.6       kimura    707:        }
1.4       kimura    708:
1.6       kimura    709:        dp_ord(2)$
1.4       kimura    710:
1.6       kimura    711:        PolyListNum=length(PolyList)$
1.4       kimura    712:
1.18      kimura    713:        OneMat=newvect(PolyListNum+1,[0])$
                    714:        ExpMat=[]$
                    715:        for(I=0;I<PolyListNum;I++){
                    716:                for(Poly=dp_ptod(PolyList[I],Vars);
                    717:                        Poly!=0;Poly=dp_rest(Poly)){
                    718:                        ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
1.8       kimura    719:                }
1.18      kimura    720:                OneMat[I+1]=length(ExpMat)$
                    721:        }
1.4       kimura    722:
1.18      kimura    723:        ExpMat=reverse(ExpMat)$
                    724:        ExpMat=newvect(length(ExpMat),ExpMat)$
1.4       kimura    725:
1.32    ! kimura    726:        TMP=unitweight1(ExpMat,Vars,PolyListNum,OneMat,FLAG)$
1.8       kimura    727:
1.18      kimura    728:        RET=append(RET,TMP[1])$
1.4       kimura    729:
1.32    ! kimura    730:        TMP=unitweight2(TMP[0],ExpMat,Vars,FLAG,3)$
1.20      kimura    731:
1.32    ! kimura    732:        RET=append(RET,TMP)$
1.4       kimura    733:
1.6       kimura    734:        ExpMat=qsort(ExpMat,junban)$
1.8       kimura    735:
1.6       kimura    736:        ExpMat2=[]$
                    737:        for(I=0;I<size(ExpMat)[0];I++)
                    738:                if(car(ExpMat2)!=ExpMat[I])
                    739:                        ExpMat2=cons(ExpMat[I],ExpMat2)$
1.4       kimura    740:
1.6       kimura    741:        if(size(ExpMat)[0]!=length(ExpMat2)){
                    742:                ExpMat=newvect(length(ExpMat2),ExpMat2)$
1.32    ! kimura    743:                RET=append(RET,unitweight2(0,ExpMat,Vars,FLAG,5))$
1.20      kimura    744:        }
                    745:        else{
                    746:                TMP0=map(ltov,TMP0)$
                    747:
                    748:                for(I=0;I<length(TMP0);I++)
1.24      kimura    749:                        if(TMP0[I][0]==3)
                    750:                                TMP0[I][0]=5$
                    751:                        else if(TMP0[I][0]==4)
                    752:                                TMP0[I][0]=6$
1.20      kimura    753:
                    754:                TMP0=map(vtol,TMP0)$
                    755:
                    756:                RET=append(RET,TMP0)$
1.4       kimura    757:        }
                    758:
1.18      kimura    759:        return RET$
1.4       kimura    760: }
                    761:
                    762: end$

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