#SYMBOLIC COMPUTATION FOR STATISTICAL INFERENCE
#
#D.F. Andrews, University of Toronto
#
#COPYRIGHT D. F. Andrews 2004
msg <- "";
inmsg <- "";
inv<- ""
invinv<- ""
inline <- "";
TypeA <- c(" ","A(","","*","","*",")","*", " ");
TypeD <- c("D["," "," "," "," ","; "," ",", ","]");
TypeM <- c(" ","E[",""," ","","*","]"," *", " ");
TypeK <- c(" ","K["," "," "," ",", ","]"," "," ");
TypeZ <- c(" ","Z["," "," "," "," ","]"," ", " ");
TypeL <-c("{ ","{ ","{ ",", "," }"," "," }"," ", " }");
TypeR <- c("R");
hasEstimated <- FALSE;
Elist <-NULL;
Clist<- NULL;
Cend <- -1;
MPE <- function(l=NULL,m=1.0){list(List = l, Mult = m)}
MPEstr <- function(mpe){
    s <- "";
    if(((mpe$Mult != 1.0) || (is.null(mpe$List)))) s <- paste(s,mpe$Mult,sep="");
    if(((mpe$Mult != 1.0) && !(is.null(mpe$List)))) s <- paste(s," * ",sep="");
    count <- 0;
    if(! is.null(mpe$List)){ia <- c(mpe$List)
        for(i in 0:(length(ia)-1)){
            if( i == 0){
                s <- paste(c(s,ia[[i+1]]),sep = "",collapse=" ");
                count <- 1;
            }
            else{
                if(ia[[i+1]] == ia[[i]]){
                    count <- count + 1;
                }
                else{
                    if( count > 1){
                        s <- paste(s , count,sep = "^",collapse="**") ;
                    }
                    s <- paste(s , ia[[i+1]],sep = "*");
                    count <- 1;
                }
            }
#print(c(i,count));
        }
        if(count > 1){
            s <-paste(s , count,sep = "^",collapse = "**") ;
        }
    }
    return(s)
}
j <- MPE(list("B","B","C","C","n","n","u"),2.0)
MPEstr(j)
MPEtimes<- function(pein1,pein2){
    al<-c(pein1$List,pein2$List);
    while(("n" %in% al) && ("1/n" %in% al)){
        alid <- 1:length(al);
        aln <- min(alid[al == "n"]);
        alu <- min(alid[al == "1/n"]);
        al <- al[-c(aln,alu)];
    }
#print("MPEstr")
#print(inv)
#print(invinv)
#print(al)
#print(c((inv %in% al), (invinv %in% al)))

    if(nchar(inv) > 1){
        while((inv %in% al) && (invinv %in% al)){
            alid <- 1:length(al);
            aln <- min(alid[al == inv]);
            alu <- 0;
            for(j in 1:length(al)){
                if(al[[j]] == invinv){
                    alu <- j
                }
            }
#			alu <- min(alid[al == invinv]);
#print(al == inv)
#print(al == invinv)
            al <- al[-c(aln,alu)];
        }
    }
    als <- NULL;
    if(!((is.null(al)) || (length(al) == 0))){
        als <- as.list(sort(unlist(al)))
    }
    return(  MPE(als, pein1$Mult * pein2$Mult));
}
j2 <- MPE(list("B","B","C","C","n"),2.0)
j1 <- MPE(list("A","1/n"),3.0)
j1m <- MPE(list("A","u"),-3.0)
j3 <- MPEtimes(j1,j2)
MPEstr(j3)
keylist <- function(l){
    if(is.null(l) || (length(l) == 0)) {
        s <- as.character(list(NULL))
    }
    else{
        s <- "";
        s <- paste(s,as.character(list(l[[1]])));
        if(length(l) > 1){
            for(i in 2:length(l)){
                s <- paste(s,as.character(list(l[[i]])))
            }
        }
    }
    s
}
MPEM<-function(mpe=NULL){
    if(is.null(mpe)) return(NULL)
#	else return(list(KeyList=list(MPEstr(mpe)),MPEList=list(mpe)))
#else return(list(KeyList=list(as.character(list(mpe$List))),MPEList=list(mpe)))
        else return(list(KeyList=list(keylist(mpe$List)),MPEList=list(mpe)))
}
MPEMaddMPE<- function(mpemp, p) {
#print("MPEMaddMPE start");
#print(mpemp)
#print(p)
    #			ia <- c(p$List);
#			nia <- sort(ia);
#			pkey <- MPEstr(p);
#			pkey <- as.character(list(p$List));
    pkey <- keylist(p$List);
    newmpem <- mpemp;
    if( pkey %in% mpemp$KeyList){
#print("match found");
#print(length(mpemp$KeyList));
#print(length(mpemp$MPEList));
        gpe <- min((1:length(mpemp$KeyList))[mpemp$KeyList == pkey]);
#print("gpe");
#print(gpe);
        gpemult <- mpemp$MPEList[[gpe]]$Mult + p$Mult;
#print(mpemp$MPEList[[gpe]]$Mult)
        if(abs(gpemult) < 10^(-8)){
            gpemult <- 0
        }
        if ( gpemult !=  0){
            newmpem$MPEList[[gpe]]$Mult <- gpemult;
        }
        				else{
                                            newmpem$KeyList <- newmpem$KeyList[-gpe];
                                            newmpem$MPEList <- newmpem$MPEList[-gpe];
                                        }
        			}
    else {
        newmpem$KeyList <- append(pkey,newmpem$KeyList);
        newmpem$MPEList <- append(list(p),newmpem$MPEList);
    }
#print("MPEMaddMPE end");
    return(newmpem)
}
m0<- MPEM(j1)
m1 <- MPEMaddMPE(m0,j1)
m2 <- MPEMaddMPE(m1,j2)
m4 <- MPEMaddMPE(m0,j1)
m4 <- MPEMaddMPE(m0,j1)
m4
MPEMaddmpem<- function(mpemin, pem) {
    newmpem <- mpemin;
#print("MPEList");
#print(pem$MPEList);
    if(!((is.null(pem$MPEList)) || (length(pem$MPEList) == 0))) {
        for(i in (1:length(pem$MPEList))){
#print("MPEMaddpem");
#print(i);
#print(pem$MPEList);
            newmpem <- MPEMaddMPE(newmpem,pem$MPEList[[i]])
        }
    }
    newmpem
}
m12 <- MPEMaddmpem(m1,m2)
m1
m2
m12
MPEMmpemtimes<- function (mpem, pemin){
    pemout <- MPEM(MPEtimes(mpem$MPEList[[1]],pemin$MPEList[[1]]));
#print(pemout);
    for( i1 in 1:length(mpem$KeyList)){
        for( j1 in 1:length(pemin$KeyList)){
            if( !((i1 ==1) && (j1 == 1)) ){
#print("MPEMmpemtimes");
#print(c(i1,j1))
#print(c(length(mpem$KeyList),length(pemin$KeyList)))
                pemout <-
                MPEMaddMPE(pemout,MPEtimes(mpem$MPEList[[i1]],pemin$MPEList[[j1]]))
            }
        }
    }
    return(pemout)
}
m1m2 <- MPEMmpemtimes(m1,m2)
m1m2
MPEMstr<- function(mpem){
    s <- "";
    kset <- mpem$KeyList;
    spe <- MPE();
    if(length(kset) > 1){
        s <- paste(s , "( ",sep="");
    }
    if( length(kset) > 0){
        spe <- mpem$MPEList[[1]];
        
#THE FOLLOWING CONDITION NEEDS TO BE FIXED
        
        if(!((length(kset) == 1) && (spe$Mult ==  1.0) &&
             (is.null(spe$List)))){
            s <- paste(s , MPEstr(spe),sep="");
        }
    }
    llen <- 0;
    if(length(kset) > 1){
        for(i in 2:length(kset)){
            spe <- mpem$MPEList[[i]];
            if(spe$Mult > 0){
                s <-paste( s , " + \n\t  ");
            }
            else{
                s <-paste( s , " + \n\t ");
            }
            s <-paste( s  , MPEstr(spe));
            llen <- llen + 1;
            if(llen > 300){
                s <- paste( s , " + \n   ");
                llen <- 0;
            }
        }
        s <- paste( s , ")",sep="");
    }
    return(s);
}
MPEMstr(m1m2)

MPEMD<- function(m1){
    mped <- MPE(m=m1);
    mped$Mult <- m1;
    return(MPEM(mped));
} 
# ELEMENTS
listsort<- function(l){
    clist <- as.character(lapply(l,strng));
    lout<- as.list(l[order(clist)]);
    lout
}
EL <- function(List = L,Type = "L"){
    if((Type == "P") ||(Type == "S")){
        return(list(List=listsort(List),Type=Type))
    }
    else{
        return(list(List=List,Type=Type))
    }
}
#ELEMNTS -> STRING
strngtype <- function(l,typelist,subs = NULL,clean=F){
    if(length(l) ==0) {
        return("")
    }
    else{
        s <- "";
        if(!((length(l) == 1) && 
             ((typelist[[1]] == "(") || (typelist[[1]] =="Dot(")) &&
             (clean))){
            s <- paste(s,typelist[[1]],sep="")
        }
        s <- paste(s,strng(l[[1]],clean=clean),sep="");
        i <- 1;
        scut <- 0;
        swidth<- 40;
        while( i < length(l)){
            i <- i + 1
            s <- paste(s, typelist[[2]],sep="")
            if( (nchar(s)) > scut + swidth) {
                scut <- scut + swidth;
                s <- paste(s,"\n\t\t")
            }
            s <- paste(s, strng(l[[i]],clean=clean),sep="")
        }
        if(!((length(l) == 1) && 
             ((typelist[[1]] == "(") || (typelist[[1]] =="Dot(")) && (clean))){
            s <- paste(s,typelist[[3]],sep="")
        }
        if(!is.null(subs)){
            s <- paste(s,"[",sep="");
            for(i in 1:length(subs)){
                s <- paste(s,subs[i],sep="");
                if( i < length(subs)){
                    s <- paste(s,",",sep="")
                }
            }
            s <- paste(s,"]",sep="");
        }
        return(s)
    }
}
strng<- function(l,clean=F){
	a <- try(l$Type,TRUE)
    if((class(a) == "try-error") | (is.null(a))){
        return(l)
    }
    else{
#print("error print")
#print(l$Type)
# print(l)

       ltype<-l$Type;
        if((length(l$List) == 1) && (l$Type == "Z")){
            ltype <- "Y"
        }
        return(
               switch(ltype,
                      "P" = strngtype(l$List,list("("," * ",")"),clean=clean),
                      "p" = strngtype(l$List,list("("," * ",")"),l$SubScripts,clean=clean),
                      "S" = strngtype(l$List,list("Dot(",", ",")"),clean=clean),
                      "Q" = strngtype(l$List,list("Dot(",", ",")"),clean=clean),
                      "q" = strngtype(l$List,list("Dot(",", ",")"),l$SubScripts,clean=clean),
                      "b" = strngtype(l$List,list("(","%*%",")^-1"),clean=clean),
                      "A" = strngtype(l$List,list("A("," * ",")"),clean=clean),
                      "K" = strngtype(l$List,list("C(",",",")"),clean=clean),
                      "D" = strngtype(l$List,list("D(",",",")"),clean=clean),
                      "M" = strngtype(l$List,list("E("," * ",")"),clean=clean),
                      "m" = strngtype(l$List,list("E("," * ",")"),l$SubScripts,clean=clean),
                      "W" = strngtype(l$List,list("IE("," * ",")"),clean=clean),
                      "Y" = strngtype(l$List,list("Z(","*",")"),l$SubScripts,clean=clean),
                      "y" = strngtype(l$List,list("z(","*",")"),l$SubScripts,clean=clean),
                      "Z" = strngtype(l$List,list("Z(Dot(",", ","))"),l$SubScripts,clean=clean),
                      "z" = strngtype(l$List,list("z(Dot(",", ","))"),l$SubScripts,clean=clean),
                      "L"= strngtype(l$List,list("{",",","}")),
                      l$List)
               )
    }
}
#PARTITION ELEMENTS
PE <- function(l=NULL,m=MPEMD(1.0)){
    return(list(List=l,MPEM=m))
}

#PE -> STRING
PEstr<- function(pe,clean=F){
    s <- "";
    if(!is.null(pe$List)){
        s <- strng(pe$List,clean=clean)
    }
    ms <- MPEMstr(pe$MPEM)
        s2 <- "";
    if(ms != ""){
        if(ms[1] != '('){
            s2 <- paste( s2 , "(",sep="");
        }
    }
    s2 <- paste( s2,  ms,sep="");
    if(ms != ""){
        if(ms[1] != '('){
            s2 <- paste( s2 , ")",sep="");
        }
    }
    if((ms != "") && (s != "")){
        s2 <- paste(s2," * ",sep="");
    }
    s2 <- paste( s2 , s,sep="");
    if(!((ms == "") && ((is.null(pe$List)) || (s2 == "")))){
        return(s2);
    }
    else{
        return("1");
    }
}

indexof<- function(item, l){
    id <- -1;
    if(item %in% l){
#print("indexof");
#print(as.vector(l));
#print(item);
        id <- min((1:length(l))[as.vector(l) == item])
    }
    return(id)
}
# PRODUCTS of PE
petimes <- function(pein1, pein2){
    al <- NULL;
    if(!(is.null(pein1$List))){
        al <- pein1$List$List;
    }
    if(!(is.null(pein2$List))){
        al <- append(al, pein2$List$List);
    }
    return( PE(l = EL(al,pein1$List$Type),
               m = MPEMmpemtimes(pein1$MPEM,pein2$MPEM))) 
}


ptimes<- function(  pein1, pein2){
    al <- NULL;
    altype <- NULL;
#identdone <- 0;
#l1 <- length(pein1);
#l2 <- length(pein2);
    #if((l1 > 0) && (l2 > 0)){
#	t1 <- pein1$List$List
#}
    if(!(is.null(pein1$List$List))){
        al <- append(al,pein1$List$List);
        altype <- pein1$List$Type;
    }
    if(!(is.null(pein2$List$List))){
        al <- append(al,pein2$List$List);
        altype <- pein2$List$Type;
    }
    if(is.null(altype)){
        return( PE( m =MPEMmpemtimes(pein1$MPEM,pein2$MPEM)))
    }
    else{	return( PE(EL(List = al,Type=altype), m =
                     MPEMmpemtimes(pein1$MPEM,pein2$MPEM)));
    } 
}

p3times<- function(  pein1, pein2){
    al <- NULL;
    if(!(is.null(pein1$List))){
        al <- append(pein1$List,list(al[[1]]));
    }
    if(!(is.null(pein2$List))){
        al <- append(pein2$List,list(al[[1]]));
    }
    return( PE(l = al, m = MPEMmpemtimes(pein1$MPEM,pein2$MPEM)));
}
#LISTS OF PES
PEM <- function(pe = NULL){
    if(is.null(pe)){
        return(list(KeyList = NULL,PEList = NULL))
    }
    else{
        return(list(KeyList=list(strng(pe$List,clean=F)),PEList=list(pe)))
    }
}

#PRINTING LISTS IF PE
PEMstr<- function(pem){
    s <- "";
    if(length(pem$KeyList) > 0){
        s <- paste(s, PEstr(pem$PEList[[1]],clean=T),sep="");
    }
    llen <- 0;
    if(length(pem$KeyList) > 1){
        for(i in 2:length(pem$KeyList)){
            s <- paste(paste(s, " + \n  "),PEstr(pem$PEList[[i]],clean=T),sep="");
            llen <- llen + 1;
#	if(llen > 2) {
#		s <- paste(s, "\n  ");
#		llen <- 0;
#	}
        }
    }
    return(s)
}

#ADD PE TO LIST
PEMaddpe<-function(pemin, p) {
    	pem <- pemin;
    tmkey <- "";
    if(!is.null(p$List)){
        tmkey <- strng(p$List);
    }
#print("tmkey");
#print(tmkey);
#print(pem$KeyList);
    if(tmkey %in% pem$KeyList){
        idk <- indexof(tmkey,pem$KeyList);
        ogpe <- pem$PEList[[idk]];
#print(idk);
#print("ogpe");
#print(ogpe);
        gpe <- PE(ogpe$List,ogpe$MPEM);
#print("gpe");
#print(gpe);
#print(p$MPEM);
        gpe$MPEM <- MPEMaddmpem(gpe$MPEM,p$MPEM);
#print("gpe");
#print(gpe);
        kset <- gpe$MPEM$KeyList;
#		if(is.null(kset)){
    if(!(length(kset) > 0)){
        pem$PEList <- as.list(pem$PEList[-idk]);
#print("LENGTH PELIST");
#print(	length(pem$PEList));
#print(pem$PEList);
        if(length(pem$PEList) == 0) {
            pem$PEList <- NULL
        }
        pem$KeyList <- as.list(pem$KeyList[-idk]);
        if( length(pem$KeyList) == 0) {
            pem$KeyList <- NULL
        }
    }
    else {
        pem$PEList[[idk]] <- gpe;
#pem$KeyList <- append(pem$KeyList,tmkey);
    }
}
 else {
     pem$PEList <- append(pem$PEList,list(p));
     pem$KeyList <- append(pem$KeyList,list(tmkey));
	}
pem
    }



# ADDING PEMS

PEMaddpem <- function(pemin, pem) {
    pemout <- pemin;
    if(!is.null(pem$KeyList)){
        for(i in 1:length(pem$KeyList)){
            pemout <- PEMaddpe(pemout,pem$PEList[[i]])
        }
    }
    class(pemout) <- class(pem);
    pemout
}

#NOT IMPLEMENTED YET

PEMpemtimes0<- function(pemin,type){
    pe <- pemin$PEList[[1]];
    pe1 <- PE(EL(list(pe$List),type),pe$MPEM);
    pemout <- PEM();
    pemout <- PEMaddpe(pemout,pe1);
    i <- 1;
    while( i < length(pemin$KeyList)){
        i <- i + 1;
        pe <- pemin$PEList[[i]];
        pe3 <- PE(EL(list(pe$List),type),pe$MPEM);
        pemout <- PEMaddpe(pemout,pe3);
    }
    return(pemout);
}

#PRODUCTS OF PEMS

PEMpemtimes<- function(pem0, pemin){
    pemout <- PEM();
    for(i in 1:length(pem0$KeyList)){
        for(j in 1:length(pemin$KeyList)){
            pemout <- PEMaddpe(pemout,petimes(pem0$PEList[[i]],pemin$PEList[[j]]))
#print(c(i,j))
#print(pemout)
        }
    }
    class(pemout) <- class(pemin)
        pemout
}


# NOT IMPLEMENTED YET

PEMpmtimes<- function(pem0, pemin){
    pemout <- PEM();
    if((length(pem0$KeyList) > 0) && (length(pemin$KeyList)> 0)) {
        for(i in 1:length(pem0$KeyList)){
            for(j in 1:length(pemin$KeyList)){
#print("PEMpmtimes")
#print(c(i,j))
                pemout <- PEMaddpe(pemout,ptimes(pem0$PEList[[i]],pemin$PEList[[j]]))
#print(pemout)
            }
        }
    }
    class(pemout) <- class(pemin)
        pemout
}

#NOT IMPLEMENTED YET

PEMp3times<- function(pem0, pemin){
    pemout <- PEM();
    for(i in 1:length(pem0$KeyList)){
        for(j in 1:length(pemin$KeyList)){
            pemout <- PEMaddpe(pemout,p3times(pem0$PEList[[i]],pemin$PEList[[j]]))
#print(c(i,j))
#print(pemout)
        }
    }
    pemout
}

TRpe<- function(f,level, pe, levelc = 0){
    if(length(pe$List$List) == 0){
        return(PEM(pe))
    }
    else{
        if(levelc > 0) {
            return(f(pe,level,levelc))
        }
        else{
            if((level == 1) ){
                return(f(pe,levelc));
            }
            else{
                ia<- pe$List$List;
                fp <- NULL;
                lm <- level - 1;
                pe4 <- PE(  ia[[1]],pe$MPEM);
#print("pe4");
#print(pe4$List);
                pem <- TRpe(f,lm,pe4,levelc);
#print("pem1");
#print(pem$KeyList);
                pem <- PEMpemtimes0(pem,pe$List$Type);
#print("pem2");
#print(pem$KeyList);
                i <- 1;
                while(i < length(ia)){
                    i <- i + 1;
                    pe5 <- PE(ia[[i]],MPEMD(1.0));
#print("pem3");
#print(pem$KeyList);
                    pem <-
                        PEMpemtimes(pem,PEMpemtimes0(TRpe(f,lm,pe5,levelc),pe$List$Type));
#print("pem4"); #print(pem$KeyList);
                }
                return(pem)
            }
        }
    }
}

TRpem<- function(f,level,pemin,levelc = 0){
    fpem <- PEM();
    if(length(pemin$KeyList) > 0){
        for(i in 1:length(pemin$KeyList)){
            pe <- pemin$PEList[[i]];
            if(!(is.null(pe$List$List))){
                if(length(pe$List$List) > 0){
                    fpem <- PEMaddpem(fpem,TRpe(f,level,pe,levelc))
                }
                else{
                    fpem <- PEMaddpe(fpem,pe)
                }
            }
            else{
                fpem <- PEMaddpe(fpem,pe)
            }
        }
    }
    return(fpem)
}


TypeA<-c("P","A","P","R")
TypeVA<-c("Q","A","Q","R")
TypeZ<-c("P","Z","P","R")
TypeM<- c("P","M","P","R")
TypeMz<- c("P","M","z","R")
TypeK<-c("P","K","P","R")
TypeD <-c("P","D","P","R")

MakeTypePEl <- function(l,type ){
    lout1 <- NULL;
    for(i in 1:length(l)){
        lout2 <- NULL;
        for(j in 1:length(l[[i]])){
            lout3 <- NULL;
            for(k in 1:length(l[[i]][[j]])){
                lout3 <- append(lout3,list(EL(l[[i]][[j]][[k]],type[[4]])))
            }
            lout2 <- append(lout2,list(EL(lout3,type[[3]])))
        }
        lout1 <- append(lout1,list(EL(lout2,type[[2]])))
    }
    PEM(PE(EL(lout1,type[[1]]),MPEMD(1.0)))
}

MakeTypePE <- function(pe,type ){
    l <- pe$List$List
    lout1 <- NULL;
    for(i in 1:length(l)){
        lout2 <- NULL;
        for(j in 1:length(l[[i]]$List)){
            lout3 <- NULL;
            for(k in 1:length(l[[i]]$List[[j]]$List)){
                lout3 <-
                append(lout3,list(EL(l[[i]]$List[[j]]$List[[k]]$List,type[[4]])))
            }
            lout2 <- append(lout2,list(EL(lout3,type[[3]])))
        }
        lout1 <- append(lout1,list(EL(lout2,type[[2]])))
    }
    PE(EL(lout1,type[[1]]),pe$MPEM)
}

MakeTypePEM<- function(pemin,type){
    fpem <- PEM();
    for(i in 1:length(pemin$KeyList)){
        pe <- pemin$PEList[[i]];
        if(!(is.null(pe$List$List))){
            if(length(pe$List$List) >0) {
                fpem <- PEMaddpem(fpem,PEM(MakeTypePE(pe,type)))
            }
            else{
                fpem <- PEMaddpe(fpem,pe)
            }
        }
        else{
            fpem <- PEMaddpe(fpem,pe)
        }
    }
    class(fpem) <- "P";
    return(fpem)
}

PEMQ <- function(pem, s, m){
    pem2 <- PEM();
    for(i in 1:length(pem$KeyList)){
        spe <- pem$PEList[[i]];
        spel2 <- spe$List$List;
        spem <- spe$MPEM;
        spel2 <- append(spel2,list(s));
        newpe <- PE(EL(spel2,spe$List$Type),spem);
        pem2 <- PEMaddpe(pem2,newpe);
        for(j in 1:length(spe$List$List)){
            spe2 <- NULL;
            for(k in 1:length(spe$List$List)){
                if(j == k){
                    spe3 <- spe$List$List[[j]];
#print("spe3");
#print(spe3);
                    spe3$List <- append(spe3$List,s$List);
#print("spe3");
#print(spe3);
                    spe3 <- listsort(spe3);
#print("spe3");
#print(spe3);
                    spe2 <- append(spe2,list(EL(spe3$List,"P")));
                }
                else{
#print("pemq")
#print(list(spe$List$List[[k]]))
                    spe2 <- append(spe2,list(spe$List$List[[k]]));
#print("spe2");
#print(spe2);
                }
            }
            newpe2 <- PE(EL(listsort(spe2),spe$List$Type),spem);
            pem2 <- PEMaddpe(pem2,newpe2);
        }
    }
    return(pem2);
}



FP<- function(  l, levelc = 0){
    ia <- l$List$List;
    lm <- l$MPEM;
    fp <- PEM();
    a1 <- NULL;
    a2 <- NULL;
    a2 <- append(a2,list(list(ia[[1]])));
    a1 <- append(a1,a2);
    pe1 <- PE(EL(list(EL(list(ia[[1]]),"P")),"P"),MPEMD(1.0));
    fp <- PEMaddpe(fp,pe1);
    i <- 1;
    while( i < length(l$List$List)){
        i <- i + 1;
#	for(i in 2:length(ia)){
#print("FP")
#print(strng(EL(list(ia[[i]]),"P")))
    fp <- PEMQ(fp,EL(list(ia[[i]]),"P"),l$MPEM);
}
	for(j in 1:length(fp$KeyList)){
            fp$PEList[[j]]$MPEM <- MPEMmpemtimes(fp$PEList[[j]]$MPEM,lm);
        }
	return(fp);
    }



U <- function(pe, levelc = 0){
    pemout <- PEM();
    lout <- NULL;
    j <- 0;
    while(j < length(pe$List$List)){
        j <- j + 1;
        lout <- append(lout,pe$List$List[[j]]$List);
#print("j");
#print(l);
    }
    pemout <- PEMaddpe(pemout,PE(l = EL(lout,"P") , m = pe$MPEM));
    pemout
}


SS <- function(pe, levelc = 0){
    pemout <- PEM();
    lout <-NULL;
    j <- 0;
    while(j < length(pe$List$List)){
        j <- j + 1;
        k <- 0;
        lout1 <- NULL;
        while( k < length(pe$List$List[[j]]$List)){
            k <- k + 1;
#print("S");
#print(c(j,k))
#print(pe$List$List[[j]]$List[[k]])
#print("EL")
#print(EL(list(pe$List$List[[j]]$List[[k]]),"P"));
            lout <- append(lout,list(EL(list(pe$List$List[[j]]$List[[k]]),"P")));
        }
#	lout <- append(lout,lout1);
    }
    pemout <- PEMaddpe(pemout,PE(l = (EL(lout,"P")) , m = pe$MPEM));
    pemout
}

len<- function( levell,  l){
    fp <- 0;
    if( levell == 1) {
        fp <- length(l$List);
    }
    else{
        llm <- levell - 1;
        i <- 0;
        while(i < length(l$List)){
            i <- i + 1;
            fp <- fp + len(llm,as.list(l$List[[i]]));
        }
    }
    return(fp);
}

cp<- function(i){
    if(i == 1) return(1)
        else return((-1)^(i-1)*gamma(i))
}
    cplc <- function(level,levelc,l){
        if(level == 1){
            if(levelc ==1 ){
                return(cp(len(levelc,l)));
            }
            else{
                cl <- 1.0;
                for(i in (1:length(l$List))){
                    cl <- cl * cp(len(levelc-1,l$List[[i]]))
                }
                return(cl);
            }
        }
        else{
            c2 <- 1;
            for(i in (1:length(l$List))){
                c2 <- c2 * cplc(level-1,levelc,l$List[[i]])
            }
            return(c2);
        }
    }
    CP<- function( pe, level, levelc){
        pemout <- PEM();
        	cout <- cplc(level,levelc,pe$List);
#print("CP");
#cat(strng(pe$List));
#print("length");
#print(len(levelc,pe$List));
#print(levelc);
#print(cout);
                pemout <-
                    PEMaddpe(pemout,PE(l=pe$List,m=MPEMmpemtimes(MPEMD(cout),pe$MPEM)));
                pemout
    }
    
    TRKM <- function(mupem){
        k0pem<-TRpem(SS,2,mupem);
        k1pem <- TRpem(FP,2,k0pem);
        k2pem<-TRpem(U,1,k1pem);
        k3pem <- MakeTypePEM(k2pem,TypeK);
        class(k3pem) <- "P";
        k3pem
    }
    
    
    TRMK<- function(k2pem){
        mk1 <- TRpem(FP,2,k2pem);
        mk2 <- TRpem(CP,1,mk1,levelc=2);
        mk3<- TRpem(U,1,mk2);
        mk4<-MakeTypePEM(mk3,TypeM);
        class(mk4) <- "P";
        mk4
    }
    
    lenc<- function( level,levell,  l){
#print("lenc, level,levell");
#print(c(level,levell));
#print(l);
        fp <- 0;
        if( level == 1) {
            fp <- len(levell,l);
        }
        else{
            llm <- level - 1;
            i <- 0;
            while(i < length(l$List)){
                i <- i + 1;
                fp <- fp + lenc(llm,levell,l$List[[i]]);
            }
        }
        return(fp);
    }
#lenc(2,1,list(list(list("a")),list(list("b"))))
    MPECA<- function(level, levelc, l){
#print("MPECA start");
        mpe <- MPE();
        i <- 0;
        li <- lenc(level,levelc,l);
        while( i < li){
            i <- i + 1;
            mpe$List <- append(mpe$List,"1/n");
        }
        return(mpe)
    }
    CA<- function( pe, level, levelc){
#print("CA start");
        pemout <- PEM();
        cout <- MPECA(level,levelc,pe$List);
        pemout <-
            PEMaddpe(pemout,PE(l=pe$List,m=MPEMmpemtimes(pe$MPEM,MPEM(cout))));
        pemout
    }
    #pab <- PE(list("a","b"),MPEMD(1.0))
    #MPECA(1,1,pab$List)
    
#t1<-TRpe(FP,1,pab)
#t1ca <- TRpem(CA,1,t1,levelc=2)
    MPECV<- function(level, levelc, l){
        mpe <- MPE();
        i <- 0;
        li <- lenc(level,levelc,l);
        while( i < li){
            i <- i + 1;
            mpe$List <- append(mpe$List,"n");
        }
        return(mpe)
    }
#pab <- PE(list("a","b"),MPEMD(1.0))
    #MPECV(1,1,pab$List)
    #pab <- PE(list("a","b"),MPEMD(1.0))
    #MPECV(1,1,pab$List)
    
    CV<- function( pe, level, levelc){
#print("CV level levelc")
#print(c(level,levelc))
        pemout <- PEM();
        cout <- MPECV(level,levelc,pe$List);
#print("CV cout");
#print(cout);
        pemout <-
            PEMaddpe(pemout,PE(l=pe$List,m=MPEMmpemtimes(pe$MPEM,MPEM(cout))));
        pemout
    }
    
    TRDA <- function(pema){
        pemd1 <- TRpem(FP,1,pema);
        pemd2 <- TRpem(U,3,pemd1);
        pemd3 <- TRpem(CA,2,pemd2,1);
        pemd4 <- TRpem(CV,1,pemd3,1);
        pemd5 <- MakeTypePEM(pemd4,TypeD);
        class(pemd5) <- "P";
        pemd5
    }
    #TRDA(mupem)
    #pemaa <- PEM()
#pemaa <-
#PEMaddpe(pemaa,PE(list(list(list("a")),list(list("b"))),MPEMD(1.0)))
    #cat(PEMstr(pemaa,TypeA))
    TRAD<- function(pemaa){
        pemd1 <- TRpem(FP,1,pemaa);
        pemd2 <- TRpem(U,3,pemd1);
        pemd3 <- TRpem(CA,2,pemd2,1);
        pemd31 <- TRpem(CP,1,pemd3,2);
        pemd4 <- TRpem(CV,1,pemd31,1);
        pemd5 <- TRpem(SS,2,pemd4);
        pemd6 <- MakeTypePEM(pemd5,TypeA);
        class(pemd6) <- "P";
        pemd6
    }
    
    cm<- function( j){
        mpem <- MPEMD( 1.0);
        i <- 2;
        while(i < (j+1)){
            i <- i + 1;
            mpem2 <- MPEMD( 1.0);
            mpe <- MPE(list("1/n"), ( -(i-2)) );
            mpem2<- MPEMaddMPE(mpem2,mpe);
            mpem <- MPEMmpemtimes(mpem,mpem2);
        }
        return(mpem);
    }
    MPECM<- function(level, levelc, l){
        mpem <- MPEM();
        if(level == 1){
#		mpem <- cm(len(levelc,l$List));
            mpem <- cm(len(levelc,l));
        }
        else{
            i <- 0;
            li <- length(l$List);
            while( i < li){
                i <- i + 1;
                if(i == 1) {
                    mpem <- cm(len(levelc,l$List[[1]]));
                }
                else{
                    mpem <- MPEMmpemtimes(mpem,cm(len(levelc,l[[i]])))
                }
            }
        }
#print(mpem);
        return(mpem)
    }
    
    CM<- function( pe, level, levelc){
        pemout <- PEM();
        cout <- MPECM(level,levelc,pe$List);
        pemout <-
            PEMaddpe(pemout,PE(l=pe$List,m=MPEMmpemtimes(pe$MPEM,cout)));
        pemout
    }
    
    TRMD <- function(pema){
        pemd1 <- TRpem(CM,1,pema,1);
        pemd2 <- TRpem(SS,2,pemd1);
        pemd3 <- MakeTypePEM(pemd2,TypeM);
        class(pemd3) <- "P";
        pemd3
    }
    
    cw<- function( j){
        mpem <- MPEMD( 1.0);
        i <- 2;
        while(i < (j+1)){
            i <- i + 1;
            mpem2 <- MPEMD( 1.0);
            mpe <- MPE(list(paste("1/(n-",i-2,")",sep = "")), ( (i-2)) );
            mpem2<- MPEMaddMPE(mpem2,mpe);
            mpem <- MPEMmpemtimes(mpem,mpem2);
        }
        return(mpem);
    }
    
    MPECW<- function(level, levelc, l){
        mpem <- MPEM();
        if(level == 1){
            mpem <- cw(len(levelc,l));
        }
        else{
            i <- 0;
            li <- length(l$List);
            while( i < li){
                i <- i + 1;
                if(i == 1) {
                    mpem <- cw(len(levelc,l[[1]]));
                }
                else{
                    mpem <- MPEMmpemtimes(mpem,cw(len(levelc,l$List[[i]])))
                }
            }
        }
#print(mpem);
        return(mpem)
    }
    CW<- function( pe, level, levelc){
        pemout <- PEM();
        cout <- MPECW(level,levelc,pe$List);
        pemout <-
            PEMaddpe(pemout,PE(l=pe$List,m=MPEMmpemtimes(pe$MPEM,cout)));
        pemout
    }
    #pemaa <- PEM()
#pemaa <-
#PEMaddpe(pemaa,PE(list(list(list("a")),list(list("b"))),MPEMD(1.0)))
    
    TRDM <- function(pema){
        pemd1 <- TRpem(CW,1,pema,1);
        pemd1 <- TRpem(SIMP,1,pemd1);
        pemd1 <- TRpem(SIMP,1,pemd1);
        pemd1 <- TRpem(SIMP,1,pemd1);
        class(pemd1) <- "D";
        pemd1
    }
    
    #cat(PEMstr(TRDA(pemaa),TypeD))
    #cat(PEMstr(TRMD(TRDA(pemaa)),TypeM))
    #cat(PEMstr(TRDM(TRMD(TRDA(pemaa))),TypeD))
    MPESIMP <- function(mpe){
        mpemsimp <- MPEM();
        if(("1/n" %in% mpe$List) && ("1/(n-1)" %in% mpe$List)){
            iall <- 1 : length(mpe$List);
            i1 <- min(iall[mpe$List == "1/n"]);
            i2 <- min(iall[mpe$List == "1/(n-1)"]);
            mpemsimp<-MPEMaddMPE(mpemsimp,MPE(mpe$List[-i1],mpe$Mult));
            mpemsimp<-MPEMaddMPE(mpemsimp,MPE(mpe$List[-i2],-1*mpe$Mult));
        }
        else{
            if(("1/n" %in% mpe$List) && ("1/(n-2)" %in% mpe$List)){
                iall <- 1 : length(mpe$List);
                i1 <- min(iall[mpe$List == "1/n"]);
                i2 <- min(iall[mpe$List == "1/(n-2)"]);
                mpemsimp<-MPEMaddMPE(mpemsimp,MPE(mpe$List[-i1],(1/2)*mpe$Mult));
                mpemsimp<-MPEMaddMPE(mpemsimp,MPE(mpe$List[-i2],-(1/2)*mpe$Mult));
            }
            else{
                if(("1/(n-1)" %in% mpe$List) && ("1/(n-2)" %in% mpe$List)){
                    iall <- 1 : length(mpe$List);
                    i1 <- min(iall[mpe$List == "1/(n-1)"]);
                    i2 <- min(iall[mpe$List == "1/(n-2)"]);
                    mpemsimp<-MPEMaddMPE(mpemsimp,MPE(mpe$List[-i1],mpe$Mult));
                    mpemsimp<-MPEMaddMPE(mpemsimp,MPE(mpe$List[-i2],-1*mpe$Mult));
                }
                
                else{
                    mpemsimp <- MPEMaddMPE(mpemsimp,mpe);
                }
            }
        }
        return(mpemsimp)
    }
    
    MPEMSIMP <- function(mpem){
        mpemsimp <- MPEM();
        i <- 0;
        while(i < length(mpem$MPEList)){
            i <- i + 1;
#print("MPEMSIMP i");
#print(i);
            mpemsimp <- MPEMaddmpem(mpemsimp,MPESIMP(mpem$MPEList[[i]]))
        }
        return(mpemsimp);
    }
    SIMP<- function( pe, level){
        pemout <- PEM();
        cout <- MPEMSIMP(pe$MPEM);
        pemout <-
            PEMaddpe(pemout,PE(l=pe$List,m=cout));
        pemout
    }
    #pemaa <- PEM()
#pemaa <-
#PEMaddpe(pemaa,PE(list(list(list("a")),list(list("a")),list(list("a"))),
#	MPEMD(1.0)))
    
    TRSIMP <- function(pema){
        pemout <- TRpem(SIMP,1,pema);
        class(pemout) <- class(pema);
        pemout
    }
#cat(PEMstr(TRDA(pemaa),TypeD))
    #cat(PEMstr(TRMD(TRDA(pemaa)),TypeM))
    #cat(PEMstr(TRDM(TRMD(TRDA(pemaa))),TypeD))
    #cat(PEMstr(TRSIMP(TRDM(TRMD(TRDA(pemaa)))),TypeD))
    #cat(PEMstr(TRSIMP(TRSIMP(TRDM(TRMD(TRDA(pemaa))))),TypeD))
    #cat(PEMstr(TRSIMP(TRSIMP(TRSIMP(TRDM(TRMD(TRDA(pemaa)))))),TypeD))
    
    TypeList <- c("P","A","L","K","M","Z")
    
    
    
    AV <- function(l){
        ape <- 	MakeTypePEl(l,TypeA)
        class(ape) <- "P";
        ape
    }
    EV <- function(l){
        ape <- 	MakeTypePEl(l,TypeM)
        class(ape) <- "P";
        ape
    }
    VAV <- function(l){
        ape <- 	MakeTypePEl(l,TypeVA)
        class(ape) <- "P";
        ape
    }
    CUM <- function(l){
        ape <- 	MakeTypePEl(l,TypeK)
        class(ape) <- "P";
        ape
    }
    print.P <- function(p){cat( PEMstr(p),"\n")}
    
    EVAL <- function(pema){
        eval(parse(text=PEMstr(pema)))
    }
    
    A <- function(x){ mean(x)}
    #Start 2004 stuff  
    #Code for univariate non asymptotic calculations
    PEMConst<- function(con){
        mpe <- MPE();
        if(is.numeric(con)){
            mpe$List<- NULL;
            mpe$Mult <- con
        }
        else{
            mpe$List <- list(con)
            mpe$Mult <- 1.0;
        }
        mpem <- MPEM();
        mpem <- MPEMaddMPE(mpem, mpe);
        pe <- PE(EL(NULL,"P"),mpem);
        pem <- PEM();
        pem <- PEMaddpe(pem, pe);
        class(pem) <- "P"
            pem
    }
    defined <- NULL
    getbit <- function(x,sarg = "F"){
#	bits <- NULL
        linear<- list("A","E","C","UE")
        if(length(x) > 1){
            if(as.character(x[[1]]) %in% linear){
                if(length(x[[2]]) >1){
                    if(x[[2]][[1]] == "+"){
                        j1 <- x
                        j2 <- x
                        j1[[2]]<- x[[2]][[2]]
                        j2[[2]]<-x[[2]][[3]]
                        bit1 <- getbit(j1,sarg)
                        bit2 <- getbit(j2,sarg)
                        if(class(bit2) == "S") {
                            bits <- SERadd(bit1,bit2)
                        }
                        else{
                            bits <- PEMaddpem(bit1,bit2)
                        }
                        gb <- bits
                        return(gb)
                    }
                }
                if(length(x[[2]]) >1){
                    if(x[[2]][[1]] == "*"){
                        if(is.numeric(x[[2]][[2]])){
                            j1 <- x
                            j2 <- x
                            j1<- x[[2]][[2]]
                            j2[[2]]<-x[[2]][[3]]
                            bits <- PEMpmtimes(getbit(j1,sarg),getbit(j2,sarg))
                            gb <- bits
                            return(gb)
                        }
                    }
                }
            }
        }
        if(length(x) == 1) {
            if(is.numeric(x)) {
                gb <- PEMConst(eval(parse(text=x)))
            }
            else {
                options(show.error.messages = FALSE)
                if(class(try(eval(x))) %in% 
                   list("P","S")) {
                    options(show.error.messages = TRUE)
                    return(eval(x))
                }
                else {
                    if(sarg == "F") {
                        gb <- PEMConst(x)
                    }
                    else {
                        gb <- list(x)
                    }
                }
                options(show.error.messages = TRUE)
            }
            undone <- F
        }
        else{
            undone <- T
            fn <- x[[1]]
            if(fn == "+"){
                arg1 <- getbit(x[[2]])
                arg2 <- getbit(x[[3]])
#			bits <- paste("add[",arg1,arg2,"]")
                if(class(arg2) == "S") {
                    bits <- SERadd(arg1,arg2)
                }
                else{
                    bits <- PEMaddpem(arg1,arg2)
                }
                class(bits) <- class(arg2)
                undone <- F
                gb <- bits
            }
            if(fn == "-"){
                if(length(x) == 3){
                    arg1 <- getbit(x[[2]])
                    arg2 <- getbit(x[[3]])
                    if(class(arg2) == "S") {
                        minus2 <- SERtimes(SERConst(-1,
                                                    length(arg2)),arg2)
                        class(minus2) <- class(arg2)
                        bits <- SERadd(arg1,minus2)
                    }
                    else{
                        minus2 <- PEMpmtimes(PEMConst(-1),arg2)
                        class(minus2) <- class(arg2)
                        bits <- TRSIMP( PEMaddpem(arg1,minus2))
                    }
                }
                if(length(x) == 2){
                    arg1 <- getbit(x[[2]])
                    if(class(arg1) == "S") {
                        bits <- SERtimes(SERConst(-1,
                                                  length(arg1)),arg1)
                        class(bits) <- class(arg1)
                    }
                    else{
                        minus2 <- PEMpmtimes(PEMConst(-1),arg1)
                        class(minus2) <- class(arg1)
                        bits <- TRSIMP(minus2)
                    }                
                }
                undone <- F
                gb <- bits
            }
            if(fn == "*"){
                arg1 <- getbit(x[[2]],sarg)
                arg2 <- getbit(x[[3]],sarg)
                if(sarg == "F"){
                    if(class(arg2) == "S") {
                        bits <- SERtimes(arg1,arg2)
                    }
                    else{
                        bits <- PEMpmtimes(arg1,arg2)
                    }
                    class(bits) <- class(arg2)
                }
                else{
                    if(length(arg1) == 1){
                        bits <- append(list(arg1),list(arg2))
                    }
                    else bits <- append(arg1,list(arg2))
                }
                    undone <- F
                    gb <- bits
            }
                if(fn == "A"){
                    if(length(x[[2]]) == 1){
                        arg1 <- getbit(x[[2]],sarg="A")
                        bits <- AV(list(list(arg1)))
                    }
                    else {
                        if((x[[2]][[1]] == "+") || (x[[2]][[1]] == "*")){
                            if(x[[2]][[1]] == "+"){
                                j1 <- x
                                j2 <- x
                                j1[[2]]<- x[[2]][[2]]
                                j2[[2]]<-x[[2]][[3]]
                                bits <- PEMaddpem(getbit(j1),getbit(j2))
                            }
                            else{
                                arg1 <- getbit(x[[2]],sarg="A")
                                bits <- AV(list(arg1))                        
                            }
                        }
                        else{
                            if(length(x[[2]]) == 2){
                                sfun <- as.character(x[[2]][[1]])
                                arg2 <- getbit(x[[2]][[2]])
                                bits <- Sfun(sfun,arg2,0,"A")
                            }
                            else{
                                if(length(x[[2]]) == 3){
                                    sfun <- as.character(x[[2]][[1]])
                                    arg2 <- getbit(x[[2]][[2]])
                                    bits <- Sfun(sfun,arg2,x[[2]][[3]],"A")                                
                                }
                                else{
                                    arg1 <- getbit(x[[2]],sarg="A")
                                    bits <- AV(list(arg1))
                                }
                            }
                        }
                    }
                    undone <- F
                    gb <- bits
                }
                if(fn == "InverseA"){
                    if(length(x[[2]]) == 2){
                        sfun <- as.character(x[[2]][[1]])
                        arg2 <- getbit(x[[2]][[2]])
                        bits <- Sroot(sfun,arg2,0,"A")
                    }
                    else{
                        if(length(x[[2]]) == 3){
                            sfun <- as.character(x[[2]][[1]])
                            arg2 <- getbit(x[[2]][[2]])
                            bits <- Sroot(sfun,arg2,x[[2]][[3]],"A")                                
                        }
                    }
                    undone <- F
                    gb <- bits
                }
                if(fn == "C"){
                    cl <- NULL
                    for(j in 2:(length(x))){
                        arg1 <- getbit(x[[j]],sarg="A")
                        if(length(x[[j]]) == 1){
                            cl <- append(cl,list(arg1))
                        }
                        else {
                            bitl <- NULL
                            for(k in 2:length(x[[j]])) bitl<- append(bitl,arg1[1][k])
                                cl <- append(cl,list(arg1))
                        }
                    }
                        bits <- CUM(list(cl))
                        undone <- F
                        gb <- bits
                }
                    if(fn == "E"){
                        arg1 <- getbit(x[[2]],sarg="A")
                        if(class(arg1) == "P"){
                            bits <- TRMD(TRDA(arg1))
                        }
                        else {
#				arg1 <- getbit(x[[2]],sarg="A")
                            if(length(x[[2]]) == 1){
                                bits <- EV(list(list(arg1)))
                            }
                            else {
                                bits <- EV(list(arg1))
                            }
                        }
                        undone <- F
                        gb <- bits
                    }
                    if(fn == "AE"){
                        arg1 <- getbit(x[[2]])
                        bits <- TRSIMP(TRSIMP(TRAD(TRDM(arg1))))
                        undone <- F
                        gb <- bits
                    }
                    if(fn == "EA"){
                        arg1 <- getbit(x[[2]])
                        bits <- TRMD(TRDA(arg1))
                        undone <- F
                        gb <- bits
                    }
                    if(fn == "EC"){
                        arg1 <- getbit(x[[2]])
                        bits <- TRMK(arg1)
                        undone <- F
                        gb <- bits
                    }
                    if(fn == "CE"){
                        arg1 <- getbit(x[[2]])
                        bits <- TRKM(arg1)
                        undone <- F
                        gb <- bits
                    }
                    if(fn == "BE"){
                        arg1 <- getbit(x[[2]])
                        bits <- MakeTypePEM(arg1,TypeA)
                        undone <- F
                        gb <- bits
                    }
                    if(fn == "EZ"){
                        arg1 <- getbit(x[[2]])
                        bits <- SexpectZ(arg1)
                        class(bits) <- class(arg1)
                        undone <- F
                        gb <- bits
                    }
                    if(fn == "APPROX"){
                        arg1 <- getbit(x[[2]],sarg)
                        x2 <- x[[2]]
                        if(length(x2) > 1) {
                            if(x2[[1]] == "-") {
                                if(is.numeric(x2[[2]])) {
                                    x2 <- -x2[[2]]
                                }
                            }
                        }
                        arg2 <- getbit(x[[3]],sarg)
                        if((class(arg1) == "P") && is.numeric(x[[3]]) && (length(x[[2]]) > 1) && !is.numeric(x2)) {
                            bits <- SERAEP(as.integer(x[[3]]+1),arg1)
                            class(bits) <- "S"
                        }
                        else{
                            if(is.numeric(x[[3]]) && (!is.numeric(x2))) {
                                bits <- SERConst(as.character(x[[2]]),x[[3]]+1)
                                class(bits) <- "S"
                            }
                            else{
                                if(is.numeric(x[[3]]) && (is.numeric(x2))) {
                                    bits <- SERConst(x2,x[[3]]+1)
                                    class(bits) <- "S"
                                }
                                else{
                                    options(show.error.messages = FALSE)
                                    if(is.numeric(try(eval(x[[3]]))) && (is.numeric(x2))) {
                                        bits <- SERConst(x2,eval(x[[3]])+1)
                                        class(bits) <- "S"
                                    }
                                    else{
                                        bits <- arg1
                                    }
                                    options(show.error.messages = TRUE)
                                }
                            }
                        }
                        undone <- F
                        gb <- bits
                    }                        
        }
                if(undone){
                    if(length(x)  > 1){
                        arg1 <- getbit(x[[2]])
                        if(class(arg1) == "S"){
                            if(fn == "/"){
                                bits <- SERtimes(arg1,Sfun("P(-1)",getbit(x[[3]]),0,"P"))
                                return(bits)
                            }
                            else{
                                if(fn == "sqrt"){ fn <- "P(1/2)"
                                }
                                if(fn == "^"){
                                    if(length(x) == 3){
                                        fn <- paste("P(",eval(x[[3]]),")",sep="")
                                    }
                                }
                            }
                            bits <- Sfun(fn,arg1,0,"P")
                            class(bits) <- class(arg1)
                            gb <- bits
                        }
                        else{
                            fnn <- fn
                        }
                    }
                    else{
                        bits <- paste(x[[1]],"(")
                        for(j in 2:length(x)) {
                            bits <- paste(bits,
                                          getbit(x[[j]]))
                        }
                        gb <- paste( bits,")")
                    }
                }
                gb
    }            
s <- function(x){
    getbit(parse(text=x)[[1]])
}
EA <- function(x) {
    TRMD(TRDA(x))
}
Const <- function(x){
    PEMConst(x)
}
Plus<- function(...){
    l <- (list(...));
    tot <- l[[1]];
    if(length(l) > 1) {
        for(i in 2:length(l)){
            tot <- PEMaddpem(tot,l[[i]])
        }
    }
    return(tot)
}

Times<- function(...){
    l <- (list(...));
    tot <- l[[1]];
    if(length(l) > 1) {
        for(i in 2:length(l)){
            tot <- PEMpmtimes(tot,l[[i]])
            class(tot)<- class(l[[i]])
        }
    }
    return(tot)
}
Eval <- function(xxx){
    if(class(xxx) == "P"){
        res <-eval(	parse(text=PEMstr(xxx)))
    }
    else{
        res <- 0
        for( i in 1:length(xxx)){
            if(!(PEMstr(xxx[[i]]) == "")) {
                res <- res + eval(	parse(text=PEMstr(xxx[[i]])))
            }
        }
    }
    res
}

#End 2004 stuff
SER <- function(i,serin = 0){
    l <- NULL;
    for(j in 1:i){
        tj <- PEM();
        if(!(serin == 0)){
            tj <- serin[[j]]
        }
        l <- append(l, list(tj));
    }
    class(l) <- "S";
    if(!(serin == 0)) {
        class(l) <- class(sein[[1]]);
    }
    return(l)
}
SERstr <- function(ser){
    s <- "";
    istart <- "";
    for(i in 1:length(ser)){
        if(!is.null(ser[[i]]$KeyList)) {
            s <- paste(s,istart,PEMstr(ser[[i]]))
            istart <- "+ \n"
        }
    }
    s <- paste(s,"\n");
    return(s);
}
print.S <- function(ser){
    cat(SERstr(ser))
}
stretchop <- function(serin,  tot,  len ){
    ser <- NULL;
    for(i in 1:tot){
        if((i <= len) && (i <= length(serin))){
            ser<- append(ser, list(serin[[i]]))
        }
        else{
            ser <- append(ser,list(PEM()))
        }
    }
    class(ser) <- "S";
    return(ser)
}
SERadd <- function(serin1,serin2){
    ser <-NULL;
    for(i in 1:length(serin1)){
        ser <- append(ser,list(PEMaddpem(serin1[[i]],serin2[[i]])))
    }
    class(ser) <- class(serin2);
    return(ser)
}
SERtimes <- function(serin1,serin2){
    ser <-NULL;
    ln <- length(serin1);
    for(i in 1:ln){
        pem <- PEM();
        for(j in 1 : i){
            if( (!is.null(serin1[[j]]$KeyList)) &&
                (!is.null(serin2[[i + 1 - j ]]$KeyList))){
                pem <- PEMaddpem(pem,PEMpmtimes(serin1[[j ]],
                                                serin2[[i + 1 - j ]]))
                class(pem) <- class(serin2[[i+1-j]]);
            }
        }
        ser <- append(ser,list(pem));
    }
    class(ser) <- class(serin1);
    return(ser)
}

SERpower <- function(serin, p){
    ser <- SER(length(serin));
    ser[[1]] <- PEMaddpe(PEM(),PE(EL(NULL,"P"),m=MPEMD(1.0)));
    if( p > 0) {
        for(i in 1:p){
            ser <- SERtimes(ser,serin)
        }
    }
    class(ser) <- class(serin);
    return(ser);
}

SERsctimes <- function(serin1,pem){
    ser <-SER(length(serin1));
    ln <- length(serin1);
    if(ln > 0){
        for(i in 1:ln){
#print("SERsctimes");
#print(serin1[[1]]);
            ser[[i]] <- PEMpmtimes(pem,serin1[[i]])
        }
    }
    class(ser) <- class(serin1);
    return(ser)
}
SERmpemtimes <- function(serin, mpem){
    ser <-SER(length(serin));
    for( i in  1:length(serin)){
        termi <-serin[[i]];
        if(length(termi$KeyList) > 0){
            for(j in 1:length(termi$KeyList)){
                if(!(is.null(termi$PEList[[j]]))){
                    termi$PEList[[j]]$MPEM <- MPEMmpemtimes(termi$PEList[[j]]$MPEM,mpem)
                }
            }
        }
        ser[[i]] <- termi;
    }
    class(ser) <- class(serin);
    return(ser);
}

SERaddpem<- function(serin,pemedm){
    ser <- serin;
#print(length(ser));
    sl <- length(ser);
    if(length(pemedm$PEList)>0){
        for(j in 1:length(pemedm$PEList)){
            pe <- pemedm$PEList[[j]];
            mult <- pe$MPEM$MPEList;
#print("mult");
#print(mult);
            if(length(mult)>0){
                for(k in 1:length(mult)){
                    mem <- mult[[k]]$List;
#print("mem");
#print(mem);
                    ordm <- ordmult(mem);
                    ord <- ordm;
                    if( (ord > -1) && (ord < sl)){
#Mar 11                        ser[[ord + 1]] <- PEMaddpe(ser[[ord+1]],pe)
                        newmpem <- MPEM()
                        newmpem<- MPEMaddMPE(newmpem,mult[[k]])
#print(PEstr(PE(l=pe$List,m=newmpem)))
                       ser[[ord + 1]] <- PEMaddpe(ser[[ord+1]],PE(l=pe$List,m=newmpem))
                    }
                }
            }
            else{
#print("addpe 1");
#print(pe);
                ser[[1]] <-  PEMaddpe(ser[[1]],pe);
            }
        }
    }
    class(ser) <- class(serin);
#print(length(ser));
    return(ser)
}
ordlist<- function( al){
    ord <- 0;
    if(length(al) > 0){
        for(i in 1:length(al)){
            ord <- ord + length(al);
        }
    }
    return(ord)
}
ordmult<- function( al){
    ord <- 0;
    if(length(al) > 0){
        for(i in 1:length(al)){
            if(al[[i]] == "1/n") {
                ord <- ord + 2
            }
        }
    }
#print("ord");
#print(ord);
            return(ord)
}
#NOTE SET UP FOR SINGLE AVERAGES ONLY
SERterma<- function( i,  pe,  type){
    ser <-SER(i);
    empe <- MPEstrtompe(pe, type);
    if( !( empe$Mult ==  0.0)){
        mpem <- MPEM();
        mpem <- MPEMaddMPE(mpem, empe);
        pe2 <- PE(m=mpem);
        ser[[1]] <- PEMaddpe(ser[[1]],pe2)
    }
    ser[[2]] <- PEMaddpe(ser[[2]],PE(pe$List,MPEMD(1.0)));
    class(ser) <- "S";
    return(ser)
}
SERterma<- function( i,  pe,  type){
    ser <-SER(i);
    empe <- pe;
    if(pe$List$List[[1]]$Type == "A"){
        empe$List$List[[1]]$Type <- "M"
    }
    ser[[1]] <- PEM(empe);
    class(ser[[1]]) <- "P";
    if(pe$List$List[[1]]$Type == "A"){
        zpe <- pe;
        zpe$List$List[[1]]$Type <- "Z"
            ser[[2]] <- PEM(zpe);
        class(ser[[2]]) <- "P";
    }
    class(ser) <- "S";
    return(ser)
}
SERterma<- function( i,  pe,  type){
    ser <-SER(i);
    empe <- pe;
    if(pe$List$List[[1]]$Type == "A"){
        pel <- pe$List$List[[1]]
        pel$Type<-"M"
        termstr <- strng(pel)
        empe<- PE(EL(NULL,"P"),MPEMmpemtimes(MPEM(MPE(list(termstr))),pe$MPEM))
    }
    ser[[1]] <- PEM(empe);
    class(ser[[1]]) <- "P";
    if(pe$List$List[[1]]$Type == "A"){
        zpe <- pe;
#Mar 6 changed Z t0 Y below
        zpe$List$List[[1]]$Type <- "Y"
            ser[[2]] <- PEM(zpe);
        class(ser[[2]]) <- "P";
    }
    class(ser) <- "S";
    return(ser)
}
SERAEP<- function( level,  pem){
    ser <-SER(level);
    if(length(pem$PEList) > 0) {
        for(i in 1:length(pem$PEList)){
            pe <- pem$PEList[[i]];
            ser0 <- SER(level);
            term0 <- PEM( PE(EL(NULL,"R"),m=pe$MPEM));
            class(term0) <- "P";
            ser0[[1]] <- term0;
            if(length(pe$List$List) > 0){
                for(j in 1:length(pe$List$List)){
                    pen <- PE(EL(list(pe$List$List[[j]]),pe$List$Type),MPEMD( 1.0));
                    ser0 <-SERtimes(ser0, SERterma(level,pen,TypeM));
                }
            }
            ser<- SERadd(ser,ser0);
        }
    }
    class(ser) <- "S";
    return(ser);
}

SERConst <- function(con,level){
    ser <- SER(level);
    mpe <- MPE();
    if(is.numeric(con)){
        mpe$List<- NULL;
        mpe$Mult <- con
    }
    else{
        mpe$List <- list(con)
        mpe$Mult <- 1.0;
    }
    mpem <- MPEM();
    mpem <- MPEMaddMPE(mpem, mpe);
    pe <- PE(EL(NULL,"P"),mpem);
    pem <- PEM();
    pem <- PEMaddpe(pem, pe);
    class(pem) <- "P"
        ser[[1]] <- pem;
    class(ser) <- "S";
    return(ser);
}

Splus<- function(...){
    l <- (list(...));
    tot <- l[[1]];
    if(length(l) > 1) {
        for(i in 2:length(l)){
            tot <- SERadd(tot,l[[i]])
        }
    }
    return(tot)
}

        Stimes<- function(...){
            l <- (list(...));
            tot <- l[[1]];
            if(length(l) > 1) {
                for(i in 2:length(l)){
                    tot <- SERtimes(tot,l[[i]])
                }
            }
            return(tot)
        }


derPE<- function( fun,  a,  i,ftype,PQ="Q"){
    derstr <- "";
    ma <- a[[1]]$PEList[[1]]$MPEM$MPEList[[1]]$Mult
        m <- 1
        ptf <- parse(text=as.character(fun))
        if(length(ptf[[1]]) > 1){
            if(ptf[[1]][[1]] == "P"){
                pow <- eval( ptf[[1]][[2]])
                if((pow-i) ==  1){
                    derstr <- paste( PEstr(a[[1]]$PEList[[1]]),sep="")
                }
                else{
                    if((pow -i) == 0) {
                        derstr <- NULL
                    }
                    else{
                        derstr <- paste( PEstr(a[[1]]$PEList[[1]]),"^(",pow-i,")",sep="")
                    }
                }
                if(i > 0) {
                    for(k in 1:i){
                        m <- m*(pow+1-k)
                    }
                }
            }
        }
    else{
        if(ptf[[1]] == "log"){
            if(i ==  0){
                derstr <- paste("log(", PEstr(a[[1]]$PEList[[1]]),")",sep="")
            }
            else{
                derstr <- paste("(",PEstr(a[[1]]$PEList[[1]]),")^-",i,sep="")
            }
            if(i > 1) {
                for(k in 2:i){
                    m <- m*(-(k-1))
                }
            }
        }
        else{
            if(ftype == "IP") {
                funstr <- paste("1/",fun)
            }
            else{
                funstr <- fun
            }
            if(i == 0){
                derstr <- paste( funstr  ,"(",
                                    PEstr(a[[1]]$PEList[[1]]),")",sep="")  ;
            }
            else{
                derstr <-paste( funstr  ,
                                "(",PEstr(a[[1]]$PEList[[1]]),",",i,")",sep="") ;
            }
        }
    }
#	type<-"P";
    if(length(a[[1]]$PEList[[1]]$List$List) >0){
        type <- a[[1]]$PEList[[1]]$List$List[[1]]$Type;
    }
    type <- ftype;
#	if(ftype == "A"){type<-"A"};
#print(type)
    if((ftype == "P") ||(ftype == "IP")){
        if(!(m ==0)){
            if(length(derstr) > 0){
                derpe <- PE(EL(NULL,"P"),MPEM(MPE(list(derstr),m)))
            }
            else{
                derpe <- PE(EL(NULL,"P"),MPEM(MPE(m=m)))
            }
        }
        else{
            derpe <- NULL
        }
    }
    else{
        all <- EL(list(EL(list(derstr),"R")),PQ);
        al <- EL(list(list(List =list(all),Type=type,Dims=i)),PQ);
        derpe <- PE(al,MPEMD(1.0));
    }
#	derpe <- PE(EL(list(al),"Q"),MPEMD(1.0));
#print(derpe);
    return(derpe);
}
Sder<- function( fun,  a,  i,ftype = "M",PQ="Q"){
    derpe <- derPE(fun,a,i,ftype,PQ=PQ);
    if((ftype == "P") || (ftype == "IP")){
        ser<- SER(length(a))
        if(!is.null(derpe)){
            ser[[1]]<- PEM(derpe)
        }
    }
    else{
        ser <- SERterma(length(a),derpe,ftype);
        if(((fun == "l") || (fun == "lo")|| (fun == "lon")) && (i == 1)){
            ser[[1]]<- PEM();	
        }
        if(( (fun == "lo")|| (fun == "lon")) && (i == 2)){
            ser[[1]]<- PEM(PE(m=MPEMD(-1)));	
        }
        if((fun == "lon") && ( i == 477)){
            ser[[2]] <-Splus(Sder("lo",a,i,ftype),
                                Stimes(SERConst(-1,length(a)), Sder("loi",a,i,ftype)))[[2]]
        }
    }
#if(length(ser[[1]]$KeyList )> 0) print(ser[[1]]$PEList[[1]]);
    return(ser);
}
Spetimes <- function(pein1, pein2){
    al <- NULL;
    if(!(is.null(pein1$List))){
        al <- pein1$List$List;
    }
    if(!(is.null(pein2$List))){
        al <- append(al, pein2$List$List);
    }
    return( PE(l = EL(al,pein1$List$Type),
                m = MPEMmpemtimes(pein1$MPEM,pein2$MPEM))) 
}
SPEMpemtimes<- function(pem0, pemin){
    pemout <- PEM();
    for(j in 1:length(pemin$KeyList)){
        pemout <- PEMaddpe(pemout,Spetimes(pem0$PEList[[i]],pemin$PEList[[j]]))
#print(c(i,j))
#print(pemout)
    }
    pemout
}
sym<- function(pein){
    symlist <- NULL;
    if(!is.null(pein$List)){
        symlist <- EL(list(pein$List),"S");
    }
    return(PE(l=symlist,pein$MPEM))
}
Ssympem<- function(pemin){
    pemout <- PEM();
    for(j in 1:length(pemin$KeyList)){
        pemout <- PEMaddpe(pemout,sym(pemin$PEList[[j]]))
#print(c(i,j))
#print(pemout)
    }
    pemout
}
termtimes1<-function(pemderiv,pemin){
    if(!is.null(pemderiv$PEList[[1]]$List)){
        el<- pemderiv$PEList[[1]]$List$List[[1]]$List[[1]]$List[[1]];
        etype <- pemderiv$PEList[[1]]$List$List[[1]]$Type;
        ele <- PE(EL(el$List,etype),pemderiv$PElist[[1]]$MPEM)
            lel <- list(EL(el$List,etype))
    }
    else{
        lel <- NULL;
        etype <- "C"
    }
            pemout <- PEM();
    for(j in 1:length(pemin$KeyList)){
        newl <- NULL;
        for(k in 1:length(pemin$PEList[[j]]$List$List)){
            itemlist <- pemin$PEList[[j]]$List$List[[k]]$List;
            for(l in 1:length(itemlist)){
                item <- itemlist[[l]];
                if(( k == 1) && (l == 1)) {
                    if((etype == "M") && (length(item$List) > 0)){
#print(el$List[[1]])
#print(item$List[[1]]$List[[1]])
#print((el$List[[1]]) ==  (item$List[[1]]$List[[1]]))
                        if((item$List[[1]]$Type == "W") && (item$List[[1]]$List[[1]] ==
                                                            el$List[[1]] )){
                            newl <- list(EL(
                                            item$List[2:length(item$List)],item$Type))
                        }
                        else{
                            newl <- list(EL(append(lel,
                                                    item$List),item$Type))
                        }
                    }
                    else{
                        newl <- list(EL(append(lel,
                                                item$List),item$Type))
                    }
                }
                else{
                    newl <- append(newl,list(item));
                }
            }
        }
        newpe <- PE(EL(newl,"Q"),
                    MPEMmpemtimes(pemderiv$PEList[[1]]$MPEM,pemin$PEList[[j]]$MPEM));
        pemout <- PEMaddpe(pemout,newpe)
    }
    pemout
}
termtimes2<-function(pemderiv,pemin){
    el<- pemderiv$PEList[[1]]$List$List[[1]];
    etype <- pemderiv$PEList[[1]]$List$List[[1]]$Type;
    ele <- PE(EL(el$List,etype),pemderiv$PElist[[1]]$MPEM)
        pemout <- PEM();
    for(j in 1:length(pemin$KeyList)){
        newl <- NULL;
        newl <- list(el);
        for(k in 1:length(pemin$PEList[[j]]$List$List)){
            itemlist <- pemin$PEList[[j]]$List$List[[k]]$List;
#				newl <- list(el);
            for(l in 1:length(itemlist)){
                item <- itemlist[[l]];
                newl <- append(newl,list(item));
            }
        }
        newpe <- PE(EL(newl,"Q"),
                    MPEMmpemtimes(pemderiv$PEList[[1]]$MPEM,pemin$PEList[[j]]$MPEM));
        pemout <- PEMaddpe(pemout,newpe)
    }
    pemout
}
Sfun <- function( fun,  a, i0=0,ftype="M"){ 
    al <- length(a);
    ser <- Sder(fun, a, i0,ftype);
    d <- SER(al);
    sa <- d;
    for(i in 1:al){
        if(length(a[[i]]$KeyList) > 0){
            if(ftype =="P"){
                sa[[i]] <- a[[i]]
            }
            else{
                sa[[i]] <- Ssympem(a[[i]])
            }
        }
    }
    for(i in 2:al){
        d[[i]]<-sa[[i]];
    }
    t<- SER(al);
# EL "S" changed ot "P" below
    t[[1]]<- PEM(PE(EL(NULL,"P"), MPEMD( 1.0)));
    class(t[[1]])<- "P";
    for(i in 2:al){
        t  <- Stimes(t,d);
        if( i > 1){
            pei <- PE(EL(NULL,"P"), MPEMD(( 1.0)/(i-1)));
            peim <- PEM(pei);
            class(peim) <- "P";
            t <- SERsctimes(t,peim);
        }
        deriv <- Sder(fun,a,i0+i-1,ftype);
        newterm <- SER(al);
        if(ftype == "P"){
            newterm <- SERtimes(deriv,t)
        }
        if(ftype != "P"){
            for( i in 2:al){
                if((length(t[[i]]$KeyList) > 0) && (length(deriv[[1]]$KeyList) > 0)){
                    if(ftype == "P"){
                        print(t)
                        print(deriv)
                        newterm[[i]] <- PEMpemtimes(deriv[[1]],t[[i]])
                    }
                    else{
                        newterm[[i]] <- termtimes1(deriv[[1]],t[[i]])
                    }
                }
                if(ftype != "P"){
                    if(length(t[[i-1]]$KeyList) > 0){
                        newterm[[i]] <- PEMaddpem(newterm[[i]],
                                                    termtimes2(deriv[[2]],t[[i-1]]))
                    }
                }
            }
        }
        ser <- SERadd(ser,newterm);
    }
    return(ser);
}
Srootinv <- function( fun,  a, i0){ 
    al <- length(a);
    derpe <- derPE(fun, a, i0 + 1);
    dpe <- (MPEstrtompe(derpe,TypeM));
    inval <- dpe$List;
    finv <- MPE();
    if(!(is.null(inval))){
        invinv <- inval[[1]];
        inv <- paste(invinv , "^(-1)",sep="");
        finv$List<- inv;
        finv$Mult <- -1.0
    }
    return(c(inv,invinv))
}

termtimes3<-function(pemderiv,pemin){
    if(!is.null(pemderiv$PEList[[1]]$List)){
        el<- pemderiv$PEList[[1]]$List$List[[1]]$List[[1]]$List[[1]];
        etype <- pemderiv$PEList[[1]]$List$List[[1]]$Type;
        ele <- PE(EL(el$List,etype),pemderiv$PElist[[1]]$MPEM);
    }
#edims <- pemderiv$PEList[[1]]$List$List[[1]]$Dims
    pemout <- PEM();
    for(j in 1:length(pemin$KeyList)){
        newl <- NULL;
        for(k in 1:length(pemin$PEList[[j]]$List$List)){
            item <- pemin$PEList[[j]]$List$List[[k]];
#print("termtimes3");
#print(pemderiv$PEList[[1]]$List$List[[1]]$Dims);
#print(item$Dims);
            if(( k == 1) &&(!is.null(pemderiv$PEList[[1]]$List)) ) {
                newl <- list(
                                #EL(append(list(EL(el$List,etype)),list(EL(item$List,"Q"))),item$Type));
                EL(append(list(EL(el$List,etype)),item$List),item$Type));
#list(List =append(list(EL(el$List,etype)),item$List),
#Type=item$Type,Dims=(edims + item$Dims - 2)));
                            }
else{
                                newl <- append(newl,list(item));
}
                    }
newpe <- PE(EL(newl,"Q"),
MPEMmpemtimes(pemderiv$PEList[[1]]$MPEM,pemin$PEList[[j]]$MPEM));
pemout <- PEMaddpe(pemout,newpe)
    }
pemout
}

Sroot <- function( fun,  a, i0,ftype){ 
    al <- length(a);
    ptype <- ftype;
    if(ftype == "A"){
        ptype <- "W"
    }
    if(ftype == "P"){
        ptype <- "b"
    }
    invser <- Sder(fun,a,i0 + 1,ptype);
    sm1 <- SER(length(a));
    sm1[[1]] <- PEM(PE(m=MPEMD(-1)));
    invser <- SERtimes(invser,sm1);
#print(invser);
    ser <- SER(1);
    ser[[1]] <- a[[1]];
    for( i in 2:al){
        ser2 <- SER(i);
        iser <- SER(i);
        iser[[1]] <- invser[[1]];
        for(j in 1:(i-1)){
            ser2[[j]] <- ser[[j]]
        }
        ser3 <- Sfun(fun,ser2,i0,ftype);
        for(j in 1:(i-1)){
            ser3[[j]]<- PEM()
        }
        ser3[[i]] <- termtimes3(invser[[1]],ser3[[i]]);
        ser2<- SERadd(ser2, ser3);
        ser <- ser2;
    }
    return(ser);
}
Sroot <- function( fun,  a, i0,ftype){ 
    al <- length(a);
    ptype <- ftype;
    if(ftype == "A"){
        ptype <- "IP"
    }
    if(ftype == "P"){
        ptype <- "IP"
    }
    invser <- Sder(fun,a,i0 + 1,ptype);
    sm1 <- SER(length(a));
    sm1[[1]] <- PEM(PE(m=MPEMD(-1)));
    invser <- SERtimes(invser,sm1);
#print(invser);
    ser <- SER(1);
    ser[[1]] <- a[[1]];
    for( i in 2:al){
        ser2 <- SER(i);
        iser <- SER(i);
        iser[[1]] <- invser[[1]];
        for(j in 1:(i-1)){
            ser2[[j]] <- ser[[j]]
        }
        ser3 <- Sfun(fun,ser2,i0,ftype);
        for(j in 1:(i-1)){
            ser3[[j]]<- PEM()
        }
        ser3[[i]] <- termtimes3(invser[[1]],ser3[[i]]);
        ser2<- SERadd(ser2, ser3);
        ser <- ser2;
    }
    return(ser);
}
Dot2 <- function(a1,a2){
    da1 <- dim(a1);
    da2 <- dim(a2);
    if(length(da1) ==1){
        c1 <- 1;
        m1 <- a1
    }
    else{	
        c1 <- da1[-length(da1)];
        m1 <- array(c(a1),dim=c(prod(c1),da1[length(da1)]));
    }
    if(length(da2) ==1){
        c2 <- 1;
        m2 <- a2
    }
    else{	
        c2 <- da2[-1];
        a2p <- aperm(a2,c(2:length(da2),1),resize=T);
        m2 <- array(c(a2p),dim=c(prod(c2),da2[1]))
    }
    if(length(da2) == 1){
        m1tm2 <- m1 %*% m2
    }
    else{
        m1tm2 <- m1 %*% t(m2);
    }
    if(c1 == 1){
        if( c2 == 1){
            c1c2 <- 1
        }
        else{
            c1c2 <- c2
        }
    }
    else{
        if( c2 == 1){
            c1c2 <- c1
        }
        else{
            c1c2 <- c(c1,c2)
        }
    }
    a1a2 <- array(c(m1tm2),c1c2);
    a1a2
}
Dot<- function(x1,x2=-37,...){
    if(x2 == -37){
        d <- x1
    }
    else{
        if(length(list(...)) > 0){
            d <-(Dot(Dot(x1,x2),...))
        }
        else{
            d <- (Dot2(x1,x2))
        }
    }
    d
}
E <- function(x){
    l <-  el;
    argstr <-deparse(substitute(x));
    eval(	parse(text=argstr))
}
Z <- function(x){
    l <-  zl;
    argstr <-deparse(substitute(x));
    zlv <-eval(	parse(text=argstr)) ;
    zlv 
}
#Compute the inverse
 IE <- function(x){
     argstr <-deparse(substitute(E(x)));
     em <-eval(	parse(text=argstr)) ;
#	em <- E(substitute(x));
     if(is.matrix(em)){
         emi <- solve(em,diag(rep(1,nrow(em))));
         return(emi)
     }
     else{
         return(1/em)
     }
 }
SEReval <- function(x){
    eval(parse(text=SERstr(x)))
}
dimstr<- function(s){
#print("dimstr");
#print(s);
    last <- nchar(s);
#print(last);
    	si <- substring(s,last-1,last-1);
#print(si);
        i <- as.numeric(si);
        if((is.na(i))){
            i <- 0
        }
#	print(i);
        i
}
dimlist <- function(l){
# print("dimlist in")
#	print(l);
    if(!(length(l)>0)){
        diml <- 0;
    }
    if(length(l) == 1){
        if(!(is.null(l[[1]]$List))){
            diml <- dimlist(l[[1]]$List)
        }
        else{
#			print("END")
#		print(l[[1]]);
            if(is.character(l[[1]])){
                diml <- dimstr(l[[1]])
            }
            else{
                diml <- 0
            }
        }
    }
    if(length(l) >1){
        ll <- length(l);
        diml <- dimlist(l[-ll]) + dimlist(l[ll])-2
    }
# print("dimlist out")
#	print(l);
#print(diml)
    return(diml)
}
indicies<-list("i","j","k","l","m","n","o",
               "p","q","r","s","t","u","v","w","x","y","z")
sscript<- function(l,used,free,d){
    if(!(d > 0)) {
        return(List = l, used = used, free = free)
    }
    else{
        if( length(used) > 0){
            useuse <- used[length(used)];
            used <- used[-length(used)];
            if(d > 1){
                usefree <- free[1:(d-1)];
                used <- append(used,usefree);
                free <- free[d: length(free)]
                    ss <- append(useuse,usefree);
            }
            else{
                free<- free[d+1:length(free)];
                ss <- useuse;
            }
        }
        else{
            if(d > 1){
                usefree <- free[1:d];
                used <- append(used,usefree);
                free <- free[d+1 : length(free)]
                    ss <- usefree;
            }
            else{
                used<-free[1];
                free<- free[d+1:length(free)];
                ss <- used[1];
            }
        }
        return(List=l,used=used,free=free,subscripts=ss)
    }
}
SubScript <- function(pein){
    plist <- pein$List$List;
    ptype <- pein$List$Type;
    if(length(plist) > 0){
        splist <- plist;
        free <-indicies;
        used <- NULL;
        for(i in 1:length(plist)){
            d <- dimlist(plist[[i]]$List);
#print("SS");
#print(d);
            splist <- sscript(plist[[i]]$List,used,free,d);
            used <- splist$used;
            free <- splist$free;
            item <- list(List = plist[[i]]$List,Type = plist[[i]]$Type,
                         SubScripts= splist$subscripts);
#print(splist$subscripts);
            plist[[i]] <- item
        }
    }
    if(ptype == "Q"){
        ptype <- "P"
    }
    elplist<- list(List=plist,Type=ptype);
    PE(elplist,pein$MPEM)
}
SubScriptPEM<- function(pemin){
    if(length(pemin$KeyList) > 0){
        pemout <- PEM();
        for(i in 1:length(pemin$KeyList)){
            pein <- pemin$PEList[[i]];
            pemout <- PEMaddpe(pemout,SubScript(pein))
        }
    }
    else{
        pemout <- pemin
    }
    return(pemout)
}

drop0 <- function(pemin){
    pemout <- PEM();
    if((length(pemin$PEList) > 0)){
        for(i in 1:length(pemin$PEList)){
            pe <- pemin$PEList[[i]];
            haszero <- F;
            plist <- pe$List$List;
            if(length(plist) > 0){
                for(j in 1:length(plist)){
                    if(length(plist[[j]]$List) == 1){
                        haszero <- T
                    }
                }
            }
            if(!haszero){
                pemout<- PEMaddpe(pemout,pe)
            }
        }
    }
    return(pemout)
}
drop1 <- function(pemin){
    pemout <- PEM();
    if((length(pemin$PEList) > 0)){
        for(i in 1:length(pemin$PEList)){
            pe <- pemin$PEList[[i]];
            haszero <- F;
            plist <- pe$List$List;
            if(length(plist)  == 1 ){
                   haszero <- T
            }
            if(!haszero){
                pemout<- PEMaddpe(pemout,pe)
            }
        }
    }
    return(pemout)
}
addl <- function(pemin){
    pemout <- PEM();
    if((length(pemin$PEList) > 0)){
        for(i in 1:length(pemin$PEList)){
            pe <- pemin$PEList[[i]];
            pemout<-PEMaddpe(pemout,
                             PE(EL(list(pe$List),"P"),pe$MPEM))
        }
    }
    return(pemout)
}


fixtype <- function(pemin,type){
    pemout <- PEM();
    if((length(pemin$PEList) > 0)){
        for(i in 1:length(pemin$PEList)){
            pe <- pemin$PEList[[i]];
            if(length(pe$List$List) >0){
                for(j in 1: length(pe$List$List)){
                    pe$List$List[[j]]$Type <- type
                }
            }
            pemout<-PEMaddpe(pemout,pe)
        }
    }
    return(pemout)
}
fixtype2 <- function(pemin,type){
    pemout <- PEM();
    if((length(pemin$PEList) > 0)){
        for(i in 1:length(pemin$PEList)){
            pe <- pemin$PEList[[i]];
            if(length(pe$List$List) >0){
                for(j in 1: length(pe$List$List)){
                    if(length(pe$List$List[[j]]$List) > 0){
                        for(k in 1:length(pe$List$List[[j]]$List)){
                            pe$List$List[[j]]$List[[k]]$Type <- type
                        }
                    }
                }
            }
            pemout<-PEMaddpe(pemout,pe)
        }
    }
    return(pemout)
}
TRMZ <- function(spem){
    tspem <- addl(spem);
    fpem <- TRpem(FP,1,tspem);
    uf <- TRpem(U,1,fpem);
    zpem <- drop0(uf);
    zpem1 <- TRpem(CA,1,zpem,2);
    zpem2 <- TRpem(CV,1,zpem1,1);
    fuzpem <- TRpem(FP,2,zpem2);
    cfuzpem <- TRpem(CP,1,fuzpem,levelc=2);
    uzcfzpem <- TRpem(U,1,cfuzpem);
    uz <- drop0(uzcfzpem);
    muz <- fixtype(uz,"M");
    muz2 <- fixtype2(muz,"z");
#Mar 6
    muz2 <- fixtype2(muz,"y");
    return(muz2)
}
#Mar 06
    TRDZ <- function(pema){
        pemd0 <- TRpem(U,2,fixtype2(pema,"y"));
        pemd2 <- TRpem(FP,1,pemd0);
#        pemd2 <- TRpem(U,3,pemd1);
        pemd3 <- TRpem(CA,2,pemd2,1);
        pemd4 <- TRpem(CV,1,pemd3,1);
        pemd5 <- MakeTypePEM(pemd4,TypeD);
        class(pemd5) <- "P";
        pemd5
    }

TRMZ <- function(spem){
    zpem <- (TRMD(TRDA(spem)));
    muz <- fixtype(zpem,"M");
    muz2 <- fixtype2(muz,"y");
    return(muz2)
}
TRMZ <- function(spem){
    zpem <- drop0(TRMD(TRDA(spem)));
    muz <- fixtype(zpem,"M");
    muz2 <- fixtype2(muz,"y");
    return(muz2)
}
TRMZ <- function(spem){
    zpem <- drop0(TRMD(TRDA(drop1(spem))));
    muz <- fixtype(zpem,"M");
    muz2 <- fixtype2(muz,"y");
    return(muz2)
}
#Mar 06
TRMZ <- function(spem){
    zpem <- drop0(TRMD(TRDZ(drop1(spem))));
    muz <- fixtype(zpem,"M");
    muz2 <- fixtype2(muz,"y");
    return(muz2)
}
#abandon MAR06 regression
TRMZ <- function(spem){
    tspem <- addl(spem);
    fpem <- TRpem(FP,1,tspem);
    uf <- TRpem(U,1,fpem);
    zpem <- drop0(uf);
    zpem1 <- TRpem(CA,1,zpem,2);
    zpem2 <- TRpem(CV,1,zpem1,1);
    fuzpem <- TRpem(FP,2,zpem2);
    cfuzpem <- TRpem(CP,1,fuzpem,levelc=2);
    uzcfzpem <- TRpem(U,1,cfuzpem);
    uz <- drop0(uzcfzpem);
    muz <- fixtype(uz,"M");
    muz2 <- fixtype2(muz,"z");
#Mar 6
    muz2 <- fixtype2(muz,"y");
    return(muz2)
}
#Apr 06
    TRMZD <- function(pema){
        pemd1 <- TRpem(CM,1,pema,1);
       # pemd2 <- TRpem(SS,2,pemd1);
        pemd3 <- MakeTypePEM(pemd1,TypeM);
        class(pemd3) <- "P";
        pemd3
    }

TRMZ <- function(spem){
    zpem <- drop0(TRMZD(TRDA(drop1(spem))));
    muz <- fixtype(zpem,"M");
    muz2 <- fixtype2(muz,"y");
    return(muz2)
}


Ssubscript<- function(serin){
    serout <- SER(length(serin));
    if(length(serin) > 0){
        for(i in 1:length(serin)){
            serout[[i]] <- SubScriptPEM(serin[[i]])
        }
    }
    class(serout) <- "S";
    return(serout)
}
SexpectZ<- function(serin){
    serout <- SER(length(serin));
    if(length(serin) > 0){
        serout[[1]] <- serin[[1]]
    }
    if(length(serin) > 1){
        for(i in 2:length(serin)){
            serouti <- TRMZ(serin[[i]]);
#print("SERZ");
#print(length(serout));
            serout <- SERaddpem(serout,serouti)
#serout[[i]] <- serouti
        }
    }
    class(serout) <- "S";
    return(serout)
}
SERterma<- function( i,  pe,  type){
    ser <-SER(i);
    empe <- pe;
    if(pe$List$List[[1]]$Type == "A"){
        pel <- pe$List$List[[1]]
        pel$Type<-"M"
        termstr <- strng(pel)
        empe<- PE(EL(NULL,"P"),MPEMmpemtimes(MPEM(MPE(list(termstr))),pe$MPEM))
    }
    ser[[1]] <- PEM(empe);
    class(ser[[1]]) <- "P";
    if(pe$List$List[[1]]$Type == "A"){
        zpe <- pe;
#Mar 6 changed Z t0 Y below
        zpe$List$List[[1]]$Type <- "Y"
            ser[[2]] <- PEM(zpe);
        class(ser[[2]]) <- "P";
    }
    class(ser) <- "S";
    return(ser)
}
derPE<- function( fun,  a,  i,ftype,PQ="Q",InvDer=F){
    sout <- SER(length(a))
    derstr <- "";
    ma <- a[[1]]$PEList[[1]]$MPEM$MPEList[[1]]$Mult
        m <- 1
        ptf <- parse(text=as.character(fun))
        if(length(ptf[[1]]) > 1){
            if(ptf[[1]][[1]] == "P"){
                pow <- eval( ptf[[1]][[2]])
                if((pow-i) ==  1){
                    derstr <- paste( PEstr(a[[1]]$PEList[[1]]),sep="")
                }
                else{
                    if((pow -i) == 0) {
                        derstr <- NULL
                    }
#                    else{
#                        if(!is.na(try(as.numeric(PEstr(a[[1]]$PEList[[1]])),silent=T))){
#                       derstr <- as.character( as.numeric(PEstr(a[[1]]$PEList[[1]]))^i)
#                       }
    else{
        derstr <- paste("(", PEstr(a[[1]]$PEList[[1]]),"^(",pow-i,")",")",sep="")
    }
#                    }
                }
                if(i > 0) {
                    for(k in 1:i){
                        m <- m*(pow+1-k)
                    }
                }
            }
        }
    else{
        if(ptf[[1]] == "log"){
            if(i ==  0){
                derstr <- paste("log(", PEstr(a[[1]]$PEList[[1]]),")",sep="")
            }
            else{
                derstr <- paste("((",PEstr(a[[1]]$PEList[[1]]),")^-",i,")",sep="")
            }
            if(i > 1) {
                for(k in 2:i){
                    m <- m*(-(k-1))
                }
            }
        }
        else{
            if(ftype == "P") {
                if( InvDer){
                    funstr <- paste("1/",fun,sep="")
                }
                else{
                    funstr <-fun
                }
            }
            else{
                funstr <- fun
            }
            if(i == 0){
                if(ftype == "A"){
                    if(InvDer){
                        derstr <- paste( "1/E(",funstr  ,"(",
                                         PEstr(a[[1]]$PEList[[1]]),")",")",sep="")
                    }
                    else{
                        derstr <- paste( "E(",funstr  ,"(",
                                         PEstr(a[[1]]$PEList[[1]]),")",")",sep="")
                    }
                }
                else{
                    if(InvDer){
                        derstr <- paste( "1/",funstr  ,"(",
                                         PEstr(a[[1]]$PEList[[1]]),")",sep="")
                    }
                    else{
                        derstr <- paste( funstr  ,"(",
                                         PEstr(a[[1]]$PEList[[1]]),")",sep="")
                    }
                    
                }
            }
            else{
                if(ftype == "A"){
                    if(InvDer){
                        derstr <- paste( "1/E(",funstr  ,"(",
                                         PEstr(a[[1]]$PEList[[1]]),",",i,")",")",sep="")
                    }
                    else{
                        derstr <- paste( "E(",funstr  ,"(",
                                         PEstr(a[[1]]$PEList[[1]]),",",i,")",")",sep="")
                    }
                }
                else{
                    if(InvDer){
                        derstr <- paste( "1/",funstr  ,"(",
                                         PEstr(a[[1]]$PEList[[1]]),",",i,")",sep="")
                    }
                    else{
                        derstr <- paste( funstr  ,"(",
                                         PEstr(a[[1]]$PEList[[1]]),",",i,")",sep="")
                    }
                }
            }
        }
    }
    
#	type<-"P";
    if(length(a[[1]]$PEList[[1]]$List$List) >0){
        type <- a[[1]]$PEList[[1]]$List$List[[1]]$Type;
    }
    type <- ftype;
#	if(ftype == "A"){type<-"A"};
#print(type)
    if((ftype == "P") ||(ftype == "IP") || (ftype == "A")){
        if(!(m ==0)){
            if(length(derstr) > 0){
                if(ftype == "A"){
                    derpe <- PE(EL(NULL,"P"),MPEM(MPE(list(paste(derstr,sep="")),m)))
                }
                else{
                    derpe <- PE(EL(NULL,"P"),MPEM(MPE(list(derstr),m)))
                }
            }
            else{
                derpe <- PE(EL(NULL,"P"),MPEM(MPE(m=m)))
            }
        }
        else{
            derpe <- NULL
        }
        sout[[1]] <- PEM(derpe)
    }
    if(ftype == "A"){
        if(i == 0){
            if(InvDer){
                derstr <- paste( "1/E(",funstr  ,"(",
                                 PEstr(a[[1]]$PEList[[1]]),")",")",sep="")
            }
            else{
                derstr <- paste( funstr  ,"(",
                                 PEstr(a[[1]]$PEList[[1]]),")",sep="")
            }
        }
        else{
            if(InvDer){
                derstr <- paste( funstr  ,"(",
                                 PEstr(a[[1]]$PEList[[1]]),",",i,")",sep="")
            }
            else{
                derstr <- paste( funstr  ,"(",
                                 PEstr(a[[1]]$PEList[[1]]),",",i,")",sep="")
            }
            
        }
        all <- EL(list(EL(list(derstr),"R")),PQ);
        al <- EL(list(list(List =list(all),Type="Y",Dims=i)),PQ);
        derpe <- PE(al,MPEMD(1.0));
        sout[[2]] <- PEM(derpe)
    }
#	derpe <- PE(EL(list(al),"Q"),MPEMD(1.0));
#print(derpe);
    return(sout);
}
Sder<- function( fun,  a,  i,ftype = "M",PQ="P",InvDer=F){
    ser <- derPE(fun,a,i,ftype,PQ=PQ,InvDer);
    if((ftype == "P") || (ftype == "IP") ||( ftype == "A")){
        ftype <- ftype
    }
#    else{
#       ser <- SERterma(length(a),derpe,ftype);
    if(((fun == "l") || (fun == "lo")|| (fun == "lon")) && (i == 1)){
        ser[[1]]<- PEM();	
    }
    if(( (fun == "lo")|| (fun == "lon")) && (i == 2)){
        ser[[1]]<- PEM(PE(m=MPEMD(-1)));	
    }
    if((fun == "lon") && ( i == 477)){
        ser[[2]] <-Splus(Sder("lo",a,i,ftype),
                         Stimes(SERConst(-1,length(a)), Sder("loi",a,i,ftype)))[[2]]
    }
#    }
#if(length(ser[[1]]$KeyList )> 0) print(ser[[1]]$PEList[[1]]);
    return(ser);
}
Spetimes <- function(pein1, pein2){
    al <- NULL;
    if(!(is.null(pein1$List))){
        al <- pein1$List$List;
    }
    if(!(is.null(pein2$List))){
        al <- append(al, pein2$List$List);
    }
    return( PE(l = EL(al,pein1$List$Type),
               m = MPEMmpemtimes(pein1$MPEM,pein2$MPEM))) 
}
SPEMpemtimes<- function(pem0, pemin){
    pemout <- PEM();
    for(j in 1:length(pemin$KeyList)){
        pemout <- PEMaddpe(pemout,Spetimes(pem0$PEList[[i]],pemin$PEList[[j]]))
#print(c(i,j))
#print(pemout)
    }
    pemout
}
Sfun <- function( fun,  a, i0=0,ftype="M"){ 
    al <- length(a);
    ser <- Sder(fun, a, i0,ftype);
    if(ftype == "A"){
        invser <- Sder(fun,a,i0 + 1,"A",InvDer=T);
        invinvser <- Sder(fun,a,i0 + 1,"A",InvDer=F)
#       inv <- invser[[1]]$PEList[[1]]$MPEM$MPEList[[1]]$List[[1]]
#       invinv <- invinvser[[1]]$PEList[[1]]$MPEM$MPEList[[1]]$List[[1]]
    }
    d <- SER(al);
    sa <- d;
    for(i in 1:al){
        if(length(a[[i]]$KeyList) > 0){
            if((ftype == "P" ) || (ftype =="A")){
                sa[[i]] <- a[[i]]
            }
            else{
                sa[[i]] <- Ssympem(a[[i]])
            }
        }
    }
    if(al > 1) {
        for(i in 2:al){
            d[[i]]<-sa[[i]];
        }
        t<- SER(al);
# EL "S" changed ot "P" below
        t[[1]]<- PEM(PE(EL(NULL,"P"), MPEMD( 1.0)));
        class(t[[1]])<- "P";
        for(i in 2:al){
            t  <- Stimes(t,d);
            if( i > 1){
                pei <- PE(EL(NULL,"P"), MPEMD(( 1.0)/(i-1)));
                peim <- PEM(pei);
                class(peim) <- "P";
                t <- SERsctimes(t,peim);
            }
            deriv <- Sder(fun,a,i0+i-1,ftype);
            if(ftype == "A"){
                invser <- Sder(fun,a,i0 +i-1,"A",InvDer=F);
                invinvser <- Sder(fun,a,i0 + i-1,"A",InvDer=T)
                    if(!is.null(invser[[1]]$PEList[[1]]$MPEM$MPEList[[1]]$List[[1]])){
                        assign("inv" , invser[[1]]$PEList[[1]]$MPEM$MPEList[[1]]$List[[1]],inherits=T)
                    }
                if(!is.null(invinvser[[1]]$PEList[[1]]$MPEM$MPEList[[1]]$List[[1]])){
                    assign("invinv" , invinvser[[1]]$PEList[[1]]$MPEM$MPEList[[1]]$List[[1]],inherits=T)
                }
            }
            newterm <- SER(al);
            if((ftype == "P" ) || (ftype =="A")){
#    print("Sfun")
#    print(inv)
#    print(invinv)
                newterm <- SERtimes(deriv,t)
            }
            if(!((ftype == "P" ) || (ftype =="A"))){
                print(ftype)
                for( i in 2:al){
                    if((length(t[[i]]$KeyList) > 0) && (length(deriv[[1]]$KeyList) > 0)){
                        if(ftype == "P"){
                            print(t)
                            print(deriv)
                            newterm[[i]] <- PEMpemtimes(deriv[[1]],t[[i]])
                        }
                        else{
                            newterm[[i]] <- termtimes1(deriv[[1]],t[[i]])
                        }
                    }
                    if(ftype != "P"){
                        if(length(t[[i-1]]$KeyList) > 0){
                            newterm[[i]] <- PEMaddpem(newterm[[i]],
                                                      termtimes2(deriv[[2]],t[[i-1]]))
                        }
                    }
                }
            }
            ser <- SERadd(ser,newterm);
        }
    }
    if(ftype == "A"){
        assign("inv",  "",inherits=T)
        assign("invinv","",inherits=T)
    }
    
    return(ser);
}
Sroot <- function( fun,  a, i0,ftype){ 
    al <- length(a);
    ptype <- ftype;
    if(ftype == "A"){
        ptype <- "IP"
    }
    if(ftype == "P"){
        ptype <- "IP"
    }
    invser <- Sder(fun,a,i0 + 1,ftype,InvDer = T);
    sm1 <- SER(length(a));
    sm1[[1]] <- PEM(PE(m=MPEMD(-1)));
    invser <- SERtimes(invser,sm1);
#print(invser);
    ser <- SER(1);
    ser[[1]] <- a[[1]];
    for( i in 2:al){
        ser2 <- SER(i);
        iser <- SER(i);
        iser[[1]] <- invser[[1]];
        for(j in 1:(i-1)){
            ser2[[j]] <- ser[[j]]
        }
        ser3 <- Sfun(fun,ser2,i0,ftype);
        for(j in 1:(i-1)){
            ser3[[j]]<- PEM()
        }
#        ser3[[i]] <- termtimes3(invser[[1]],ser3[[i]]);
        ser3[[i]] <- PEMpmtimes(invser[[1]],ser3[[i]]);
        ser2<- SERadd(ser2, ser3);
        ser <- ser2;
    }
    return(ser);
}
Z <- function(x){
    mean(x - E(x))
}
z<- function(x){
    x - E(x)
}

makefn <- function(xxx){
	    if(class(xxx) == "P"){
		pt <- parse(text=PEMstr(xxx))
        res <-function(){eval(pt)}
    }
    else{
        res <- 0
        for( i in 1:length(xxx)){
            if(!(PEMstr(xxx[[i]]) == "")) {
                res <- res + eval(	parse(text=PEMstr(xxx[[i]])))
            }
        }
    }
    res
}
S <- function(x){
    l <-  el;
    argstr <-deparse(substitute(x));
s(argstr)
}


lprod <- function(x){
		if(length(x) == 1) lp <- 1
		if(length(x) > 1){
			 lp <- 1 + lprod(x[[2]])
			}
		lp
	}
gm <- function(i){
	m <- 0
	if(!(i %%2 )) {
		if( i == 2) m <- 1
		if( i > 3)  m <- (i-1) * gm( i -2)
		}
	m
}

	
GaussMoment <- function(xxx){
	gm(lprod(substitute(xxx)))
	}

Sfn <- function(xxx,x,...){
	    if(class(xxx) == "P"){
		pt <- parse(text=PEMstr(xxx))
        res <-function(x,...){n <- length(x);eval(pt,list(X = x))}
    }
    else{
        res <- 0
        for( i in 1:length(xxx)){
            if(!(PEMstr(xxx[[i]]) == "")) {
                res <- res + eval(	parse(text=PEMstr(xxx[[i]])))
            }
        }
    }
    res
}