//	NEBCheckerboardTilings.swift
//
//	© 2026 by Jeff Weeks
//	See TermsOfUse.txt

import Metal


extension NEBRenderer {

	//	Starting from the reflection groups generated by
	//	a right-angled spherical tetrahedron,
	//	a right-angled Euclidean cube and
	//	a right-angled hyperbolic dodecahedron,
	//	we'll compute only the non-mirror-reversed elements.
	//	In each tiling, the mirror-reversed and non-mirror-reversed
	//	tiles alternate checkboard style, so if the renderer
	//	draws the walls of the non-mirror-reversed copies of
	//	the tetrahedron, cube or dodecahedron with face culling disabled,
	//	that will show the walls of the mirror-reversed copies "for free".
	//
	//	We'll adjust the poses so that the central cell
	//	in each tiling sits directly over the pool table.
	//	That central cell will be the first element in each list;
	//	its 30 nearest neighbors will follow immediately after it.

	func makeCheckerboardPoseMatrices() -> [simd_double4x4] {

		switch spaceGeometry.geometryType {
		case .spherical:	return sphericalCheckerboardPoses()
		case .euclidean:	return euclideanCheckerboardPoses()
		case .hyperbolic:	return hyperbolicCheckerboardPoses()
		}
	}

	func fogCutoff(
		for poseMatrices: [simd_double4x4]
	) -> Double {	//	in radians

		let theMaxTranslationDistanceInRadians: Double
		switch spaceGeometry.geometryType {
		
		case .spherical:

			//	Caution: If the maximum translation distance is
			//	set to noticeably less than π, walls heading
			//	toward the south pole experience severe banding.
			//	Not with bands of similar colors, but
			//	with adjacent bands of very different colors.
			//	I don't understand why that is. As a practical
			//	matter, though, it makes sense to set
			//	the maximum translation distance to exactly π.
			//
			theMaxTranslationDistanceInRadians = Double.pi
			
		case .euclidean:
			let theMaxSquaredDistanceInRadians = poseMatrices.reduce(0.0,
			{ maxSquaredDistanceSoFar, matrix in
				max(
					  matrix[3][0] * matrix[3][0]
					+ matrix[3][1] * matrix[3][1]
					+ matrix[3][2] * matrix[3][2],
					maxSquaredDistanceSoFar
				)
			})
			theMaxTranslationDistanceInRadians = sqrt(theMaxSquaredDistanceInRadians)
		
		case .hyperbolic:
			
			let theMaxW = poseMatrices.reduce(1.0,
			{ maxWSoFar, matrix in
				max( matrix[3][3], maxWSoFar )
			})
			theMaxTranslationDistanceInRadians = safeAcosh(theMaxW)
		}

		//	We want theFogCutoff to be somewhat greater than
		//	theMaxTranslationDistanceInRadians, for two reasons:
		//
		//	1.	Roughly half of the translated cell will, of course,
		//		be more distant than its center.
		//
		//	2.	We want the most distant parts of the lattice
		//		to be at least dimly visible.
		//
		let theFogCutoff = 1.15 * theMaxTranslationDistanceInRadians
		
		return theFogCutoff
	}

	func makeCheckerboardPoseBuffer(
		device: MTLDevice,
		checkerboardPoseMatrices: [simd_double4x4]
	) -> MTLBuffer {	//	buffer of [simd_float4x4]

		let thePosesAsFloat = checkerboardPoseMatrices.map() { pose in
			convertDouble4x4toFloat4x4(pose)
		}
		
		guard let thePoseBuffer = device.makeBuffer(
				bytes: thePosesAsFloat,
				length: thePosesAsFloat.count *
							MemoryLayout<simd_float4x4>.stride,
				options: [.storageModeShared]
		) else {
			preconditionFailure("Couldn't make thePoseBuffer in makeCheckerboardPoseBuffer()")
		}
		
		return thePoseBuffer
	}

	// MARK: -
	// MARK: Spherical tiling

	private func sphericalCheckerboardPoses() -> [simd_double4x4] {

		let theRawPoses = [
			simd_diagonal_matrix(SIMD4<Double>(  1.0,  1.0,  1.0,  1.0 )),
			simd_diagonal_matrix(SIMD4<Double>( -1.0, -1.0,  1.0,  1.0 )),
			simd_diagonal_matrix(SIMD4<Double>( -1.0,  1.0, -1.0,  1.0 )),
			simd_diagonal_matrix(SIMD4<Double>( -1.0,  1.0,  1.0, -1.0 )),
			simd_diagonal_matrix(SIMD4<Double>(  1.0, -1.0, -1.0,  1.0 )),
			simd_diagonal_matrix(SIMD4<Double>(  1.0, -1.0,  1.0, -1.0 )),
			simd_diagonal_matrix(SIMD4<Double>(  1.0,  1.0, -1.0, -1.0 )),
			simd_diagonal_matrix(SIMD4<Double>( -1.0, -1.0, -1.0, -1.0 ))
		]

		//	The base cell is a right-angled tetrahedron
		//	with a corner at the north pole (0,0,0,1)
		//	and edges running in the x, y and z directions.
		//	First rotate it face the player in his/her default position...
		let theRotation = axisAlignedRotation(
							from: 0,
							towards: 2,
							angle: 0.25 * Double.pi)
		
		//	...and then slide it away from the player so that
		//	the tetrahedron's bottom face coincides with the pool table.
		//	To compute the necessary distance, we apply
		//	the Spherical Law of Cosines for Angles
		//
		//		cos γ = -cos α cos β + sin α sin β cos c
		//
		//	to a triangle with angles π/2, π/3 and π/4, giving
		//
		//		cos(π/2) = -cos(π/3) cos(π/4) + sin(π/3) sin(π/4) cos c
		//	 =>
		//		0 = -√¼ √½ + √¾ √½ cos c
		//	 =>
		//		0 = -1 + √3 cos c
		//	 =>
		//		c = acos(1/√3)
		//
		let c = acos( 1.0 / sqrt(3.0) )
		let theTranslation = axisAlignedTranslation(
								axis: 2,
								distanceInRadians: -c,
								geometryType: .spherical)
		
		let theAdjustment = theTranslation * theRotation	//	right-to-left matrix action
		
		let thePoses = theRawPoses.map() { rawPose in
			theAdjustment * rawPose	//	right-to-left matrix action
		}
		
		return thePoses
	}


	// MARK: -
	// MARK: Euclidean tiling

	private func euclideanCheckerboardPoses() -> [simd_double4x4] {

		let n = 8
		
		let theEdgeLength = GeometryType.euclidean.tableEdgeLengthInRadians

		//	Adjustment (in radians) to raise the central cell's
		//	bottom face to table height.
		let Δy = 0.5 * theEdgeLength
		
		var thePose = matrix_identity_double4x4
		
		var thePoses: [simd_double4x4] = []
		
		//	Let the identity matrix be the first one in the array,
		//	so that the renderer will know to omit its bottom face.
		thePose[3][0] = 0.0			//	redundant but clear
		thePose[3][1] = 0.0 + Δy
		thePose[3][2] = 0.0			//	redundant but clear
		thePoses.append(thePose)
		
		for x in -n ... +n {
			for y in -n ... +n {
				for z in -n ... +n {
				
					//	Omit alternate cells.
					//
					//		Note: '%' is a 'remainder' operator
					//		not a 'modulus' operator -- it handles
					//		negative numbers differently.
					//
					if abs(x + y + z) % 2 == 1 {
						continue
					}
					
					//	Omit cells that are too far away.
					if x*x + y*y + z*z > n*n {
						continue
					}
					
					//	We've already put the identity pose onto the array.
					if x == 0 && y == 0 && z == 0 {
						continue
					}
					
					thePose[3][0] = theEdgeLength * Double(x)
					thePose[3][1] = theEdgeLength * Double(y) + Δy
					thePose[3][2] = theEdgeLength * Double(z)
					
					thePoses.append(thePose)
				}
			}
		}
		
		return thePoses
	}


	// MARK: -
	// MARK: Hyperbolic tiling

	private func hyperbolicCheckerboardPoses() -> [simd_double4x4] {

		//	Conceptually our starting point is the group generated
		//	by reflections across the 12 faces of a right-angled
		//	hyperbolic dodecahedron, but we exclude the mirror-reversed
		//	elements and keep only non-reversed ones. The renderer
		//	will draw the non-reversed copies of the dodecahedron
		//	with face culling disabled, so in effect the "back faces"
		//	of the non-reversed dodecahedra provide the "front faces"
		//	for the omitted mirror-reversed dodecahedra.
		//	By using only the non-reversed group elements,
		//	our tiling algorithm should run more quickly,
		//	and then at render time we'll be requesting only
		//	half as many instances of the dodecahedral mesh
		//	as we would need if we were drawing mirror reversed
		//	instances as well.
		//
		//	In practical terms, each generator for our group will be
		//	a 180° rotation about one of the right-angled dodecahedron's edges.
		//	The dodecahedron has 30 edges, so we'll have 30 generators.
		//
		let theGenerators = hyperbolicGenerators()
		
		//	Computationally each (non-reversed!) group element is a pose,
		//	that is, a 4×4 matrix in SO(3,1), and nothing more.
		//	But conceptually it's convenient to visualize each pose
		//	as a copy of the dodecahedron. Neighboring poses meet
		//	along the shared edge of the two dodecahedron.
		//
		//	It's easy to generate all such poses by multiplying generators
		//	together. The challenging part is to make sure we produce
		//	each pose once and only once. As we tile outwards, we think
		//	of each new pose as the "child" of the "parent" pose
		//	that produced it. That is
		//
		//		child = parent · generator		(right-to-left matrix action)
		//
		//	To avoid duplications, each child will have a "preferred parent":
		//	the child will let itself get generated by its preferred parent,
		//	but will reject attempts from other potential parents to generate it.
		//	Each pose has 30 potential parents, one for each of the dodecahedron's edges.
		//
		//		Definition. Among a pose's 30 potential parents,
		//		the preferred parent is the one that's closest to the origin,
		//		that is, the one that translates the origin by the least amount.
		//		If more than one potential parent translates the origin
		//		by the same amount, the preferred parent is the one
		//		that appears first in the array theGenerators.
		//
		//	The only tricky part will be deciding whether one distance
		//	is less than, equal to, or greater than another.
		//	We can set some tolerance for deciding equality,
		//	but we'll need to keep in mind that if we were to ever
		//	tile out to progressively larger distances, the matrix entries
		//	in the O(3,1) matrices would grow exponentially fast,
		//	and the numerical errors would also grow exponentially fast.
		//	For details, see
		//
		//		The Achilles' heel of O(3,1)?,
		//		Experimental Mathematics 11 (2002) 91-97
		//
		//	So if we ever want to tile substantially further,
		//	we might need to reconsider our tolerances.

		//	Start with a single "seed" pose alone on thePoseList.
		//	The seed pose aligns a face of the dodecahedron
		//	with the pool table.
		let theSeedPose = hyperbolicSeedPose()
		var thePoseList = [theSeedPose]
		var theTotalPoseCount = 1
		
		//	Our algorithm will be to walk along thePoseList,
		//	examining each pose in turn and adding each of its children
		//	to thePoseList if and only if the given pose is the child's
		//	preferred parent and the depth hasn't exceeded theMaxDepth.
		//
		//		Performance note: Setting theMaxDepth = 5.3
		//		provide aesthetically pleasing results,
		//		no obvious flicker once fog is enabled,
		//		and it uses only ~50% of available GPU time
		//		-- no noticeable overheating during my modest
		//		testing.
		//
		let theMaxDepth = 5.3	//	in radians
		let theMaxW = cosh(theMaxDepth)
		
		precondition(
			theGenerators.count == 30,
			"Internal error: we expected 30 generators, one for each edge of the dodecahedron.")

		var theProcessedPoseCount = 0
		while theProcessedPoseCount < theTotalPoseCount {

			let theParent = thePoseList[theProcessedPoseCount]
			
			for i in 0...29 {
			
				let theGenerator = theGenerators[i]
				let theChild = theParent * theGenerator	//	right-to-left matrix action
				
				//	Does this child lie beyond our maximum allowed depth?
				if theChild[3][3] > theMaxW {
					continue
				}
				
				//	Are we not this child's preferred parent?
				if preferredParentIndex(for: theChild, generators: theGenerators) != i {
					continue
				}
				
				//	Add this child to thePoseList.
				thePoseList.append(theChild)
				theTotalPoseCount += 1
			}
			
			theProcessedPoseCount += 1
		}
		
		return thePoseList
	}

	private func hyperbolicGenerators() -> [simd_double4x4] {

		//	Recall from hyperbolicCheckerboardPoses() above that
		//	each generator for our group will be a 180° rotation
		//	about one of the right-angled dodecahedron's edges.
		//	We'll construct each such 180° rotation as the product
		//	of the reflections in the two faces of the dodecahedron
		//	adjacent to the given edge.

		//	To keep the coordinate system simple, it's convenient
		//	to construct a dodecahedron by starting with a cube and then
		//	put a little "roof" onto each of the cube's six faces, like this
		//
		//		°°°°°°°°°°°°°°°°°
		//		°\             /°
		//		° \           / °
		//		°  \         /  °
		//		°   ---------   °
		//		°  /         \  °
		//		° /           \ °
		//		°/             \°
		//		°°°°°°°°°°°°°°°°°
		//
		//	By adjusting the ridgeline's length and its height above
		//	the face of the cube, each trapezoidal facet of one roof
		//	will align with a triangular facet of an adjacent roof
		//	to give a regular pentagon. The underlying cube has 12 edges,
		//	so we get 12 such pentagonal faces, which together form
		//	a regular dodecahedron.

		//	When writing the Curved Spaces app, I worked out an exact matrix
		//	for the reflection in each of the right-angled dodecahedron's 12 faces,
		//	in terms of the golden ratio φ ≈ 1.61803… which is a root
		//	of the quadratic polynomial
		//
		//		φ² - φ - 1 = 0
		//

		//	To get started, let's compute the matrix for a 180° rotation
		//	about the ridgeline on the "positive x" face of the underlying cube.
		//	To get this matrix, we may compose the reflections across
		//	the two adjacent pentagonal faces, in either order (they commute):

		//	Matrix A
		//		(    -φ      φ   0 -√(2φ + 1) )		//	left-to-right action, row-major format
		//		(     φ      0   0     √φ     )
		//		(     0      0   1      0     )
		//		( √(2φ + 1) -√φ  0     φ+1    )
		//
		let matrixA = simd_double4x4(
			SIMD4<Double>( -1.61803398874989484820,  1.61803398874989484820,  0.00000000000000000000, -2.05817102727149225032 ),
			SIMD4<Double>(  1.61803398874989484820,  0.00000000000000000000,  0.00000000000000000000,  1.27201964951406896425 ),
			SIMD4<Double>(  0.00000000000000000000,  0.00000000000000000000,  1.00000000000000000000,  0.00000000000000000000 ),
			SIMD4<Double>(  2.05817102727149225032, -1.27201964951406896425,  0.00000000000000000000,  2.61803398874989484820 )
		)
		
		//	Matrix B
		//		(    -φ     -φ   0 -√(2φ + 1) )	//	left-to-right action, row-major format
		//		(    -φ      0   0     -√φ    )
		//		(     0      0   1      0     )
		//		( √(2φ + 1)  √φ  0     φ+1    )
		//
		let matrixB = simd_double4x4(
			SIMD4<Double>( -1.61803398874989484820, -1.61803398874989484820,  0.00000000000000000000, -2.05817102727149225032 ),
			SIMD4<Double>( -1.61803398874989484820,  0.00000000000000000000,  0.00000000000000000000, -1.27201964951406896425 ),
			SIMD4<Double>(  0.00000000000000000000,  0.00000000000000000000,  1.00000000000000000000,  0.00000000000000000000 ),
			SIMD4<Double>(  2.05817102727149225032,  1.27201964951406896425,  0.00000000000000000000,  2.61803398874989484820 )
		)
		
		let ridgelineRotation = matrixA * matrixB
		
		//	Next let's compute the matrix for a 180° rotation
		//	about one of the "side ridges" of the roof
		//	on the "positive x" face of the underlying cube.
		//	It will be the product of matrixB and a new matrixC that
		//	represents a reflection across one of the "end triangles"
		//	of that roof. Again the order in which we multiply
		//	the two matrices doesn't matter, because they commute.
		//
		//	Matrix C
		//		(  0  0     -φ        -√φ    )	//	left-to-right action, row-major format
		//		(  0  1      0         0     )
		//		( -φ  0     -φ    -√(2φ + 1) )
		//		( √φ  0  √(2φ + 1)    φ+1    )
		//
		let matrixC = simd_double4x4(
			SIMD4<Double>(  0.00000000000000000000,  0.00000000000000000000, -1.61803398874989484820, -1.27201964951406896425 ),
			SIMD4<Double>(  0.00000000000000000000,  1.00000000000000000000,  0.00000000000000000000,  0.00000000000000000000 ),
			SIMD4<Double>( -1.61803398874989484820,  0.00000000000000000000, -1.61803398874989484820, -2.05817102727149225032 ),
			SIMD4<Double>(  1.27201964951406896425,  0.00000000000000000000,  2.05817102727149225032,  2.61803398874989484820 )
		)
		
		let theSideRotation00 = matrixB * matrixC
		
		//	The dodecahedron's symmetry gives us the 180° rotations
		//	about the three remaining "side ridges" of the roof
		//	on the "positive x" face of the underlying cube.
		//	We just need to conjugate theSideRotation00
		//	by a reflection in y and/or z direction.

		let theYReflection = simd_diagonal_matrix(SIMD4<Double>( 1.0, -1.0,  1.0, 1.0))
		let theSideRotation10 = theYReflection * theSideRotation00 * theYReflection

		let theZReflection = simd_diagonal_matrix(SIMD4<Double>( 1.0,  1.0, -1.0, 1.0))
		let theSideRotation01 = theZReflection * theSideRotation00 * theZReflection

		let theYZRotation  = simd_diagonal_matrix(SIMD4<Double>( 1.0, -1.0, -1.0, 1.0))
		let theSideRotation11 = theYZRotation  * theSideRotation00 * theYZRotation
		
		//	We've now got a group generator (i.e. a 180° rotation)
		//	for each of the five edges belonging to the roof
		//	on the "positive x" face of the underlying cube.
		let theGeneratorsXPos = [
			ridgelineRotation,
			theSideRotation00,
			theSideRotation01,
			theSideRotation10,
			theSideRotation11
		]
		
		//	To get the group generators for the cube's
		//	other five faces, we can take theGeneratorsXPos
		//	and conjugate them by an appropriate rotation
		//	of the cube.
		
		let theXYZRotation = simd_double4x4(
			SIMD4<Double>( 0.0, 1.0, 0.0, 0.0 ),
			SIMD4<Double>( 0.0, 0.0, 1.0, 0.0 ),
			SIMD4<Double>( 1.0, 0.0, 0.0, 0.0 ),
			SIMD4<Double>( 0.0, 0.0, 0.0, 1.0 )
		)
		let theZYXRotation = simd_double4x4(
			SIMD4<Double>( 0.0, 0.0, 1.0, 0.0 ),
			SIMD4<Double>( 1.0, 0.0, 0.0, 0.0 ),
			SIMD4<Double>( 0.0, 1.0, 0.0, 0.0 ),
			SIMD4<Double>( 0.0, 0.0, 0.0, 1.0 )
		)
		
		let theGeneratorsYPos = theGeneratorsXPos.map() { generator in
			theXYZRotation * generator * theZYXRotation	//	right-to-left matrix action
		}
		let theGeneratorsZPos = theGeneratorsXPos.map() { generator in
			theZYXRotation * generator * theXYZRotation	//	right-to-left matrix action
		}
		
		let theXReflection = simd_diagonal_matrix(SIMD4<Double>( -1.0,  1.0,  1.0,  1.0 ))
		//	theYReflection is already defined above
		//	theZReflection is already defined above

		let theGeneratorsXNeg = theGeneratorsXPos.map() { generator in
			theXReflection * generator * theXReflection
		}
		let theGeneratorsYNeg = theGeneratorsYPos.map() { generator in
			theYReflection * generator * theYReflection
		}
		let theGeneratorsZNeg = theGeneratorsZPos.map() { generator in
			theZReflection * generator * theZReflection
		}
		
		let theGenerators
			= theGeneratorsXPos
			+ theGeneratorsYPos
			+ theGeneratorsZPos
			+ theGeneratorsXNeg
			+ theGeneratorsYNeg
			+ theGeneratorsZNeg

		return theGenerators
	}

	private func hyperbolicSeedPose() -> simd_double4x4 {

		let π = Double.pi

		//	The function hyperbolicGenerators() above uses
		//	a cube-with-roofs construction of the dodecahedron,
		//	which places one of the dodecahedron's edges
		//	along the bottom. In order for our seed pose
		//	to align with the pool table, we need to rotate
		//	the dodecahedron to place a face at the bottom instead.
		//	To do this, we rotate the z axis towards the y axis
		//	through an angle φ equal to the angular distance
		//	from the center of desired face to the midpoint
		//	of the edge that's currently at the bottom.
		//	This angle is easily found using
		//	the Spherical Law of Cosines for angles
		//
		//		cos γ = -cos α cos β + sin α sin β cos φ
		//
		//	which in the present case gives
		//
		//		cos(π/3) = -cos(π/5) cos(π/2) + sin(π/5) sin(π/2) cos(φ)
		//	=>
		//		1/2 = sin(π/5) cos(φ)
		//	=>
		//		φ = acos( 1/2 / sin(π/5) )
		//
		let φ = acos( 0.5 / sin( 0.2 * Double.pi ) )
		let theRotation = axisAlignedRotation(
							from: 2,
							towards: 1,
							angle: φ)
		
		//	We now need to raise the dodecahedron through a distance
		//	(in radians) equal to the dodecahedron's inradius, which
		//	we may find using the Hyperbolic Law of Cosines for Angles
		//	as explained in the comments at the beginning
		//	of gRightAngledDodecahedron in NEBCheckerboardCells.swift.
		//	Those comments also provide the values for the angles
		//	θ₀, θ₂ and θ₃:

		let cosθ₀ = sqrt(2.0/3.0)
	//	let sinθ₀ = sqrt(1.0/3.0)

	//	let cosθ₂ = 0.0	//	cos(π/2)
	//	let sinθ₂ = 1.0	//	sin(π/2)

		let cosθ₃ = ( 0.5 * cos(π/5) ) / ( sqrt(0.75) * sin(π/5) )
		let sinθ₃ = sqrt( 1.0 - cosθ₃ * cosθ₃ )

		//	We may now apply the Hyperbolic Law of Cosines for Angles
		//
		//		cos γ = -cos α cos β + sin α sin β cosh c
		//
		//	but with a different assignment of the angles
		//	θ₀, θ₂ and θ₃ to the variables α, β and γ than
		//	we use in gRightAngledDodecahedron, because here
		//	we want to solve for the inradius r instead of the outradius.
		//
		//		cos(θ₀) = -cos(θ₂) cos(θ₃) + sin(θ₂) sin(θ₃) cosh(r)
		//	=>
		//		cos(θ₀) = sin(θ₃) cosh(r)
		//	=>
		//		cosh(r) = cos(θ₀) / sin(θ₃)
		//
		let r = acosh( cosθ₀ / sinθ₃ )
		let theTranslation = axisAlignedTranslation(
								axis: 1,
								distanceInRadians: r,
								geometryType: .hyperbolic)

		let theSeedPose = theTranslation * theRotation	//	right-to-left matrix action
		
		return theSeedPose
	}

	private func preferredParentIndex(
		for child: simd_double4x4,
		generators: [simd_double4x4]
	) -> Int? {

		let theChildW = child[3][3]
		
		//	thePreferredIndex is the index of the generator
		//	that connects the child to its preferred parent.
		//
		//	Exception: The seed pose, which sits close to
		//	the north pole (0,0,0,1), doesn't have a preferred parent.
		//	Indeed, because hyperbolicCheckerboardPoses() uses it
		//	to seed thePoseList, it doesn't need a parent.
		//	Morever, it's also exceptional in that it's the only pose
		//	all of whose neighboring poses have larger w values
		//	than it does. So when the caller passes us the seed pose,
		//	we return nil.
		//
		var thePreferredIndex: Int? = nil
		var theSmallestW = theChildW

		for i in 0...29 {
		
			let theGenerator = generators[i]
		
			//	A full matrix multiplication would give
			//
			//		theParent = child * theGenerator	(right-to-left matrix action)
			//
			//	but really all we need to the theParent[3][3],
			//	which we can compute with just a dot product
			//	(no need for the full matrix multiplication).
			//
			let thePotentialParentW
				= dot(
					SIMD4<Double>(child[0][3], child[1][3], child[2][3], child[3][3]),
					theGenerator[3] )
		
			//	Recall from a comment in hyperbolicCheckerboardPoses()
			//	that the definition of the "preferred parent" requires
			//	that in the case of two potential parents with the same w value,
			//	we must prefer the one that comes first in the list
			//	of generators. To satisfy that requirement, we accept
			//	a potential parent iff its w value is at least ε less
			//	than our previous best candidate.
			//
			//	Picking a suitable value for ε is very easy. Even when
			//	tiling out to theMaxDepth = 6.0 (which is as deep as
			//	we'd ever want to go), the largest error between
			//	two w values that are judged to be equal is ~1.83e-12,
			//	while the smallest error between two w values that are
			//	judged to be unequal is ~0.74. So we can set ε to any value
			//	that's well away from those two extremes.
			//
			//		Note: The reason that the 0.74 value is as large
			//		as it is is that for poses near the tiling's center,
			//		neighboring tiles sit at either exactly the same
			//		distance from the center or at substantially different
			//		distances. More distant poses might sit at almost
			//		the same distance, but once we're far from the center
			//		the w coordinate increases exponentially fast
			//		as a function of the distance from the center,
			//		so small differences in the distance would yield
			//		large differences in w.
			//
			let ε = 1.0e-6
			
			if thePotentialParentW < theSmallestW - ε {
				theSmallestW = thePotentialParentW
				thePreferredIndex = i
			}
		}
		
		return thePreferredIndex
	}
}
