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

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

1.2     ! kimura      1: load("solve")$
        !             2: 
        !             3: def nonposdegchk(Res){
        !             4: 
        !             5:        for(I=0;I<length(Res);I++)
        !             6:                if(Res[I][1]<=0)
        !             7:                        return 0$
        !             8: 
        !             9:        return 1$
        !            10: }
        !            11: 
        !            12: def resvars(Res,Vars){
        !            13: 
        !            14:        ResVars=newvect(length(Vars),Vars)$
        !            15: 
        !            16:        for(I=0;I<length(Res);I++){
        !            17:        
        !            18:                for(J=0;J<size(ResVars)[0];J++)
        !            19:                        if(Res[I][0]==ResVars[J])
        !            20:                                break$
        !            21: 
        !            22:                ResVars[J]=Res[I][1]$
        !            23:        }
        !            24: 
        !            25:        return(ResVars)$
        !            26: }
        !            27: 
        !            28: def makeret(Res,Vars,B){
        !            29: 
        !            30:        VarsNum=length(Vars)$
        !            31: 
        !            32:        ResMat=newvect(VarsNum)$
        !            33:        for(I=0;I<VarsNum;I++)
        !            34:                ResMat[I]=newvect(2)$
        !            35: 
        !            36:        for(I=0;I<VarsNum;I++){
        !            37:                ResMat[I][0]=Vars[I]$
        !            38:                ResMat[I][1]=Vars[I]$
        !            39:        }
        !            40:        
        !            41:        for(F=0,L=1,D=0,M=1,I=0;I<length(Res);I++){
        !            42: 
        !            43:                for(J=0;J<size(ResMat)[0];J++)
        !            44:                        if(Res[I][0]==ResMat[J][0])
        !            45:                                break$
        !            46: 
        !            47:                if(J<VarsNum){
        !            48:                        K=Res[I][1]$
        !            49:                        ResMat[J][1]=K$
        !            50: 
        !            51:                        if(F==0 && type(K)==1){
        !            52:                                if(B==2){
        !            53:                                        L=ilcm(L,dn(K))$
        !            54:                                        D=igcd(D,nm(K))$
        !            55:                                }
        !            56:                                else{
        !            57:                                        if(K<M)
        !            58:                                                M=K$
        !            59:                                }
        !            60:                        }
        !            61:                        else
        !            62:                                F=1$
        !            63: 
        !            64:                }
        !            65:        }
        !            66: 
        !            67:        if(F==0)
        !            68:                if(B==2)
        !            69:                        for(I=0;I<VarsNum;I++)
        !            70:                                ResMat[I][1]=ResMat[I][1]*L/D$
        !            71:                else
        !            72:                        for(I=0;I<VarsNum;I++)
        !            73:                                ResMat[I][1]=ResMat[I][1]/M*1.0$
        !            74: 
        !            75:        for(I=0;I<VarsNum;I++)
        !            76:                for(J=0;J<length(Vars);J++)
        !            77:                        ResMat[I][1]=subst(ResMat[I][1],Vars[J],
        !            78:                                strtov(rtostr(Vars[J])+"_deg"))$
        !            79: 
        !            80:        ResMat=map(vtol,ResMat)$
        !            81:        return(vtol(ResMat))$
        !            82: 
        !            83: }
        !            84: 
        !            85: def afo(A,B){
        !            86: 
        !            87:        for(I=0;I<size(A)[0];I++){
        !            88:                if(A[I]<B[I])
        !            89:                        return 1$
        !            90:                
        !            91:                if(A[I]>B[I])
        !            92:                        return -1$
        !            93:        }
        !            94: 
        !            95:        return 0$
        !            96: }
        !            97: 
        !            98: def weight(PolyList,Vars){
        !            99: 
        !           100:        dp_ord(2)$
        !           101: 
        !           102:        PolyListNum=length(PolyList)$
        !           103: 
        !           104:        ExpMat=[]$
        !           105:        for(I=0;I<PolyListNum;I++)
        !           106:                for(Poly=dp_ptod(PolyList[I],Vars);Poly!=0;Poly=dp_rest(Poly))
        !           107:                        ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
        !           108: 
        !           109:        ExpMat=reverse(ExpMat)$
        !           110:        ExpMat=newvect(length(ExpMat),ExpMat)$
        !           111: 
        !           112:        ExpMatRowNum=size(ExpMat)[0]$
        !           113:        ExpMatColNum=size(ExpMat[0])[0]$
        !           114:        ExtMatColNum=ExpMatColNum+PolyListNum$
        !           115: 
        !           116:        OneMat=newvect(PolyListNum+1,[0])$
        !           117:     for(I=0,SUM=0;I<PolyListNum;I++){
        !           118:                SUM+=nmono(PolyList[I])$
        !           119:                OneMat[I+1]=SUM$
        !           120:        }
        !           121: 
        !           122:        RevOneMat=newvect(ExpMatRowNum)$
        !           123:        for(I=0;I<PolyListNum;I++)
        !           124:                for(J=OneMat[I];J<OneMat[I+1];J++)
        !           125:                        RevOneMat[J]=I$
        !           126: 
        !           127:        NormMat=newmat(ExpMatColNum,ExtMatColNum)$
        !           128: 
        !           129:        for(I=0;I<ExpMatColNum;I++)
        !           130:                for(J=0;J<ExpMatColNum;J++)
        !           131:                        for(K=0;K<ExpMatRowNum;K++)
        !           132:                                NormMat[I][J]+=ExpMat[K][I]*ExpMat[K][J]$
        !           133: 
        !           134:        for(I=0;I<ExpMatColNum;I++)
        !           135:                for(J=0;J<PolyListNum-1;J++)
        !           136:                        for(K=OneMat[J];K<OneMat[J+1];K++)
        !           137:                                NormMat[I][J+ExpMatColNum]-=ExpMat[K][I]$
        !           138: 
        !           139:        for(I=0;I<ExpMatColNum;I++)
        !           140:                for(J=OneMat[PolyListNum-1];J<OneMat[PolyListNum];J++)
        !           141:                        NormMat[I][ExtMatColNum-1]+=ExpMat[J][I]$
        !           142: 
        !           143:        NormMat2=newmat(PolyListNum-1,ExpMatColNum+1)$
        !           144: 
        !           145:        for(I=0;I<PolyListNum-1;I++)
        !           146:                for(J=0;J<ExpMatColNum;J++)
        !           147:                        for(K=OneMat[I];K<OneMat[I+1];K++)
        !           148:                                NormMat2[I][J]-=ExpMat[K][J]$
        !           149: 
        !           150:        for(I=0;I<PolyListNum-1;I++)
        !           151:                NormMat2[I][ExpMatColNum]=OneMat[I+1]-OneMat[I]$
        !           152: 
        !           153:        ExtVars=Vars$
        !           154:        for(I=0;I<PolyListNum-1;I++)
        !           155:                ExtVars=append(ExtVars,[uc()])$
        !           156: 
        !           157:        SolveList=[]$
        !           158:        for(I=0;I<ExpMatColNum;I++){
        !           159:                TMP=0$
        !           160:                for(J=0;J<ExtMatColNum-1;J++)
        !           161:                        TMP+=NormMat[I][J]*ExtVars[J]$
        !           162: 
        !           163:                TMP-=NormMat[I][ExtMatColNum-1]$
        !           164:                SolveList=cons(TMP,SolveList)$
        !           165:        }
        !           166: 
        !           167:        for(I=0;I<PolyListNum-1;I++){
        !           168:                TMP=0$
        !           169:                for(J=0;J<ExpMatColNum;J++)
        !           170:                        TMP+=NormMat2[I][J]*ExtVars[J]$
        !           171: 
        !           172:                TMP+=NormMat2[I][ExpMatColNum]*ExtVars[I+ExpMatColNum]$
        !           173: 
        !           174:                SolveList=cons(TMP,SolveList)$
        !           175:        }
        !           176: 
        !           177:        ReaVars=vars(SolveList)$
        !           178:        Res=solve(SolveList,reverse(ExtVars))$
        !           179: 
        !           180:        RET=[]$
        !           181:        if(nonposdegchk(Res)){
        !           182: 
        !           183:                ResVars=resvars(Res,ExtVars)$
        !           184:                ResVars=append(vtol(ResVars),[1])$
        !           185: 
        !           186:                for(I=0;I<ExpMatRowNum;I++){
        !           187:                        TMP=0$
        !           188:                        for(J=0;J<ExpMatColNum;J++)
        !           189:                                TMP+=ExpMat[I][J]*ResVars[J]$
        !           190: 
        !           191:                        TMP-=ResVars[RevOneMat[I]+ExpMatColNum]$
        !           192: 
        !           193:                        if(TMP!=0)
        !           194:                                break$
        !           195:                }
        !           196: 
        !           197:                if(I==ExpMatRowNum){
        !           198:                        print("complitely homogenized")$
        !           199:                        return([makeret(Res,Vars,2)])$
        !           200:                }
        !           201:                else
        !           202:                        RET=cons(makeret(Res,Vars,1),RET)$
        !           203:        }
        !           204: 
        !           205:        ExpMat=qsort(ExpMat,afo)$
        !           206:        ExpMat2=[]$
        !           207:        for(I=0;I<size(ExpMat)[0];I++)
        !           208:                if(car(ExpMat2)!=ExpMat[I])
        !           209:                        ExpMat2=cons(ExpMat[I],ExpMat2)$
        !           210: 
        !           211:        ExpMat=newvect(length(ExpMat2),ExpMat2)$
        !           212:        ExpMatRowNum=size(ExpMat)[0]$
        !           213:        ExpMatColNum=size(ExpMat[0])[0]$
        !           214: 
        !           215:        NormMat=newmat(ExpMatColNum,ExpMatColNum+1)$
        !           216: 
        !           217:        for(I=0;I<ExpMatColNum;I++)
        !           218:                for(J=0;J<ExpMatColNum;J++)
        !           219:                        for(K=0;K<ExpMatRowNum;K++)
        !           220:                                NormMat[I][J]+=ExpMat[K][I]*ExpMat[K][J]$
        !           221: 
        !           222:        for(I=0;I<ExpMatColNum;I++)
        !           223:                for(J=0;J<ExpMatRowNum;J++)
        !           224:                        NormMat[I][ExpMatColNum]+=ExpMat[J][I]$
        !           225: 
        !           226:        SolveList=[]$
        !           227:        for(I=0;I<ExpMatColNum;I++){
        !           228:                TMP=0$
        !           229:                for(J=0;J<ExpMatColNum;J++)
        !           230:                        TMP+=NormMat[I][J]*Vars[J]$
        !           231: 
        !           232:                TMP-=NormMat[I][ExpMatColNum]$
        !           233:                SolveList=cons(TMP,SolveList)$
        !           234:        }
        !           235: 
        !           236:        Res=solve(SolveList,Vars)$
        !           237: 
        !           238:        if(nonposdegchk(Res))
        !           239:                RET=cons(makeret(Res,Vars,1),RET)$
        !           240: 
        !           241:        return reverse(RET)$
        !           242: }
        !           243: 
        !           244: end$

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