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).