//		R o u t i n g   A l g o r i t h m s

//	Part 1 - Combinatorics and auxiliary functions

//	E. McCreight
//	last editedJune 15, 1977  12:14 PM  by emm

external
	[ 
	MoveBlock		// Defined by OS
	Zero
	CallSwat

	HeuristicNet	// Defined by other modules
	bestTotalNetLength
	Route		// Locally defined
	InitRandom
	Random
	GetOrderedRandomSet
	ArcLength
	ManhattanDistFn
	EuclideanDistFn

	n
	exhaustThresh
	Perm
	forceFirstNodeToEnd
	]


manifest
	[ memoTableEntries = 229	// prime
	]

static
	[ n
	Perm
	X
	Y
	forceFirstNodeToEnd

	randomTable
	randomIndex
	randomTrailer

	exhaustThresh = 7

	TrialPerm
	bestTotalNetLength

	randomInitialized = false

	DistFn
	memoTable
	]


structure MEMO:
	[ key word
	value word
	]


let Route(nNodes, localX, localY, ResultPerm, fFNTE, distanceMetricFn;
		numargs na) be

	[ if na ls 5 then forceFirstNodeToEnd = false

	test (na ls 6) % (distanceMetricFn eq 0)
		ifso	DistFn = ManhattanDistFn
		ifnot	DistFn = distanceMetricFn

	unless randomInitialized do InitRandom()

	n = nNodes
	X = localX
	Y = localY
	Perm = ResultPerm
	forceFirstNodeToEnd = fFNTE

	let mt = vec memoTableEntries*(size MEMO/16)
	memoTable = mt
	Zero(memoTable, memoTableEntries*(size MEMO/16))

	test nNodes le exhaustThresh

	ifso	[ bestTotalNetLength = #77776	// infinity
		let tP = vec 200
		TrialPerm = tP

		for i=1 to nNodes do TrialPerm!i = (i rem n)+1

		TryAllPermsRecursively(nNodes-1, 0)

		unless forceFirstNodeToEnd do
			for i=1 to nNodes-1 do
				[ TrialPerm!i = 1
				TrialPerm!nNodes = i+1

				TryAllPermsRecursively(nNodes-1, 0)

				TrialPerm!i = i+1
				]
		]

	ifnot HeuristicNet()

	if fFNTE then unless ResultPerm!nNodes eq 1 do
		CallSwat("Wrong node at end")
	]




//	Inter-node distance calculating functions.

and ArcLength(i, j) = valof

	[ if i eq j then resultis 0

	if i gr j then
		[ let t = i; i = j; j = t ]	// exchange i and j

	let memoKey = (i lshift 8)+j
	let memo = memoTable+
			((memoKey rem memoTableEntries) lshift 1)

	if memo>>MEMO.key eq memoKey then
		resultis memo>>MEMO.value

	memo>>MEMO.key = memoKey

	let value = DistFn(X!i, Y!i, X!j, Y!j)
	memo>>MEMO.value = value
	resultis value
	]


and ManhattanDistFn(X1, Y1, X2, Y2) = valof

	[ let deltaX = X1-X2
	let deltaY = Y1-Y2
	resultis ((deltaX ls 0)? -deltaX, deltaX)+
		((deltaY ls 0)? -deltaY, deltaY)
	]


and EuclideanDistFn(X1, Y1, X2, Y2) = valof

	[ let deltaX =  X1-X2
	let deltaY = Y1-Y2
	deltaX = (deltaX ls 0)? -deltaX, deltaX
	deltaY = (deltaY ls 0)? -deltaY, deltaY

	let bigDelta = nil
	let smallDelta = nil

	test deltaX ge deltaY
	ifso	[ bigDelta = deltaX
		smallDelta = deltaY
		]
	ifnot	[ bigDelta = deltaY
		smallDelta = deltaX
		]

	let ftIndex = (smallDelta ge #3777)?
		smallDelta/(bigDelta rshift 4),
		(smallDelta lshift 4)/bigDelta

	let fractionTable = table
		[ // entryK = 32*(sec(atan(K/16))-1)
		0; 0; 0; 1; 1; 1; 2; 3; 4; 5; 6; 7; 8; 9; 11; 12; 13
		]

	let fraction = fractionTable!ftIndex

	resultis bigDelta+(fraction*((bigDelta+16) rshift 4))
	]


//	1. The combinatorial algorithm results in a true
//	optimum but is prohibitively expensive for large nets.
//	The idea is that there are two vectors X and Y holding the
//	X and Y co-ordinates of a set of nodes (1...n), and a vector
//	TrialPerm (1...n) containing a permutation of the integers 1...n.
//	The outer program sets bestTotalNetLength to infinity
//	and then calls TryAllPerms(n-1, 0) with each node in
//	last place in TrialPerm.


and TryAllPermsRecursively(nNodesRemaining, netLengthSoFar) be

	[ if netLengthSoFar ge bestTotalNetLength then return

	if nNodesRemaining eq 0 then
		RecordBetterNet(netLengthSoFar)

	for i=nNodesRemaining to 1 by -1 do
		[ let t = TrialPerm!nNodesRemaining
		TrialPerm!nNodesRemaining = TrialPerm!i
		TrialPerm!i = t

		TryAllPermsRecursively(nNodesRemaining-1,
			netLengthSoFar+
			ArcLength(TrialPerm!(nNodesRemaining+1),
				TrialPerm!nNodesRemaining))

		t = TrialPerm!nNodesRemaining
		TrialPerm!nNodesRemaining = TrialPerm!i
		TrialPerm!i = t
		]
	]


and RecordBetterNet(length) be

	[ bestTotalNetLength = length
	MoveBlock(lv (Perm!1), lv (TrialPerm!1), n)
	]



//	This random number generator derives from the answer to
//	exercise 3.2.2-11 in the first edition of the second volume
//	of Knuth's "Art of Computer Programming." From Appendix C
//	in Peterson & Weldon's "Error-Correcting Codes" we learn
//	that x^33+x^13+1 is a primitive polynomial over GF(2). Thus
//	the sequence X(n) = (X(n-33)+X(n-13)) mod 2^16 has a period
//	length greater than 2^33 if the first 33 elements are not all
//	even.

and InitRandom() be

	[ manifest
		[ degree = 33
		midPower = 13
		]

	let foo = table [ 0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0; ]
	randomTable = foo
	randomTable!0 = #101011
	randomIndex = 0
	randomTrailer = degree-midPower

	for i=1 to 2000 do Random(1)
	randomInitialized = true
	]


and Random(max) = valof

	[ manifest
		[ degree = 33
		midPower = 13
		]

	let result = randomTable!randomIndex+
			randomTable!randomTrailer
	randomTable!randomIndex = result

	test randomIndex eq degree-1

	ifso	[ randomIndex = 0
		randomTrailer = degree-midPower
		]

	ifnot	[ randomIndex = randomIndex+1
		test randomTrailer eq degree-1
		ifso	randomTrailer = 0
		ifnot	randomTrailer = randomTrailer+1
		]

	resultis (result & #77777) rem (max+1)
			// This "rem" introduces a slight non-
			//	randomness.
	]


and GetOrderedRandomSet(n, vector, lowerLimit, upperLimit) be

	[ for i=1 to n do
		[ let newValue = Random(upperLimit+1-
				lowerLimit-i)+lowerLimit
		for j=1 to i-1 do
			test vector!j le newValue
			ifso newValue = newValue+1
			ifnot	[ let t = newValue
				newValue = vector!j
				vector!j = t
				]

		vector!i = newValue
		]
	]