Euler Problem 75 asks:
It turns out that 12 cm is the smallest length of wire that can be bent to form an integer sided right angle triangle in exactly one way, but there are many more examples.
12 cm: (3,4,5)
24 cm: (6,8,10)
30 cm: (5,12,13)
36 cm: (9,12,15)
40 cm: (8,15,17)
48 cm: (12,16,20)In contrast, some lengths of wire, like 20 cm, cannot be bent to form an integer sided right angle triangle, and other lengths allow more than one solution to be found; for example, using 120 cm it is possible to form exactly three different integer sided right angle triangles.
120 cm: (30,40,50), (20,48,52), (24,45,51)
Given that L is the length of the wire, for how many values of L <= 2,000,000 can exactly one integer sided right angle triangle be formed?
Everything discussed here can be found at http://en.wikipedia.org/wiki/Pythagorean_triple
Integer Pythagorean Triples (keep [3,4,5] in mind–it’s the first example of one) are values that form right triangles. Triples that are not multiples of other triples are called primitive. [3,4,5] is the first integer triple, and it’s obviously primitive. [6,8,10] is a triple, but as a 2x multiple of [3,4,5], it is not primitive. Similarly for [9,12,15].
Using Euclid’s formula (upcoming), integer primitive triples can be computed from a pair of numbers m and n with the following properties:
- m greater than or equal to 2
- n less than m
- n and m of opposite parity
- m and n are relatively prime (GCD = 1)
With this understanding then, Euclid found Side a is m^2 – n^2. Side b is 2*m*n and Side c is m^2 + n^2. For [3,4,5], m = 2 and n = 1.
- 3 = 2^2 – 1^2
- 4 = 2*2*1
- 5 = 2^2 + 1^2
If you add (m^2 + n^2) + (2*m*n) + (m^2 – n^2) to get the perimeter, it simplifies to 2*m*(m + n).
The [3,4,5] triangle’s perimeter is
- 3 + 4 + 5 = 12 = 2*2*(2 +1)
and 12 is the minimum value for any integer triple.
Some things fall out from the above properties.
- With a value of 2*m*(m + n) , all right triangle perimeters are even.
- For a primitive triangle, one side is even, and the other two are odd (to add to an even)
- For a primitive triangle, Side c is always odd
- For multiples, either one side is even, with the other two are odd, or they are all even.
- And since area = 1/2 * a * b, and b is even, the area of an integer right triangle is always a whole number.
Euler asks us for perimeters less or equal to 2,000,000. We can use that to find some bounds on c and m. For all triangles, one way is: a + b > c. Then a + b + c > c + c. 2000000>= a + b +c means 2000000 > 2*c, and 1000000 > c. Take c = m^2 +n^2 = 1000000, then for n = 1, m^2 + 1 = 1000000 and m is approximately the square root of 10^6, or 10^3.
Alternatively, if 2*m*(m + n) = 2000000, then m*(m + n) = 1000000. For n = 1, m^2 + 1 = 1000000, and we are at the same place.
We can build loops with m = 2 to Sqr(1000000) and n = 1 + (m mod 2) to m – 1 step 2 (for opposite parity) to build all primitives. If we multiply the primitive perimeter P by an incrementing K until K*P is greater than 2,000,000, and increment the appropriate counter L(K*P) by 1 for each case, we have accounted for and binned all the perimeters. Loop through L() from 12 to 2000000 Step 2, and tally the L(k) equal to 1.
Here is the code that does all that, with the GCD2 function. It runs in about eight-tenths of a second, and also reports out the number of primitives (PPT).
Const Max_C As Currency = 1000000
Const Max_P As Currency = 2000000
Dim a As Currency
Dim b As Currency
Dim c As Currency
Dim P As Currency
Dim m As Long
Dim n As Long
Dim K As Long
Dim PPT As Long ‘Primitive Pythagorean Triple
Dim L(Max_P) As Currency
Dim Answer As Long, T As Single
‘ Dim Squared(Max_C) As Currency
T = Timer
‘ For a = 1 To Max_C
‘ Squared(a) = a * a
‘ Next a
For m = 2 To Sqr(Max_C)
For n = 1 + m Mod 2 To m – 1 Step 2
If GCD2(m, n) != 1 Then GoTo Next_n
a = m ^ 2 – n ^ 2
b = 2 * m * n
c = m ^ 2 + n ^ 2
If c GT Max_C Then Exit For
‘ If Squared(a) + Squared(b) = Squared(c) Then
P = a + b + c ‘ also = 2*m*(m+n)
PPT = PPT + 1
K = 1
If P LTE Max_P then
Do
L(K * P) = L(K * P) + 1
K = K + 1
Loop Until K * P GT Max_P
End If
‘ End If
Next_n:
Next n
Next m
For K = 12 To Max_P Step 2 ‘ all perimeters are even
If L(K) = 1 Then
Answer = Answer + 1
End If
Next K
Debug.Print Answer; ” Time:”; Timer – T; PPT
End Sub
Function GCD2(ByVal a As Long, ByVal b As Long)
Dim T As Long
While b != 0
T = b
b = a Mod b
a = T
Wend
GCD2 = a
End Function
When I first wrote the code I included a test (commented out above) for “right-triangle-ness” using precomputed squares up to 1000000^2. This test used currency variables to hold 10^12. The final loop used a currency-type counter to total cases of L(P) = 1. That loop never went above 1717986, and the wrong answer from a right method resulted. That was fixed by using k, a long. See the Code Snipit thread. That right-angled test was never failed, or Euclid’s formula would have failed, and yours truly would be in wikipedia with a dissertation topic ;-). It’d only take one example for fame and glory, but it’s just not needed. But then I wouldn’t have found the bug. Some irony, I suppose.
All the currency types can be changed to longs, and if you use the alternative calculation for P, a and b need not be computed. Though this isn’t how I did them, this also applies to Problems 9 and 39. There are the usual angle bracket adjustments.
…mrt
Iteration is again the key.
It can be shown that for a .LT. b .LT. c, c – b = 1 or 2 only for primitive triples. When c – b = 1, this means the difference between sequential squares b2 and c2 must also be a square. Since n2 is the sum of the first n positive odd integers, the square of every odd number is in a primitive triple, so
a = 3, b = 4, c = 5, a + b + c = 12 = 3 * 4 = a(a + 1)
a = 5, b = c – 1 = 12, c = (a2 + 1)/2 = 13, a + b + c = 30 = 5 * 6 = a(a + 1)
a = 7, b = c – 1 = 24, c = (a2 + 1)/2 = 25, a + b + c = 56 = 7 * 8 = a(a + 1)
a = 9, b = c – 1 = 40, c = (a2 + 1)/2 = 41, a + b + c = 90 = 9 * 10 = a(a + 1)
So sufficient to iterate a = 3 to a = 1413 step 2. For c – b = 2, a can be every multiple of 4 starting with 8 (a = 4 corresponds to b = 3, but that violates a .LT. b).
a = 8, b = 15, c = 17, a + b + c = 40 = 8 * 5 = a(a/2 + 1)
a = 12, b = 35, c = 37, a + b + c = 84 = 12 * 7 = a(a/2 + 1)
a = 16, b = 63, c = 65, a + b + c = 144 = 16 * 9 = a(a/2 + 1)
So sufficient to iterare a = 8 to 1996 step 4.
The a(a + 1) and a(a/2 + 1) terms give the perimeters of primitive triples. Their multiples are perimeters of nonprimitive triples. Easiest to iterate from 1 to 2E6 primitive perimeter storing the multiples of the primitive perimeters as keys in a WSH Dictionary object with the count as the key’s values. Then get the result by counting the resulting .Item values equal to 1.
Hi fzz –
Different ways to get to P. If one uses P = 2m(m+n), then you don’t need a and b. I used c only as a kick-out to shave some msecs that I might waste in calculating it. Taking m to 999 and n to 2 does give c > 1000000.
Significantly less loops your way, though. I put in a counter and I reach Next n an even 250000 times.
My only experience with dictionaries, actually, is with this problem, and I found stuffing a series/array faster.
What do you get for run time?
…mrt
I also worked the problem using the T1/2/3 transforms mentioned further down the Wikipedia page, which allow recursion to be applied. I try to take any opportunity to code a recursive function, since I still have to think them through very carefully. That solution ran in about half the time of the Euclid-based one.
I can’t believe I’d missed the stupid GCD function, by the way, although I didn’t use it for this problem: I just used the formulae to crank out primitives until they reached the maximum, and added the multiples with a loop.
I don’t think dictionaries have much of an advantage – if any at all – over arrays in this case. They’re very handy as a substitute for sparsely-populated arrays (vectors, really) of large size. When, as in this case, over 10% of the possible values are of interest (and many more are used but not really interesting) then the additional overhead of the hashing in the Dictionary lookup probably isn’t worth it.
Comparing Dictionaries against Collections:
On the one hand, you lose “For Each” when you move from a Collection.
On the other hand, you gain a lot of performance, plus the “.Exist” method.
On the gripping hand they both suck against Arrays when you load them up with non-trivial quantities of non-primitive objects and then destroy them.
Hi Mike –
As a recursive kind of guy ;-) you no doubt have this, but for those who don’t, here’s the recursive GCD version, GCD1().
If b = 0 Then
GCD1 = a
Else
GCD1 = GCD1(b, a Mod b)
End If
End Function
I’ve found GCD2() above to be faster when I swap ’em out.
…mrt
Michael, with such a tight routine, all the variant/long conversions will probably have quite an impact. If you declare a and b ‘As Long’, you should see it improve.
Good morning, Stephen –
Thank you. I left them as currency here to tie in with the Snipit thread. Exactly as you say, they are very properly “longs” since I no longer had a need to store large square values.
Have you seen the “Code Snipit to test” thread? That came about because of a need I didn’t really have in this problem. When Dick allowed me to author, I didn’t think I’d become an insectologist.
I also didn’t quite realize whose notice I’d attract.
Thank you.
…mrt