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

Annotation of OpenXM_contrib2/asir2018/lib/weight, Revision 1.1

1.1     ! noro        1: load("solve")$
        !             2: load("gr")$
        !             3:
        !             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:
        !           114: def interval2value(A,Vars){
        !           115:
        !           116:        B=atl(A)$
        !           117:
        !           118:        if(length(B)>2){
        !           119:                print("bug")$
        !           120:                return []$
        !           121:        }
        !           122:        else if(length(B)==0){
        !           123:                if(fop(A)==0)
        !           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]$
        !           138:                else if(fop(B[0])==8)
        !           139:                        return [Vars,E[0][1]]$
        !           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:
        !           201: def junban(A,B){
        !           202:        return (A<B ? 1:(A>B ? -1:0))$
        !           203: }
        !           204:
        !           205: def bsort(A){
        !           206:
        !           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:
        !           222: def wsort(A,B,C,ID){
        !           223:
        !           224:        D=newvect(length(B))$
        !           225:        for(I=0;I<length(B);I++)
        !           226:                D[I]=[A[I],B[I],C[I]]$
        !           227:
        !           228:        D=bsort(D)$
        !           229:
        !           230:        for(E=[],I=0;I<length(B);I++)
        !           231:                E=cons(D[I][1],E)$
        !           232:        E=reverse(E)$
        !           233:
        !           234:        for(F=[],I=0;I<length(B);I++)
        !           235:                F=cons(D[I][2],F)$
        !           236:        F=reverse(F)$
        !           237:
        !           238:        return [[ID,E,F]]$
        !           239: }
        !           240:
        !           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:
        !           250: def getgcd(A,B){
        !           251:
        !           252:        Anum=length(A)$
        !           253:
        !           254:        TMP=[]$
        !           255:        for(I=0;I<length(B);I++){
        !           256:
        !           257:                for(J=0;J<Anum;J++)
        !           258:                        if(B[I]==A[J][0])
        !           259:                                break;
        !           260:
        !           261:                if(J==Anum)
        !           262:                        TMP=cons([B[I],B[I]],TMP)$
        !           263:        }
        !           264:        A=append(A,TMP)$
        !           265:
        !           266:        Anum=length(A)$
        !           267:        A=map(ltov,A)$
        !           268:
        !           269:        for(D=0,I=0;I<Anum;I++)
        !           270:                D=gcd(D,A[I][1])$
        !           271:
        !           272:        if(D!=0){
        !           273:                for(I=0;I<Anum;I++)
        !           274:                        A[I][1]=red(A[I][1]/D)$
        !           275:        }
        !           276:
        !           277:        for(L=1,D=0,I=0;I<Anum;I++){
        !           278:                if(type(TMP=dn(A[I][1]))==1)
        !           279:                        L=ilcm(L,TMP)$
        !           280:
        !           281:                if(type(TMP=nm(A[I][1]))==1)
        !           282:                        D=igcd(D,TMP)$
        !           283:        }
        !           284:
        !           285:        for(I=0;I<Anum;I++)
        !           286:                A[I][1]=A[I][1]*L$
        !           287:
        !           288:        if(D!=0){
        !           289:
        !           290:                for(I=0;I<Anum;I++)
        !           291:                        A[I][1]=A[I][1]/D$
        !           292:        }
        !           293:
        !           294:        return map(vtol,A)$
        !           295: }
        !           296:
        !           297: def makeret(Res,Vars,FLAG){
        !           298:
        !           299:        ResNum=length(Res)$
        !           300:        VarsNum=length(Vars)$
        !           301:
        !           302:        ResVec=newvect(ResNum)$
        !           303:
        !           304:        if(FLAG)
        !           305:                M=0$
        !           306:        else
        !           307:                M=-1$
        !           308:
        !           309:        for(I=0;I<ResNum;I++){
        !           310:        if(member(Res[I][0],Vars)){
        !           311:                ResVec[I]=Res[I][1]$
        !           312:
        !           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:                }
        !           325:        }
        !           326:
        !           327:        if(M!=-1)
        !           328:                ResVec=map(red,ResVec/M)$
        !           329:
        !           330:        RET=newvect(VarsNum,Vars)$
        !           331:
        !           332:        for(I=0;I<ResNum;I++){
        !           333:                for(J=0;J<VarsNum;J++)
        !           334:                        if(Vars[J]==Res[I][0])
        !           335:                                break$
        !           336:
        !           337:                if(J<VarsNum)
        !           338:                        RET[J]=ResVec[I]$
        !           339:        }
        !           340:
        !           341:        for(I=0;I<VarsNum;I++)
        !           342:                if(type(RET[I])!=1)
        !           343:                        return [1,vtol(RET)]$
        !           344:
        !           345:        return [0,vtol(RET)]$
        !           346: }
        !           347:
        !           348: def roundret(V){
        !           349:
        !           350:        VN=length(V)$
        !           351:
        !           352:        K=1$
        !           353:        RET0=V$
        !           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++){
        !           360:                RET0=I*ltov(V)$
        !           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:                }
        !           371:        }
        !           372:
        !           373:        return map(drint,vtol(K*ltov(V)))$
        !           374: }
        !           375:
        !           376: def chkou(L,ExpMat,CHAGORD){
        !           377:
        !           378:        for(P=1,I=0;I<L;I++){
        !           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:
        !           403: def qcheckmain(PolyList,Vars){
        !           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)$
        !           419:                for(;L!=VarsNum && Poly!=0;Poly=dp_rest(Poly)){
        !           420:                        ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
        !           421:                        L=chkou(L,ExpMat,CHAGORD)$
        !           422:                        if(L==VarsNum-1)
        !           423:                                return [L,CHAGORD,ExpMat]$
        !           424:                }
        !           425:        }
        !           426:
        !           427:        return [L,CHAGORD,ExpMat]$
        !           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:
        !           457: def value2(Vars,Ans,Ba,FLAG){
        !           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]$
        !           464:                Res[I][1]=Ba*Ans[I]$
        !           465:        }
        !           466:        Res=map(vtol,Res)$
        !           467:        Res=vtol(Res)$
        !           468:
        !           469:        Res=getgcd(Res,Vars)$
        !           470:
        !           471:        if(nonposdegchk(Res))
        !           472:                return makeret(Res,Vars,FLAG)$
        !           473:        else
        !           474:                return []$
        !           475: }
        !           476:
        !           477: def qcheck(PolyList,Vars,FLAG){
        !           478:
        !           479:        RET=[]$
        !           480:        Res=qcheckmain(PolyList,Vars)$
        !           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:
        !           496:        Rea=vars(SolveList)$
        !           497:
        !           498:        VarsList=[]$
        !           499:        for(I=0;I<VarsNum;I++)
        !           500:                if(member(TMP0=Vars[CHAGORD[I]],Rea))
        !           501:                        VarsList=cons(Vars[CHAGORD[I]],VarsList)$
        !           502:
        !           503:        Res=solve(reverse(SolveList),reverse(VarsList))$
        !           504:        Res=getgcd(Res,Rea)$
        !           505:
        !           506:        if(nonposdegchk(Res)){
        !           507:
        !           508:                TMP1=makeret(Res,Vars,0)$
        !           509:
        !           510:                if(checktd(PolyList,Vars,TMP1[1])==1){
        !           511:
        !           512:                        if(FLAG==0){
        !           513:
        !           514:                                if(TMP1[0]==0){
        !           515:                                        RET=append(RET,wsort(TMP1[1],Vars,TMP1[1],0))$
        !           516:                                }
        !           517:                                else{
        !           518:
        !           519:                                        TMP=TMP1[1]$
        !           520:                                        RET1=[]$
        !           521:                                        if((TMP0=fixedpoint(TMP,0))!=[]){
        !           522:
        !           523:                                                for(I=0;I<length(TMP0);I++)
        !           524:                                                        TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
        !           525:
        !           526:                                                RET0=value2(Vars,TMP,1,0)$
        !           527:
        !           528:                                                if(RET0!=[])
        !           529:                                                        RET1=wsort(RET0[1],Vars,RET0[1],-1)$
        !           530:                                        }
        !           531:
        !           532:                                        TMP=TMP1[1]$
        !           533:                                        if(RET1==[] && (TMP0=fixedpoint(TMP,1))!=[]){
        !           534:
        !           535:                                                for(I=0;I<length(TMP0);I++)
        !           536:                                                        TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
        !           537:
        !           538:                                                RET0=value2(Vars,TMP,-1,0)$
        !           539:
        !           540:                                                if(RET0!=[])
        !           541:                                                        RET1=wsort(RET0[1],Vars,RET0[1],-1)$
        !           542:                                        }
        !           543:
        !           544:                                        if(RET1!=[])
        !           545:                                                RET=append(RET,RET1)$
        !           546:
        !           547:                                }
        !           548:
        !           549:                        }
        !           550:                        else if(FLAG==1)
        !           551:                                RET=append(RET,[[0,Vars,TMP1[1]]])$
        !           552:                }
        !           553:        }
        !           554:
        !           555:        return RET$
        !           556: }
        !           557:
        !           558: def unitweight1(ExpMat,Vars,PolyListNum,OneMat,FLAG){
        !           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:
        !           572:        NormMat0=newvect(ExpMatColNum+1)$
        !           573:        for(I=0;I<ExpMatColNum;I++)
        !           574:                NormMat0[I]=newvect(ExpMatColNum+1)$
        !           575:
        !           576:        for(I=0;I<ExpMatColNum;I++)
        !           577:                for(J=I;J<ExpMatColNum;J++)
        !           578:                        for(K=0;K<ExpMatRowNum;K++)
        !           579:                                NormMat0[I][J]+=
        !           580:                                        ExpMat[K][I]*
        !           581:                                        ExpMat[K][J]$
        !           582:
        !           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:
        !           595:        for(I=0;I<ExpMatColNum;I++)
        !           596:                for(J=0;J<PolyListNum;J++)
        !           597:                        for(K=OneMat[J];K<OneMat[J+1];K++)
        !           598:                                NormMat1[I][J+ExpMatColNum]-=ExpMat[K][I]$
        !           599:
        !           600:        for(I=0;I<PolyListNum;I++)
        !           601:                NormMat1[I+ExpMatColNum][I+ExpMatColNum]=OneMat[I+1]-OneMat[I]$
        !           602:
        !           603:        if(jacobi(ExtMatColNum,NormMat1,WorkMat)){
        !           604:
        !           605:                Res=newvect(ExtMatColNum)$
        !           606:                for(I=0;I<ExtMatColNum;I++){
        !           607:                        Res[I]=newvect(2)$
        !           608:                        Res[I][0]=ExtVars[I]$
        !           609:                        Res[I][1]=WorkMat[ExtMatColNum-1][I]$
        !           610:                }
        !           611:
        !           612:                if(nonposdegchk(Res)){
        !           613:
        !           614:                        TMP1=makeret(Res,Vars,1)$
        !           615:
        !           616:                        if(FLAG==0){
        !           617:                                RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]),1))$
        !           618:
        !           619:                                TMP=roundret(TMP1[1])$
        !           620:                                if(TMP!=[])
        !           621:                                        RET=append(RET,wsort(TMP1[1],Vars,TMP,2))$
        !           622:
        !           623:                        }
        !           624:                        else if(FLAG==1)
        !           625:                                RET=append(RET,[[1,Vars,TMP1[1]]])$
        !           626:                }
        !           627:        }
        !           628:
        !           629:        return [NormMat0,RET]$
        !           630: }
        !           631:
        !           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){
        !           685:                                RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]),ID))$
        !           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:
        !           694:                                TMP=TMP1[1]$
        !           695:                                RET1=[]$
        !           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])$
        !           700:                                        RET0=value2(Vars,TMP,1,1)$
        !           701:
        !           702:                                        if(RET0!=[])
        !           703:                                                RET1=wsort(RET0[1],Vars,map(drint,RET0[1]),-ID)$
        !           704:
        !           705:                                }
        !           706:
        !           707:                                TMP=TMP1[1]$
        !           708:                                if(RET1==[] && (TMP0=fixedpoint(TMP,1))!=[]){
        !           709:
        !           710:                                        for(I=0;I<length(TMP0);I++)
        !           711:                                                TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
        !           712:                                        RET0=value2(Vars,TMP,-1,1)$
        !           713:
        !           714:                                        if(RET0!=[])
        !           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)))$
        !           723:                                }
        !           724:
        !           725:                        }
        !           726:
        !           727:                }
        !           728:                else if(FLAG==1)
        !           729:                        RET=append(RET,[[ID,Vars,TMP1[1]]])$
        !           730:        }
        !           731:
        !           732:        return [NormMat0,RET]$
        !           733: }
        !           734:
        !           735: def weight(PolyList,Vars,FLAG){
        !           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:
        !           747:        TMP=qcheck(PolyList,Vars,FLAG)$
        !           748:
        !           749:        if(TMP!=[]){
        !           750:                RET=append(RET,TMP)$
        !           751:                return RET$
        !           752:        }
        !           753:
        !           754:        dp_ord(2)$
        !           755:
        !           756:        PolyListNum=length(PolyList)$
        !           757:
        !           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)$
        !           764:                }
        !           765:                OneMat[I+1]=length(ExpMat)$
        !           766:        }
        !           767:
        !           768:        ExpMat=reverse(ExpMat)$
        !           769:        ExpMat=newvect(length(ExpMat),ExpMat)$
        !           770:
        !           771:        TMP=unitweight1(ExpMat,Vars,PolyListNum,OneMat,FLAG)$
        !           772:        if(TMP[1]!=[])
        !           773:                RET=append(RET,TMP[1])$
        !           774:
        !           775:        TMP=leastsq(0,ExpMat,Vars,FLAG,3)$
        !           776:        if(TMP[1]!=[])
        !           777:                RET=append(RET,TMP[1])$
        !           778:
        !           779:        ExpMat=qsort(ExpMat,junban)$
        !           780:
        !           781:        ExpMat2=[]$
        !           782:        for(I=0;I<size(ExpMat)[0];I++)
        !           783:                if(car(ExpMat2)!=ExpMat[I])
        !           784:                        ExpMat2=cons(ExpMat[I],ExpMat2)$
        !           785:
        !           786:        if(size(ExpMat)[0]!=length(ExpMat2)){
        !           787:                ExpMat=newvect(length(ExpMat2),ExpMat2)$
        !           788:                TMP=leastsq(0,ExpMat,Vars,FLAG,5)$
        !           789:                if(TMP[1]!=[])
        !           790:                        RET=append(RET,TMP[1])$
        !           791:        }
        !           792:        else{
        !           793:                TMP=map(ltov,TMP[1])$
        !           794:
        !           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$
        !           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:
        !           805:                }
        !           806:
        !           807:                TMP=map(vtol,TMP)$
        !           808:
        !           809:                RET=append(RET,TMP)$
        !           810:        }
        !           811:
        !           812:        return RET$
        !           813: }
        !           814:
        !           815: end$
        !           816:
        !           817:

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