> 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)
}    #}0
> 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
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,]   -1    1    1    1   -1   -1   -1    1    1     1    -1    -1    -1     1
[2,]    1   -1    1    1   -1    1    1   -1   -1     1    -1    -1     1    -1
[3,]    1    1   -1    1    1   -1    1   -1    1    -1    -1     1    -1    -1
[4,]    1    1    1   -1    1    1   -1    1   -1    -1     1    -1    -1    -1