Euler Problem 102 asks:
‘-1000 LTE x, y LTE 1000, such that a triangle is formed.
‘
‘Consider the following two triangles:
‘
‘A(-340,495), B(-153,-910), C(835,-947)
‘
‘X(-175,41), Y(-421,-714), Z(574,-645)
‘
‘It can be verified that triangle ABC contains the origin, whereas
‘triangle XYZ does not.
‘
‘Using triangles.txt (right click and ‘Save Link/Target As…’), a 27K text
‘file containing the co-ordinates of one thousand “random” triangles, find
‘the number of triangles for which the interior contains the origin.
‘
‘NOTE: The first two examples in the file represent the triangles in the
‘example given above.
In the above, LTE is “less than or equal”, to outsmart the html bugs.
Any triangle that has all-positive x-coordinates, or all-negative x-coordinates, or similarly all-positive or all-negative y-coordinates, can not contain the origin. A pre-screen will throw all those out. When I first solved this problem, I looked for triangles with both a positive and a negative y-intercept, coupled with both a positive and negative x-intercept, and counted those triangles. That was a successful strategy. LINEST() will return x-intercepts it you swap known-xs and known-ys in the formula. The code looked like this:
Application.WorksheetFunction.LinEst(Known_Ys, Known_Xs), 2)
XIntercept(3, 1) = Application.WorksheetFunction.Index( _
Application.WorksheetFunction.LinEst(Known_Xs, Known_Ys), 2)
It’s right out of the Help for LINEST(), or at least the y-intercept part is. However, after I checked my answer in I read another strategy that was just so flat-out neat that I coded it up. It uses Heron’s Law which calculates the area of a triangle based on a calculation of the semi-perimeter, or 1/2 the sum of the sides: The area of a triangle, given the lengths a,b,c of the sides, is A = sqrt s*(s-a)*(s-b)*(s-c), where s is the semiperimeter 0.5*(a+b+c). If we calculate the area of the suspect triangle, and then the area of the three sub-triangles with a common vertex of the origin, if the sum of the sub-triangle areas equals the area of the whole, then the origin must be within the suspect triangle. That code looked like this:
Option Base 1
Sub Problem_102B()
Dim Triangle(1000) As Variant
Dim i As Long
Dim j As Long
Dim Most As Long
Dim Impossible As Boolean
Dim IsTest As Boolean
Dim x(3) As Long ‘X Coordinates
Dim Y(3) As Long ‘Y Coordinates
Dim d(3) As Double ‘Distance between points
Dim O(3) As Double ‘Distance from points to (0,0)
Dim Area As Double ‘Triangle area
Dim a(3) As Double ‘Sub-triangle area
Dim T As Single
Dim Answer As Long
Const text As String = “D:DownloadsEuler riangles.txt”
T = Timer
IsTest = False
If IsTest Then
Most = 2
Else
Most = 1000
End If
i = 1
Open text For Input As #1 ‘1000 lines, comma delimited
Do While Not EOF(1)
Line Input #1, Triangle(i)
Triangle(i) = Split(Triangle(i), “,”) ‘zero-based array
i = i + 1
Loop
Close #1
For i = 1 To Most
Impossible = False
x(1) = Triangle(i)(0) ‘zero-based array
Y(1) = Triangle(i)(1)
x(2) = Triangle(i)(2)
Y(2) = Triangle(i)(3)
x(3) = Triangle(i)(4)
Y(3) = Triangle(i)(5)
‘For Triangles all above or all below the X-axis
If Y(1) > 0 And Y(2) > 0 And Y(3) > 0 Then Impossible = True
If Y(1) < 0 And Y(2) < 0 And Y(3) < 0 Then Impossible = True
‘For Triangles all to left of or all to right of the Y-axis
If Not Impossible Then ‘yet
If x(1) > 0 And x(2) > 0 And x(3) > 0 Then Impossible = True
If x(1) < 0 And x(2) < 0 And x(3) < 0 Then Impossible = True
End If
If Impossible Then GoTo Next_i
d(1) = TriSides(CDbl(x(1) – x(2)), CDbl(Y(1) – Y(2)))
d(2) = TriSides(CDbl(x(2) – x(3)), CDbl(Y(2) – Y(3)))
d(3) = TriSides(CDbl(x(1) – x(3)), CDbl(Y(1) – Y(3)))
Area = Heron(d(1), d(2), d(3))
O(1) = TriSides(CDbl(x(1)), CDbl(Y(1)))
O(2) = TriSides(CDbl(x(2)), CDbl(Y(2)))
O(3) = TriSides(CDbl(x(3)), CDbl(Y(3)))
a(1) = Heron(d(1), O(1), O(2))
a(2) = Heron(d(2), O(2), O(3))
a(3) = Heron(d(3), O(1), O(3))
If CSng(Area) = CSng(a(1) + a(2) + a(3)) Then
Answer = Answer + 1
End If
Next_i:
Next i
Debug.Print Answer; ” Time:”; Timer – T
End Sub
Those are the escape code for > and < . Can’t figure out how to get it right inside the vb blocks. We can escape it outside, but we truncate if we use the angle brackets, and don’t render the brackets right if we escape. It’s been a conundrum for a while.
I used two functions, TriSides(), which returns the square root of the sum or difference of two squares, so I could do Pythagorean math, and Heron(), which implements Heron’s law:
If IsMissing(Sum) Then Sum = True
If Sum Then
TriSides = Sqr(a ^ 2 + b ^ 2)
Else
TriSides = Sqr(Abs(a ^ 2 – b ^ 2)) ‘ order doesn’t matter
End If
End Function
Function Heron(a As Double, b As Double, C As Double) As Double
‘Use Heron ‘s formula for the area of a triangle
‘given the lengths a,b,c of the sides
‘A = sqrt s*(s-a)*(s-b)*(s-c)
‘where s is the semiperimeter 0.5*(a+b+c).
Dim s As Double
s = 0.5 * (a + b + C)
Heron = Sqr(s * (s – a) * (s – b) * (s – C))
End Function
I admire those who saw the application of Heron going in. One thing that I don’t understand is why I had to convert doubles to singles at the end for the areas to total per Heron. I presume is has to do with the square root routine, but I invite comment. Code ran in one-tenth the time of 102A, at about .01 seconds.
…mrt
This one for me was the easiest one yet. I already had a UDF that checks if a specified 2D point is inside any polygon specified by a series of points, so it was just a matter of:
Open the text file as a comma delimited file
Use the Index function to get the triange coordinates into a 3×2 array
Write a quick and dirty sub to step the index from 1 to 1000, and count the number of times the Inside() function returned 1.
The Inside() function can be downloaded from: http://newtonexcelbach.wordpress.com/2008/08/10/intersections-interpolations-and-rotations/
along with a load of other geometric and interpolation functions.
The sub was:
Sub Triangles()
Dim i As Long, Sum As Long
For i = 1 To 1000
[G3] = i
Sum = Sum + [L3]
[M3] = Sum
Next i
End Sub
Maybe not very pretty or very fast, but the overall solution time (including coding and thinking) was excellent, and that’s my personal goal with these things.
That said, I’m interested in the area approach and will be having a closer look at that.
Michael – if you calculate two “equal” values by two different methods the last figure in the results is likely to be different due to rounding errors. You can subtract the results and check that the absolute value of the difference is less than some very small number (like 1e-12), or you can round them to some lower precision and check that the rounded values are equal, which is effectively what you did by converting from double to single.
Hi Doug –
Thank you two times. I’ll go get your routine. I found out I didn’t know how to write one.
…mrt
“I found out I didn’t know how to write one”
It’s one of those things that is harder than it looks.
Probably the easiest way (which I used) is to project a horizontal line from the point you are checking to infinity and count how many times it crosses the boundary of the shape.
An odd number of crossings indicates inside, and zero or an even number indicates outside.
But you get situations like intersecting a corner, or a segment of the shape exactly overlying the projected line from the point, so you have to have a consistent way of dealing with those things.
Having said I was more interested in the total workin time than in shaving a few milliseconds of the computer time, I decided to see if I could shave a few milliseconds the computer time.
The function below uses the area method, but finds the areas without finding the side lengths, thus saving a lot of squaring and square rooting. On my machine it runs in 7.8 milliseconds, compared with about 15 milliseconds for the routine potsed by Michael.
Replace .LT. with the Less Than symbol, and enter as an array function in two adjacent cells to get the solution time as well as the result.
Function Triangles2(Coords As Variant) As Variant
Dim NumTri As Long, i As Long, j As Long, Area1 As Double, Area2 As Double
Dim XDiff(1 To 3) As Double, Ysum(1 To 3) As Double, Area2inc As Double
Dim count As Long, ResA(1 To 1, 1 To 2) As Double
Const MinDiff As Double = 0.0000000001
ResA(1, 2) = Timer
Coords = Coords.Value2
NumTri = UBound(Coords)
For i = 1 To NumTri
Area1 = 0
Area2 = 0
Area2inc = 0
XDiff(1) = Coords(i, 3) – Coords(i, 1)
XDiff(2) = Coords(i, 5) – Coords(i, 3)
XDiff(3) = Coords(i, 1) – Coords(i, 5)
Ysum(1) = (Coords(i, 2) + Coords(i, 4)) / 2
Ysum(2) = (Coords(i, 4) + Coords(i, 6)) / 2
Ysum(3) = (Coords(i, 6) + Coords(i, 2)) / 2
For j = 1 To 3
Area1 = Area1 + XDiff(j) * Ysum(j)
Next j
Area1 = Abs(Area1)
XDiff(2) = -Coords(i, 3)
XDiff(3) = Coords(i, 1)
Ysum(2) = (Coords(i, 4)) / 2
Ysum(3) = (Coords(i, 2)) / 2
For j = 1 To 3
Area2inc = Area2inc + XDiff(j) * Ysum(j)
Next j
Area2 = Abs(Area2inc)
XDiff(1) = Coords(i, 3)
XDiff(2) = Coords(i, 5) – Coords(i, 3)
XDiff(3) = -Coords(i, 5)
Ysum(1) = (Coords(i, 4)) / 2
Ysum(2) = (Coords(i, 4) + Coords(i, 6)) / 2
Ysum(3) = (Coords(i, 6)) / 2
Area2inc = 0
For j = 1 To 3
Area2inc = Area2inc + XDiff(j) * Ysum(j)
Next j
Area2 = Area2 + Abs(Area2inc)
XDiff(1) = -Coords(i, 1)
XDiff(2) = Coords(i, 5)
XDiff(3) = Coords(i, 1) – Coords(i, 5)
Ysum(1) = (Coords(i, 2)) / 2
Ysum(2) = (Coords(i, 6)) / 2
Ysum(3) = (Coords(i, 6) + Coords(i, 2)) / 2
Area2inc = 0
For j = 1 To 3
Area2inc = Area2inc + XDiff(j) * Ysum(j)
Next j
Area2 = Area2 + Abs(Area2inc)
If Abs(Area2 – Area1) .LT. MinDiff Then count = count + 1
Next i
ResA(1, 1) = count
ResA(1, 2) = Timer – ResA(1, 2)
Triangles2 = ResA
End Function
Good morning, Doug –
I like it when you’re workin’. I take it you used 1/2(base*height), and while I think I see 1/2 in there, I can’t vision the area computation loop. Could I impose upon you to talk thru one of the XDiff and YSum/2 loops? I do see how you employed the mindiff concept. Thanks for that.
…mrt
Before searching the ‘Net, I thought of two ways to check if the origin was in the triangle. One was to ensure that each axis was intercepted at least once. And, since we were dealing with straight lines, there would be only 1 intercept.
The other method was to ensure that the angles formed around the origin with each pair of vertices must equal 360 degrees.
But, I rejected both as being sufficiently compled to compute given just the vertex coordinates that there had to be a better approach.
A Google search pointed to an approach using vector products. This would be easy to compute.
In fact, the computations are sufficiently straightforward that the problem can be solved in Excel without any code support.
http://www.tushar-mehta.com/misc_tutorials/project_euler/euler102.html
Michael – This blog post:
http://newtonexcelbach.wordpress.com/2008/03/09/section-properties-from-coordinates/
shows how the area calculation works for an irregular quadrilateral, but the same principle applies to a triangle (or any other polygon). It’s actually based on the trapezoid rule (area = base x average height) rather than the 1/2(base * height).
I took a more visual (Venn diagram) approach. I created two Base Classes: Coordinates and Triangle. The Triangle base class was made up of three Coordinates properties called Vertices1 through Vertices3. Then, I created a method called InTheTriangle which received an arbitrary Coordinate as a parameter. In this case, it is always the Origin. Then, I checked to see if the Origin was either on one of the sides of the triangle or if it was on the side of the line connecting two of the Vertices as the third Vertex. Surprisingly (to me), this approach yielded the correct answer.
Here is the main function as well as the CalcDirection helper function:
Public Function InTheTriangle(Pt As Coordinates) As Boolean
Dim InOut As Directioning
InTheTriangle = False
m = CalcSlope(Vertices1, Vertices2): b = CalcIntercept(Vertices1, m)
InOut = CalcDirection(m, b, Pt)
If InOut = Equal Then
InTheTriangle = True
Else
d = CalcDirection(m, b, Vertices3)
If d = InOut Then
m = CalcSlope(Vertices2, Vertices3): b = CalcIntercept(Vertices2, m)
InOut = CalcDirection(m, b, Pt)
If InOut = Equal Then
InTheTriangle = True
Else
d = CalcDirection(m, b, Vertices1)
If d = InOut Then
m = CalcSlope(Vertices3, Vertices1): b = CalcIntercept(Vertices3, m)
InOut = CalcDirection(m, b, Pt)
If InOut = Equal Then
InTheTriangle = True
Else
d = CalcDirection(m, b, Vertices2)
InTheTriangle = (d = InOut)
End If
End If
End If
End If
End If
End Function
Private Function CalcDirection(Slope As Variant, Intercept As Double, Vertex3 As Coordinates) As Directioning
Dim Result As Double
If CStr(Slope) = Infinity Then
Result = Vertex3.X + Intercept
Else
Result = Vertex3.Y – (Slope * Vertex3.X) – Intercept
End If
If Abs(Result) 0 Then
CalcDirection = Greater
Else
CalcDirection = Less
End If
End Function
Me and a friend established that the origin is enclosed by the triangle if and only if one of the edges crossed the Y axis at a negative Y value and another edge crossed it at a positive Y value. Quite a simple rule that seems to make sense.
Now I need to figure out the supporting maths…
I’d figure the quickest approach would be identification first for quick results, then calculations when needed.
1. If any of the points is the origin (0,0), done — TRUE.
2. Check quadrants: if either all x or all y values are all either positive or negative, then done — FALSE.
3. At this point one point must have an x value with a different sign than the other 2 points, and it’s the lines between this point and each of the other 2 points that must have y intercepts on either side of the origin. Let the single point p0 (x0, y0) and the other points p1 (x1, y1) and p2 (x2, y2). The y intercepts between p0 and each of p1 and p2 are given by
(y0-x0*(y1-y0)/(x1-x0))
and
(y0-x0*(y2-y0)/(x2-x0))
and these must have different signs (one or the other could be exactly 0).