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

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

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:
1.37      kimura    309:        for(I=0;I<ResNum;I++){
                    310:        if(member(Res[I][0],Vars)){
                    311:                ResVec[I]=Res[I][1]$
1.16      kimura    312:
1.37      kimura    313:                        if(FLAG){
1.29      kimura    314:                                if(type(ResVec[I])==1){
1.40    ! kimura    315:                                if(M==0)
1.37      kimura    316:                                M=ResVec[I]$
                    317:                        else
                    318:                                if(ResVec[I]<M)
                    319:                                M=ResVec[I]$
1.29      kimura    320:                                }
                    321:                                else
                    322:                                        M=-1$
                    323:                        }
                    324:                }
1.40    ! kimura    325:        }
1.29      kimura    326:
                    327:        if(M!=-1)
1.40    ! kimura    328:                ResVec=map(red,ResVec/M)$
1.4       kimura    329:
1.17      kimura    330:        RET=newvect(VarsNum,Vars)$
1.4       kimura    331:
1.8       kimura    332:        for(I=0;I<ResNum;I++){
                    333:                for(J=0;J<VarsNum;J++)
                    334:                        if(Vars[J]==Res[I][0])
                    335:                                break$
1.4       kimura    336:
1.8       kimura    337:                if(J<VarsNum)
                    338:                        RET[J]=ResVec[I]$
                    339:        }
                    340:
1.4       kimura    341:        for(I=0;I<VarsNum;I++)
1.8       kimura    342:                if(type(RET[I])!=1)
1.40    ! kimura    343:                        return [1,vtol(RET)]$
1.4       kimura    344:
1.40    ! kimura    345:        return [0,vtol(RET)]$
1.4       kimura    346: }
                    347:
                    348: def roundret(V){
                    349:
1.40    ! kimura    350:        VN=length(V)$
1.4       kimura    351:
1.32      kimura    352:        K=1$
1.8       kimura    353:        RET0=V$
1.32      kimura    354:        RET1=map(drint,RET0)$
                    355:        S=0$
                    356:        for(J=0;J<VN;J++)
                    357:                S+=(RET0[J]-RET1[J])^2$
                    358:
                    359:        for(I=2;I<10;I++){
1.40    ! kimura    360:                RET0=I*ltov(V)$
1.32      kimura    361:                RET1=map(drint,RET0)$
                    362:
                    363:                T=0$
                    364:                for(J=0;J<VN;J++)
                    365:                        T+=(RET0[J]-RET1[J])^2$
                    366:
                    367:                if(T<S){
                    368:                        K=I$
                    369:                        S=T$
                    370:                }
1.4       kimura    371:        }
                    372:
1.40    ! kimura    373:        return map(drint,vtol(K*ltov(V)))$
1.4       kimura    374: }
                    375:
                    376: def chkou(L,ExpMat,CHAGORD){
                    377:
1.8       kimura    378:        for(P=1,I=0;I<L;I++){
1.4       kimura    379:                Q=ExpMat[L][CHAGORD[I]]$
                    380:                for(J=0;J<size(ExpMat[0])[0];J++){
                    381:                        ExpMat[L][CHAGORD[J]]=red((ExpMat[I][CHAGORD[I]]
                    382:                                *ExpMat[L][CHAGORD[J]]-
                    383:                                        Q*ExpMat[I][CHAGORD[J]])/P)$
                    384:                }
                    385:
                    386:                P=ExpMat[I][CHAGORD[I]]$
                    387:        }
                    388:
                    389:        for(J=0;J<size(ExpMat[0])[0];J++)
                    390:                if(ExpMat[L][CHAGORD[J]]!=0)
                    391:                        break$
                    392:
                    393:        if(J==size(ExpMat[0])[0])
                    394:                return L$
                    395:        else{
                    396:                TMP=CHAGORD[L]$
                    397:                CHAGORD[L]=CHAGORD[J]$
                    398:                CHAGORD[J]=TMP$
                    399:                return (L+1)$
                    400:        }
                    401: }
                    402:
1.8       kimura    403: def qcheckmain(PolyList,Vars){
1.4       kimura    404:
                    405:        RET=[]$
                    406:        PolyListNum=length(PolyList)$
                    407:        VarsNum=length(Vars)$
                    408:
                    409:        ExpMat=newvect(VarsNum)$
                    410:        CHAGORD=newvect(VarsNum)$
                    411:        for(I=0;I<VarsNum;I++)
                    412:                CHAGORD[I]=I$
                    413:
                    414:        L=0$
                    415:        for(I=0;I<PolyListNum;I++){
                    416:                Poly=dp_ptod(PolyList[I],Vars)$
                    417:                BASE0=dp_etov(dp_ht(Poly))$
                    418:                Poly=dp_rest(Poly)$
1.39      kimura    419:                for(;L!=VarsNum && Poly!=0;Poly=dp_rest(Poly)){
1.4       kimura    420:                        ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
                    421:                        L=chkou(L,ExpMat,CHAGORD)$
1.8       kimura    422:                        if(L==VarsNum-1)
                    423:                                return [L,CHAGORD,ExpMat]$
1.4       kimura    424:                }
                    425:        }
                    426:
1.8       kimura    427:        return [L,CHAGORD,ExpMat]$
1.4       kimura    428: }
                    429:
                    430: def inner(A,B){
                    431:
                    432:        SUM=0$
                    433:        for(I=0;I<size(A)[0];I++)
                    434:                SUM+=A[I]*B[I]$
                    435:
                    436:        return SUM$
                    437: }
                    438:
                    439: def checktd(PolyList,Vars,ResVars){
                    440:
                    441:        PolyListNum=length(PolyList)$
                    442:        VarsNum=length(Vars)$
                    443:
                    444:        L=0$
                    445:        for(I=0;I<PolyListNum;I++){
                    446:                Poly=dp_ptod(PolyList[I],Vars)$
                    447:                J0=inner(dp_etov(dp_ht(Poly)),ResVars)$
                    448:                Poly=dp_rest(Poly)$
                    449:                for(;Poly!=0;Poly=dp_rest(Poly))
                    450:                        if(J0!=inner(dp_etov(dp_ht(Poly)),ResVars))
                    451:                                return 0$
                    452:        }
                    453:
                    454:        return 1$
                    455: }
                    456:
1.38      kimura    457: def value2(Vars,Ans,Ba,FLAG){
1.24      kimura    458:
                    459:        N=length(Vars)$
                    460:        Res=newvect(N)$
                    461:        for(I=0;I<N;I++){
                    462:                Res[I]=newvect(2)$
                    463:                Res[I][0]=Vars[I]$
1.33      kimura    464:                Res[I][1]=Ba*Ans[I]$
1.24      kimura    465:        }
1.31      kimura    466:        Res=map(vtol,Res)$
                    467:        Res=vtol(Res)$
1.24      kimura    468:
                    469:        Res=getgcd(Res,Vars)$
                    470:
1.40    ! kimura    471:        if(nonposdegchk(Res))
        !           472:                return makeret(Res,Vars,FLAG)$
1.24      kimura    473:        else
                    474:                return []$
                    475: }
                    476:
1.10      kimura    477: def qcheck(PolyList,Vars,FLAG){
1.4       kimura    478:
1.10      kimura    479:        RET=[]$
1.8       kimura    480:        Res=qcheckmain(PolyList,Vars)$
1.4       kimura    481:        VarsNum=length(Vars)$
                    482:
                    483:        IndNum=Res[0]$
                    484:        CHAGORD=Res[1]$
                    485:        ExpMat=Res[2]$
                    486:
                    487:        SolveList=[]$
                    488:        for(I=0;I<IndNum;I++){
                    489:                TMP=0$
                    490:                for(J=0;J<VarsNum;J++)
                    491:                        TMP+=ExpMat[I][CHAGORD[J]]*Vars[CHAGORD[J]]$
                    492:
                    493:                SolveList=cons(TMP,SolveList)$
                    494:        }
                    495:
1.8       kimura    496:        Rea=vars(SolveList)$
                    497:
1.4       kimura    498:        VarsList=[]$
                    499:        for(I=0;I<VarsNum;I++)
1.31      kimura    500:                if(member(TMP0=Vars[CHAGORD[I]],Rea))
1.8       kimura    501:                        VarsList=cons(Vars[CHAGORD[I]],VarsList)$
1.4       kimura    502:
                    503:        Res=solve(reverse(SolveList),reverse(VarsList))$
1.8       kimura    504:        Res=getgcd(Res,Rea)$
1.4       kimura    505:
                    506:        if(nonposdegchk(Res)){
                    507:
1.26      kimura    508:                TMP1=makeret(Res,Vars,0)$
1.4       kimura    509:
1.26      kimura    510:                if(checktd(PolyList,Vars,TMP1[1])==1){
1.18      kimura    511:
1.26      kimura    512:                        if(FLAG==0){
1.24      kimura    513:
1.40    ! kimura    514:                                if(TMP1[0]==0){
1.26      kimura    515:                                        RET=append(RET,wsort(TMP1[1],Vars,TMP1[1],0))$
1.40    ! kimura    516:                                }
1.26      kimura    517:                                else{
1.24      kimura    518:
1.40    ! kimura    519:                                        TMP=TMP1[1]$
        !           520:                                        RET1=[]$
1.26      kimura    521:                                        if((TMP0=fixedpoint(TMP,0))!=[]){
1.40    ! kimura    522:
1.26      kimura    523:                                                for(I=0;I<length(TMP0);I++)
                    524:                                                        TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.40    ! kimura    525:
1.38      kimura    526:                                                RET0=value2(Vars,TMP,1,0)$
1.37      kimura    527:
1.27      kimura    528:                                                if(RET0!=[])
1.40    ! kimura    529:                                                        RET1=wsort(RET0[1],Vars,RET0[1],-1)$
1.27      kimura    530:                                        }
1.40    ! kimura    531:
        !           532:                                        TMP=TMP1[1]$
        !           533:                                        if(RET1==[] && (TMP0=fixedpoint(TMP,1))!=[]){
        !           534:
1.26      kimura    535:                                                for(I=0;I<length(TMP0);I++)
                    536:                                                        TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.40    ! kimura    537:
1.38      kimura    538:                                                RET0=value2(Vars,TMP,-1,0)$
1.26      kimura    539:
1.27      kimura    540:                                                if(RET0!=[])
1.40    ! kimura    541:                                                        RET1=wsort(RET0[1],Vars,RET0[1],-1)$
1.26      kimura    542:                                        }
1.40    ! kimura    543:
        !           544:                                        if(RET1!=[])
        !           545:                                                RET=append(RET,RET1)$
        !           546:
1.26      kimura    547:                                }
1.40    ! kimura    548:
1.10      kimura    549:                        }
1.26      kimura    550:                        else if(FLAG==1)
1.40    ! kimura    551:                                RET=append(RET,[[0,Vars,TMP1[1]]])$
1.4       kimura    552:                }
                    553:        }
                    554:
1.26      kimura    555:        return RET$
1.4       kimura    556: }
                    557:
1.32      kimura    558: def unitweight1(ExpMat,Vars,PolyListNum,OneMat,FLAG){
1.8       kimura    559:
                    560:        RET=[]$
                    561:
                    562:        ExpMatRowNum=size(ExpMat)[0]$
                    563:        ExpMatColNum=size(ExpMat[0])[0]$
                    564:        ExtMatColNum=ExpMatColNum+PolyListNum$
                    565:
                    566:        ExtVars=reverse(Vars)$
                    567:        for(I=0;I<PolyListNum;I++)
                    568:                ExtVars=cons(uc(),ExtVars)$
                    569:
                    570:        ExtVars=reverse(ExtVars)$
                    571:
1.32      kimura    572:        NormMat0=newvect(ExpMatColNum+1)$
1.18      kimura    573:        for(I=0;I<ExpMatColNum;I++)
1.32      kimura    574:                NormMat0[I]=newvect(ExpMatColNum+1)$
1.8       kimura    575:
                    576:        for(I=0;I<ExpMatColNum;I++)
                    577:                for(J=I;J<ExpMatColNum;J++)
                    578:                        for(K=0;K<ExpMatRowNum;K++)
1.18      kimura    579:                                NormMat0[I][J]+=
1.8       kimura    580:                                        ExpMat[K][I]*
                    581:                                        ExpMat[K][J]$
                    582:
1.18      kimura    583:        NormMat1=newvect(ExtMatColNum)$
                    584:        for(I=0;I<ExtMatColNum;I++)
                    585:                NormMat1[I]=newvect(ExtMatColNum)$
                    586:
                    587:        WorkMat=newvect(ExtMatColNum)$
                    588:        for(I=0;I<ExtMatColNum;I++)
                    589:                WorkMat[I]=newvect(ExtMatColNum)$
                    590:
                    591:        for(I=0;I<ExpMatColNum;I++)
                    592:                for(J=I;J<ExpMatColNum;J++)
                    593:                        NormMat1[I][J]=NormMat0[I][J]$
                    594:
1.8       kimura    595:        for(I=0;I<ExpMatColNum;I++)
                    596:                for(J=0;J<PolyListNum;J++)
                    597:                        for(K=OneMat[J];K<OneMat[J+1];K++)
1.26      kimura    598:                                NormMat1[I][J+ExpMatColNum]-=ExpMat[K][I]$
1.8       kimura    599:
                    600:        for(I=0;I<PolyListNum;I++)
1.18      kimura    601:                NormMat1[I+ExpMatColNum][I+ExpMatColNum]=OneMat[I+1]-OneMat[I]$
                    602:
                    603:        if(jacobi(ExtMatColNum,NormMat1,WorkMat)){
1.8       kimura    604:
1.30      kimura    605:                Res=newvect(ExtMatColNum)$
                    606:                for(I=0;I<ExtMatColNum;I++){
1.18      kimura    607:                        Res[I]=newvect(2)$
1.30      kimura    608:                        Res[I][0]=ExtVars[I]$
1.18      kimura    609:                        Res[I][1]=WorkMat[ExtMatColNum-1][I]$
                    610:                }
1.8       kimura    611:
                    612:                if(nonposdegchk(Res)){
                    613:
1.10      kimura    614:                        TMP1=makeret(Res,Vars,1)$
1.16      kimura    615:
1.26      kimura    616:                        if(FLAG==0){
1.38      kimura    617:                                RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]),1))$
1.18      kimura    618:
1.40    ! kimura    619:                                TMP=roundret(TMP1[1])$
1.26      kimura    620:                                if(TMP!=[])
                    621:                                        RET=append(RET,wsort(TMP1[1],Vars,TMP,2))$
1.40    ! kimura    622:
1.26      kimura    623:                        }
1.27      kimura    624:                        else if(FLAG==1)
1.40    ! kimura    625:                                RET=append(RET,[[1,Vars,TMP1[1]]])$
1.4       kimura    626:                }
1.18      kimura    627:        }
                    628:
                    629:        return [NormMat0,RET]$
1.6       kimura    630: }
                    631:
1.34      kimura    632: def leastsq(NormMat,ExpMat,Vars,FLAG,ID){
                    633:
                    634:        RET=[]$
                    635:
                    636:        ExpMatRowNum=size(ExpMat)[0]$
                    637:        ExpMatColNum=size(ExpMat[0])[0]$
                    638:
                    639:        if(NormMat==0){
                    640:                NormMat=newmat(ExpMatColNum,ExpMatColNum)$
                    641:
                    642:                for(I=0;I<ExpMatColNum;I++)
                    643:                        for(J=I;J<ExpMatColNum;J++)
                    644:                                for(K=0;K<ExpMatRowNum;K++)
                    645:                                        NormMat[I][J]+=
                    646:                                                ExpMat[K][I]*ExpMat[K][J]$
                    647:        }
                    648:
                    649:        BVec=newvect(ExpMatColNum)$
                    650:
                    651:        for(I=0;I<ExpMatColNum;I++)
                    652:                for(J=0;J<ExpMatRowNum;J++)
                    653:                        BVec[I]+=ExpMat[J][I]$
                    654:
                    655:        SolveList=[]$
                    656:        for(I=0;I<ExpMatColNum;I++){
                    657:                TMP=0$
                    658:                for(J=0;J<I;J++)
                    659:                        TMP+=NormMat[J][I]*Vars[J]$
                    660:
                    661:                for(J=I;J<ExpMatColNum;J++)
                    662:                        TMP+=NormMat[I][J]*Vars[J]$
                    663:
                    664:                TMP-=BVec[I]$
                    665:                SolveList=cons(TMP,SolveList)$
                    666:        }
                    667:
                    668:        Rea=vars(SolveList)$
                    669:
                    670:        VarsList=[]$
                    671:        for(I=0;I<length(Vars);I++)
                    672:                if(member(Vars[I],Rea))
                    673:                        VarsList=cons(Vars[I],VarsList)$
                    674:
                    675:        Res=solve(SolveList,VarsList)$
                    676:        Res=getgcd(Res,Rea)$
                    677:
                    678:        if(nonposdegchk(Res)){
                    679:
                    680:                TMP1=makeret(Res,Vars,1)$
                    681:
                    682:                if(FLAG==0){
                    683:
                    684:                        if(TMP1[0]==0){
1.40    ! kimura    685:                                RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]),ID))$
1.34      kimura    686:
                    687:                                TMP=roundret(TMP1[1])$
                    688:
                    689:                                if(TMP!=[])
                    690:                                        RET=append(RET,wsort(TMP1[1],Vars,TMP,ID+1))$
                    691:                        }
                    692:                        else{
                    693:
1.40    ! kimura    694:                                TMP=TMP1[1]$
        !           695:                                RET1=[]$
1.34      kimura    696:                                if((TMP0=fixedpoint(TMP,0))!=[]){
                    697:
                    698:                                        for(I=0;I<length(TMP0);I++)
                    699:                                                TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.38      kimura    700:                                        RET0=value2(Vars,TMP,1,1)$
1.37      kimura    701:
1.34      kimura    702:                                        if(RET0!=[])
1.40    ! kimura    703:                                                RET1=wsort(RET0[1],Vars,map(drint,RET0[1]),-ID)$
        !           704:
1.34      kimura    705:                                }
                    706:
1.40    ! kimura    707:                                TMP=TMP1[1]$
        !           708:                                if(RET1==[] && (TMP0=fixedpoint(TMP,1))!=[]){
1.34      kimura    709:
                    710:                                        for(I=0;I<length(TMP0);I++)
                    711:                                                TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.38      kimura    712:                                        RET0=value2(Vars,TMP,-1,1)$
1.34      kimura    713:
                    714:                                        if(RET0!=[])
1.40    ! kimura    715:                                                RET1=wsort(RET0[1],Vars,map(drint,RET0[1]),-ID)$
        !           716:                                }
        !           717:
        !           718:                                if(RET1!=[]){
        !           719:                                        RET=append(RET,RET1)$
        !           720:                                        TMP=roundret(RET0[1])$
        !           721:                                        if(TMP!=[])
        !           722:                                                RET=append(RET,wsort(RET0[1],Vars,TMP,-(ID+1)))$
1.34      kimura    723:                                }
                    724:
                    725:                        }
                    726:
                    727:                }
                    728:                else if(FLAG==1)
1.40    ! kimura    729:                        RET=append(RET,[[ID,Vars,TMP1[1]]])$
1.34      kimura    730:        }
                    731:
1.36      kimura    732:        return [NormMat0,RET]$
1.34      kimura    733: }
                    734:
1.18      kimura    735: def weight(PolyList,Vars,FLAG){
1.6       kimura    736:
                    737:        Vars0=vars(PolyList)$
                    738:        Vars1=[]$
                    739:        for(I=0;I<length(Vars);I++)
                    740:                if(member(Vars[I],Vars0))
                    741:                        Vars1=cons(Vars[I],Vars1)$
                    742:
                    743:        Vars=reverse(Vars1)$
                    744:
                    745:        RET=[]$
                    746:
1.18      kimura    747:        TMP=qcheck(PolyList,Vars,FLAG)$
1.6       kimura    748:
1.8       kimura    749:        if(TMP!=[]){
                    750:                RET=append(RET,TMP)$
1.18      kimura    751:                return RET$
1.6       kimura    752:        }
1.4       kimura    753:
1.6       kimura    754:        dp_ord(2)$
1.4       kimura    755:
1.6       kimura    756:        PolyListNum=length(PolyList)$
1.4       kimura    757:
1.18      kimura    758:        OneMat=newvect(PolyListNum+1,[0])$
                    759:        ExpMat=[]$
                    760:        for(I=0;I<PolyListNum;I++){
                    761:                for(Poly=dp_ptod(PolyList[I],Vars);
                    762:                        Poly!=0;Poly=dp_rest(Poly)){
                    763:                        ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
1.8       kimura    764:                }
1.18      kimura    765:                OneMat[I+1]=length(ExpMat)$
                    766:        }
1.4       kimura    767:
1.18      kimura    768:        ExpMat=reverse(ExpMat)$
                    769:        ExpMat=newvect(length(ExpMat),ExpMat)$
1.4       kimura    770:
1.32      kimura    771:        TMP=unitweight1(ExpMat,Vars,PolyListNum,OneMat,FLAG)$
1.34      kimura    772:        if(TMP[1]!=[])
                    773:                RET=append(RET,TMP[1])$
1.4       kimura    774:
1.37      kimura    775:        TMP=leastsq(0,ExpMat,Vars,FLAG,3)$
1.34      kimura    776:        if(TMP[1]!=[])
                    777:                RET=append(RET,TMP[1])$
1.4       kimura    778:
1.6       kimura    779:        ExpMat=qsort(ExpMat,junban)$
1.8       kimura    780:
1.6       kimura    781:        ExpMat2=[]$
                    782:        for(I=0;I<size(ExpMat)[0];I++)
                    783:                if(car(ExpMat2)!=ExpMat[I])
                    784:                        ExpMat2=cons(ExpMat[I],ExpMat2)$
1.4       kimura    785:
1.6       kimura    786:        if(size(ExpMat)[0]!=length(ExpMat2)){
                    787:                ExpMat=newvect(length(ExpMat2),ExpMat2)$
1.37      kimura    788:                TMP=leastsq(0,ExpMat,Vars,FLAG,5)$
1.34      kimura    789:                if(TMP[1]!=[])
                    790:                        RET=append(RET,TMP[1])$
1.20      kimura    791:        }
                    792:        else{
1.35      kimura    793:                TMP=map(ltov,TMP[1])$
1.20      kimura    794:
1.35      kimura    795:                for(I=0;I<length(TMP);I++){
                    796:                        if(TMP[I][0]==3)
                    797:                                TMP[I][0]=5$
                    798:                        else if(TMP[I][0]==4)
                    799:                                TMP[I][0]=6$
1.40    ! kimura    800:                        else if(TMP[I][0]==-3)
        !           801:                                TMP[I][0]=-5$
        !           802:                        else if(TMP[I][0]==-4)
        !           803:                                TMP[I][0]=-6$
        !           804:
1.35      kimura    805:                }
1.20      kimura    806:
1.35      kimura    807:                TMP=map(vtol,TMP)$
1.20      kimura    808:
1.35      kimura    809:                RET=append(RET,TMP)$
1.4       kimura    810:        }
                    811:
1.18      kimura    812:        return RET$
1.4       kimura    813: }
                    814:
                    815: end$
1.34      kimura    816:
                    817:

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