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