gc<-function(nn,r,s){ # gets total counts. nn = vector of n's #{0 lnn <- length(nn) # nn is a vector of n's. x's have s1 = r*n, s2 = s*n^2 gcout <- matrix(0,lnn,3) for(i in 1:lnn){ #{1 n <- nn[i] xi <- getxs(n,r,s) # matrix containing x's (x is a 4-ple with s1=rn,s2=sn^2) if(is.null(xi)) next if(length(xi)==4)xi <- matrix(xi,1,4) dxi <- dim(xi)[1] ci <- getcount(xi) # gets the sum of the counts for each x in xi gcout[i,] <- c(n,dxi,ci) } # c(n, number of x's, count) } #}1 gcout # matrix nnx3; rows are c(n,number of primitive sorted x's, total count) } #} getxs<-function(n,r=2,s=2,eps = 1e-6){ # gets primitive x's #{0 { eps2 <- eps^2 out <- NULL xs <- NULL rn <- r*n sn2 <- s*n^2 for(h in round(sqrt(sn2/4)):round(sqrt(sn2))) { #{1 sn21 <- sn2-h^2 if(sn21 < 0) next for(i in round(sqrt(sn21/3)):round(sqrt(sn21))){ #{2 # line 10 sn22 <- sn21 - i^2 if(sn22 <0 ) next for(j in round(sqrt(sn22/2)):round(sqrt(sn22))){ #{3 sn23 <- sn22 - j^2 if(sn23 < 0) next k <- sqrt(sn23+eps2) if(k%%1 < 2*eps) {k <- round(k)} else next x <- sort(c(h,i,j,k)) # x is a 4-ple with s2 = s*n^2 if(sum(x^2)!= sn2) next # just in case if(gcd(x)>1) next # retains only primitive x's # line 20 s1 <- sum(x) if(s1 < rn) next # no sign-flips will help xs <- rbind(xs,x) }}} #}1}2}3 # all x's have been found if(is.null(xs)) return(NULL) if(length(xs)==4)xs <- matrix(xs,1,4) xs <- unique(xs) # eliminate duplicates lxs <- length(xs)/4 if(lxs == 1)xs <- matrix(xs,1,4) for(i in 1:lxs) { #{4 outi <- NULL # line 30 xi <- xs[i,] # explore sign-flips if(sum(xi)==rn){outi<-matrix(xi,1,4)}else{ #{5 # xi is OK, no flips possible flipxi <- xi %*% flipmat # flipmat is 4x14 # line 30 for(j in 1:14) { #{6 #check all possible flips if(flipxi[j]==rn) outi <- rbind(outi,sort(xi*flipmat[,j])) }} #)5}6 if(is.null(outi)) next outi <- unique(outi) # eliminate duplicates out <- rbind(out,outi) } #}4 if(length(out) == 4) out <- matrix(out,1,4) out # matrix (4 cols) containing unique sorted primitive x's } # line 41 getcount<-function(X){ # gets total count from output of getxs(n,r,s) g <- c(0,24,12,6,4,0,0,1) tc <- 0 # total count if(is.null(X)) return(NULL) dX <- dim(X)[1] for(i in 1:dX){ #{1 xi <- X[i,] tabi <- table(xi) tc <- tc+g[sum(tabi^2)/2]} #}1 tc } gcd<-function(v){ # finds the g.c.d. of the elements of vector v. lv <- length(v) if(lv==0) return(0) if(lv==1)return(v) g <- gcd2(v[1],v[2]) if(lv==2)return(g) for(i in 3:lv) g <- gcd2(g,v[i]) g } gcd2<-function(a,b){ # gets the g.c.d. of two integers. a <- abs(a) b <- abs(b) if(a==0)return(b) x <- 1 while(x > 0){ x <- b%%a b <- a a <- x } b } flipmat<-matrix(c(-1,1,1,1,-1,-1,-1,1,1,1,-1,-1,-1,1,1,-1,1,1,-1,1,1,-1,-1,1,-1,-1,1,-1,1,1,-1,1,1,-1,1,-1,1,-1,-1,1,-1,-1,1,1,1,-1,1,1,-1,1,-1,-1,1,-1,-1,-1),nrow=4,byrow=T)