A simple R script containing functions to generate solutions to monic polynomials

Monic Polynomial Solver 1.0.r

MonicSolve = function(

	order = 5,				# the order of the monic polynomials being solved
	n = 5,					# n, in (n,n-1,...,-n), the acceptable coefficient values
	iterations = 1000,			# the number of equations to solve
	random = FALSE				# if m!= FALSE, use randomly generated polynomial coefficients
	
)
{

	# Ryan Grannell, 2012
	# A function to generate a set of eigenvalues

	pos = 1
	size = order - 1
	polynomial = matrix(0, size, size)

	vector = rep(-n, size)
	
	storeX = storeY = rep(0,
		size * order * iterations
	)
	
	for(i in 1:(size - 1)){
	
		polynomial[i+1,i] = 1
	
	}
	
	# 1. Generate each  set of coefficients for the
	# monic polynomials
	
	if(random == FALSE){
	
		polyCoefficients = AutoLoop(
			vector, 
			n, 
			iterations
		)
	
	}else{
	
		polyCoefficients = matrix(0, iterations, size)
	
		for(i in 1:iterations){
		
			polyCoefficients[i,] = sample(-n:n, 
				size, 
				replace = TRUE
			)
		
		}
	
	}
	
	startTime = Sys.time()
					
	for(i in 1:nrow(polyCoefficients)){
	
		#2. basic time estimation function
	
		if(i %/% 1000 == i/1000){
		
			timeTaken = difftime(Sys.time(), startTime, unit = "mins")
			estimate = round((timeTaken/(i)) * ((iterations - i)+1),3)
		
			print(
				paste("iteration", i, "of", iterations, "; ~ ",
				estimate, "minutes left")
			)
		
		}
		
		polynomial[,size] = -1* polyCoefficients[i,]
		
		#3. solve the polynomial and store its
		#   eigenvalues z. storeX = Re(z), storeY = Im(z)
		
		roots = eigen(polynomial,
			only.values = TRUE)$values
		
		storePositions = pos:(pos + (size - 1))
		
		storeX[storePositions] = Re(roots)
		storeY[storePositions] = Im(roots)
				
		pos = pos + order
	
	}
	
	storeX = storeX[1:(pos + (size-1))]
	storeY = storeY[1:(pos + (size-1))]
	
	return(
		rbind(storeX, storeY)
	)
	
}

AutoLoop = function(

	vector = c(-5,-5,-5,-5),	# vector, a vector representing an integer, each position containing a digit
	m = 5,				# m, in (m,m-1, ... -m), the acceptable values for a coefficient
	iterations = 1000		# iterations, the amount to be added to vector in base m arithmetic

)
{

	# Ryan Grannell, 2012
	# A function to generate the coefficients for MonicSolve()

	range = (2*m) + 1
	startVector = vector
	
	low <- rep(-m, length(vector))
	iterations = min(range^length(vector), iterations)
	
	store = matrix(0, iterations, length(vector))
	
	for(i in 1:iterations){
	
		j = length(vector)
	
		repeat{
		
			vector[j] = vector[j] + 1
		
			if(vector[j] > m && j != 1){
			
				vector[j] = low[j]
				j = j - 1
				
			}else break
		
		}
		
		store[i,] = vector

	}
	
	if(store[nrow(store),1] > m){
	
		store = store[-nrow(store),]
	
	}
	
	store = rbind(startVector, store)
	return(store)

}