qs <- function(cont.table)
  {
    #author Ludwig Heigenhauser, Department of Statistics,
    #University of Florida, Gainesville
    #email: heigen@stat.ufl.edu
    time1 <- Sys.time()
    cont.table <- as.matrix(cont.table)
    if(dim(cont.table)[1] != dim(cont.table)[2])
      stop("This is not a square contingency table.")
    k <- dim(cont.table)[1]
    y <- as.vector(t(cont.table))

    sum.row <- rowSums(cont.table)
    sum.col <- colSums(cont.table)
    
    #Creating the design-matrix X
    #Parameters for the rows
    X.help <- kronecker(diag(1,k),matrix(1,k))
    X <- X.help

    #Parameters for the columns
    X.help <- kronecker(matrix(1,k),diag(1,k))
    X <- cbind(X,X.help[,-k])


    #Parameters for the interaction (There's probably an easier way to do it;-))
    X.help <- matrix(0,k*k,k*(k+1)/2)
    j <- 1
    to.extract <- NULL
    for(i in 1:k) 
      {
        to.extract <- c(to.extract,j)
        X.help[((i-1)*k+i):(i*k),j:(j+k-i) ] <- diag(k-i+1)
        j <- (j+k-i)+1
      }
    l.row <- NULL
    for(i in 1:(k-1))
      {
        for(j in 0:(k-i-1))
          l.row <- c(l.row,i*(k+1)+k*j)
      }
    l.col <- 1:(k*(k+1)/2)
    l.col <- l.col[-to.extract[-k]]

    for(i in 1:(k*(k-1)/2))
      {
         X.help[l.row[i],l.col[i]] <- 1
      }
    
    #take away the ones for the diagonal elements to make X of full column rank
    X.help <- X.help[,-(to.extract)] 
    X <- cbind(X,X.help)

    #fit a loglinear model for the data
    quasi.glm <- glm.fit(x=as.table(X), y=y, family=poisson())

    #Substract the diagonal elements from the margins
    sum.row <- sum.row-diag(cont.table)
    sum.col <- sum.col - diag(cont.table)

    #Take away the diagonal elements from the y and mean vector
    mu <- t(matrix(quasi.glm$fitted.values,k))
    y <- cont.table


    
    #get y and mu in the correct form for the C-function (row1~col1~row2 ...)
    y.c <- NULL
    mu.c <- NULL
    for(i in 1:(k-1))
    {
      y.c <- c(y.c,y[1,-1], y[-1,1])
      mu.c <- c(mu.c,mu[1,-1], mu[-1,1])

      mu <- mu[-1,-1]
      y <- y[-1,-1]
    }

    #sum of transposed opposite cells
    oppcell <- cont.table+t(cont.table)
    oppcell <- as.vector(oppcell)[as.vector(lower.tri(oppcell))]

    
    #Pearson Statistic for the observed table
    pears <- sum((y.c-mu.c)^2 / mu.c)

    #proportional probability for the table
    logfac <- function(x)
      {
        x[x==0] <- 1
        return(sum(-log(1:x)))
      }
    logfac.list <- function(x)
    {
      x <- as.list(x)
      lapply(x,logfac)
    }

    #table for computing the factorials
    fac.table <- c(0,unlist(logfac.list(1:max(oppcell))))

    reftab <- matrix(0,,k*(k-1))
    pos <- 0
    devup <- 0.0
    pearsup <- 0.0
    fisherup <- 0.0
    sumfac <- 0.0
    counts <- 0
    dcounts <- 0
    pcounts <- 0
    fcounts <- 0
   
    fisher <- sum(fac.table[y.c +1])

    # for deviance 
    dev <- y.c*(log(y.c) - log(mu.c)) + mu.c - y.c
    dev[dev=="NaN"] <- mu.c[dev=="NaN"]
    dev <- sum(dev)

    #Now we are ready to call the C-function
    Ccalled <- .C("subrec",
       as.integer(k),
       as.integer(sum.row),
       as.integer(sum.col),
       as.integer(oppcell),
       as.integer(reftab),
       as.integer(pos),
       as.double(dev),
       as.double(pears),
       as.double(fisher),
       as.double(devup),
       as.double(pearsup),
       as.double(fisherup),
       as.double(sumfac),
       as.integer(counts),
       as.double(mu.c),
       as.double(fac.table),
       as.integer(dcounts),
       as.integer(pcounts),
       as.integer(fcounts))
    
    time2 <- Sys.time()
    print(time2-time1)

    return(list(deviance.pval = Ccalled[[10]]/ Ccalled[[13]],
                pearson.pval = Ccalled[[11]]/ Ccalled[[13]],
                fisher.pval = Ccalled[[12]]/ Ccalled[[13]],
                counts = Ccalled[[14]],
                deviance.counts = Ccalled[[17]],
                pearson.counts = Ccalled[[18]],
                fisher.counts = Ccalled[[19]],
                deviance.asymptotic = 1 - pchisq(2*dev, k*k-dim(X)[2]),
                pears.asymptotic = 1 - pchisq(pears, k*k-dim(X)[2]) ))
  }