Nice Vectors
I wanted to be nice this semester in computer graphics. In my exercises, I ask students to do a lot of 3D vector math. These vectors must often be normalized, which means their magnitude is 1:
If that's the case, so is their squared magnitude:
I needed vectors for which this equation is true, but I didn't want to have a lot of digits in \(a\), \(b\), and \(c\). The more digits, the greater the chance for transcription and roundoff errors. If we lose some of those digits, the vectors lose their normality.
One way to generate nice vectors is to guess short values for \(a\), \(b\), and \(c\) and see if they satisfy the magnitude equation. But this is scattershot. Suppose I choose \(a = 0.5\) and \(b = 0.6\). Then I end up with \(c \approx 0.6244997998398398\), which is not nice at all.
I decided instead to iterate through all triplets of integers of a small magnitude, plug them into the equation, and keep the good ones. For instance, all single-digit triplets would have \(1 \leq a, b, c \leq 9\). These integers need to be normalized when plugging them in:
We can manipulate the equation to stick with integers throughout. For single-digit numbers, this is the equation we need to satisfy:
As it turns out, there are no integer triplets that satisfy this equation. For two-digit numbers, this is the equation we need to satisfy:
Exactly two vectors satisfy this equation. In normalized form, they are \([0.36, 0.48, 0.8]\) and \([0.48, 0.6, 0.64]\). Those are nice and tidy. But there are only two of them. To get more, I needed some code to generate all possible triplets with an arbitrary number of digits. This just felt like a job for Haskell. Here was my first draft solution:
magnitudeTriplets :: Int -> [[Int]]
magnitudeTriplets digitCount = roots
where
-- For 1 digit, the divisor is 10. For 2, 100.
-- For 3, 1000. But we really need its square.
divisor = 10 ^ digitCount
divisor2 = divisor ^ 2
-- What values can a^2, b^2, and c^2 take on?
-- Any perfect square less than divisor2.
squares = takeWhile (< divisor2) perfectSquares
-- Let's find all possible 3-combinations of
-- these squares and filter out the ones that
-- sum to divisor2.
triplets = combos 3 squares
matches = sumTriplets divisor2 triplets
-- Turn each [a^2, b^2, c^2] into [a, b, c].
roots = map (map sqrtFloor) matches
magnitudeTriplets :: Int -> [[Int]]
magnitudeTriplets digitCount = roots
where
-- For 1 digit, the divisor is 10. For 2, 100.
-- For 3, 1000. But we really need its square.
divisor = 10 ^ digitCount
divisor2 = divisor ^ 2
-- What values can a^2, b^2, and c^2 take on?
-- Any perfect square less than divisor2.
squares = takeWhile (< divisor2) perfectSquares
-- Let's find all possible 3-combinations of
-- these squares and filter out the ones that
-- sum to divisor2.
triplets = combos 3 squares
matches = sumTriplets divisor2 triplets
-- Turn each [a^2, b^2, c^2] into [a, b, c].
roots = map (map sqrtFloor) matchesThis function relies on the utilities below. Haskell probably has some module that can generate combinations, but I enjoy writing these.
-- [1, 4, 9, 16, 25, ...]
perfectSquares = map (^ 2) [1..]
-- Select all the triplets that sum to the target n.
-- For example:
-- sumTriplets 100 [(36, 48, 80), (30, 40, 50)] == [(36, 48, 80)]
sumTriplets n triplets = filter ((n ==) . sum) triplets
-- Find all the n-combinations of a list. Take each element
-- and prepend it on all the (n - 1)-combinations of its
-- successors. For example:
-- combos 1 [1, 2, 3] = [[1], [2], [3]]
-- combos 3 [1, 2, 3, 4] == [[1, 2, 3], [1, 2, 4], [1, 3, 4], [2, 3, 4]]
combos 0 _ = []
combos 1 xs = map singleton xs
combos _ [] = []
combos n (x : xs) = map (x :) (combos (n - 1) xs) ++ combos n xs
-- Compute the integral floor of the square root of
-- a number. For example:
-- sqrtFloor 25 == 5
-- sqrtFloor 103 == 10
sqrtFloor :: Int -> Int
sqrtFloor = floor . sqrt . fromIntegral
-- [1, 4, 9, 16, 25, ...] perfectSquares = map (^ 2) [1..] -- Select all the triplets that sum to the target n. -- For example: -- sumTriplets 100 [(36, 48, 80), (30, 40, 50)] == [(36, 48, 80)] sumTriplets n triplets = filter ((n ==) . sum) triplets -- Find all the n-combinations of a list. Take each element -- and prepend it on all the (n - 1)-combinations of its -- successors. For example: -- combos 1 [1, 2, 3] = [[1], [2], [3]] -- combos 3 [1, 2, 3, 4] == [[1, 2, 3], [1, 2, 4], [1, 3, 4], [2, 3, 4]] combos 0 _ = [] combos 1 xs = map singleton xs combos _ [] = [] combos n (x : xs) = map (x :) (combos (n - 1) xs) ++ combos n xs -- Compute the integral floor of the square root of -- a number. For example: -- sqrtFloor 25 == 5 -- sqrtFloor 103 == 10 sqrtFloor :: Int -> Int sqrtFloor = floor . sqrt . fromIntegral
Unfortunately, this strategy is really slow. We are generating all possible triplets, making this an \(O(n^3)\) algorithm, where \(n\) is the number of possible values of each component. With 3 digits, \(n = 1000\) and \(n^3 = 1,000,000,000\).
A faster stategy is to generate all possible combinations of just \(a\) and \(b\). That's an \(O(n^2)\) algorithm. For each combination, we calculate \(c\) to satisfy the magnitude equation. If \(c\) is an integer, then we found a triplet. This solution runs much faster:
magnitudeTriplets :: Int -> [[Int]]
magnitudeTriplets digitCount = foldr mix [] pairs
where
divisor = 10 ^ digitCount
divisor2 = divisor * divisor
squares = takeWhile (< divisor2) perfectSquares
pairs = combos 2 squares
mix ([a, b]) accum
| c > b && isSquare c = [sqrtFloor a, sqrtFloor b, sqrtFloor c] : accum
| otherwise = accum
where
c = divisor2 - a - b
magnitudeTriplets :: Int -> [[Int]]
magnitudeTriplets digitCount = foldr mix [] pairs
where
divisor = 10 ^ digitCount
divisor2 = divisor * divisor
squares = takeWhile (< divisor2) perfectSquares
pairs = combos 2 squares
mix ([a, b]) accum
| c > b && isSquare c = [sqrtFloor a, sqrtFloor b, sqrtFloor c] : accum
| otherwise = accum
where
c = divisor2 - a - bOur combination generator yields combinations that are in sorted order, like \([48, 60, 64]\). Since c is derived separately, we add in the c > b check. Without it, we get vectors like \([36, 80, 48]\). This is a duplicate of \([36, 48, 80]\), which we will also generate. We don't need both.
We need one additional utility to check if \(c\) is a square:
-- Determine whether a number is a square. For example:
-- isSquare 36 == true
-- isSquare 42 == false
isSquare :: Int -> Bool
isSquare x = (sqrtFloor x) ^ 2 == x
-- Determine whether a number is a square. For example: -- isSquare 36 == true -- isSquare 42 == false isSquare :: Int -> Bool isSquare x = (sqrtFloor x) ^ 2 == x
This new function finds these 14 vectors with 3-digit components:
- \([0.024, 0.64, 0.768]\)
- \([0.096, 0.36, 0.928]\)
- \([0.096, 0.48, 0.872]\)
- \([0.152, 0.48, 0.864]\)
- \([0.168, 0.224, 0.96]\)
- \([0.168, 0.576, 0.8]\)
- \([0.192, 0.48, 0.856]\)
- \([0.192, 0.64, 0.744]\)
- \([0.224, 0.6, 0.768]\)
- \([0.28, 0.576, 0.768]\)
- \([0.352, 0.36, 0.864]\)
- \([0.36, 0.48, 0.8]\)
- \([0.424, 0.48, 0.768]\)
- \([0.48, 0.6, 0.64]\)
And these 76 vectors with 4-digit components:
- \([0.0192, 0.168, 0.9856]\)
- \([0.0192, 0.224, 0.9744]\)
- \([0.0240, 0.0512, 0.9984]\)
- \([0.0240, 0.2304, 0.9728]\)
- \([0.0240, 0.64, 0.768]\)
- \([0.0512, 0.6384, 0.768]\)
- \([0.0656, 0.192, 0.9792]\)
- \([0.0656, 0.4992, 0.864]\)
- \([0.0768, 0.424, 0.9024]\)
- \([0.0768, 0.4976, 0.864]\)
- \([0.0784, 0.2688, 0.96]\)
- \([0.0960, 0.1392, 0.9856]\)
- \([0.0960, 0.2688, 0.9584]\)
- \([0.0960, 0.36, 0.928]\)
- \([0.0960, 0.4096, 0.9072]\)
- \([0.0960, 0.48, 0.872]\)
- \([0.0960, 0.5264, 0.8448]\)
- \([0.1216, 0.48, 0.8688]\)
- \([0.1216, 0.6288, 0.768]\)
- \([0.1296, 0.192, 0.9728]\)
- \([0.1296, 0.6272, 0.768]\)
- \([0.1344, 0.152, 0.9792]\)
- \([0.1344, 0.4992, 0.856]\)
- \([0.1392, 0.3456, 0.928]\)
- \([0.152, 0.4032, 0.9024]\)
- \([0.152, 0.48, 0.864]\)
- \([0.1664, 0.4752, 0.864]\)
- \([0.168, 0.224, 0.96]\)
- \([0.168, 0.2944, 0.9408]\)
- \([0.168, 0.3968, 0.9024]\)
- \([0.168, 0.576, 0.8]\)
- \([0.168, 0.6336, 0.7552]\)
- \([0.1808, 0.6144, 0.768]\)
- \([0.192, 0.2112, 0.9584]\)
- \([0.192, 0.3968, 0.8976]\)
- \([0.192, 0.48, 0.856]\)
- \([0.192, 0.64, 0.744]\)
- \([0.2112, 0.2816, 0.936]\)
- \([0.2112, 0.4416, 0.872]\)
- \([0.2112, 0.5616, 0.8]\)
- \([0.224, 0.2544, 0.9408]\)
- \([0.224, 0.4416, 0.8688]\)
- \([0.224, 0.6, 0.768]\)
- \([0.224, 0.6672, 0.7104]\)
- \([0.2304, 0.2928, 0.928]\)
- \([0.2304, 0.352, 0.9072]\)
- \([0.2304, 0.6272, 0.744]\)
- \([0.2368, 0.36, 0.9024]\)
- \([0.2368, 0.576, 0.7824]\)
- \([0.2688, 0.28, 0.9216]\)
- \([0.2688, 0.4416, 0.856]\)
- \([0.28, 0.576, 0.768]\)
- \([0.2816, 0.6, 0.7488]\)
- \([0.2928, 0.4096, 0.864]\)
- \([0.2928, 0.5696, 0.768]\)
- \([0.2928, 0.64, 0.7104]\)
- \([0.3264, 0.3648, 0.872]\)
- \([0.3264, 0.424, 0.8448]\)
- \([0.3456, 0.5392, 0.768]\)
- \([0.3456, 0.576, 0.7408]\)
- \([0.352, 0.36, 0.864]\)
- \([0.352, 0.4752, 0.8064]\)
- \([0.352, 0.5616, 0.7488]\)
- \([0.36, 0.48, 0.8]\)
- \([0.36, 0.6336, 0.6848]\)
- \([0.3648, 0.5264, 0.768]\)
- \([0.3968, 0.48, 0.7824]\)
- \([0.3968, 0.5376, 0.744]\)
- \([0.4032, 0.4976, 0.768]\)
- \([0.424, 0.48, 0.768]\)
- \([0.4416, 0.6288, 0.64]\)
- \([0.4464, 0.48, 0.7552]\)
- \([0.4464, 0.576, 0.6848]\)
- \([0.48, 0.5696, 0.6672]\)
- \([0.48, 0.6, 0.64]\)
- \([0.5392, 0.576, 0.6144]\)
There are 388 normalized vectors with 5-digit components. I won't list them here.
The faster version of magnitudeTriplets takes 11 seconds on my MacBook Pro to find the vectors with 4-digit components. I gave up on the slower solution after five minutes.