# Happy New Year problem.....

1. Jan 1, 2016

### Cosmos

Let N be the no. of right angled triangles with integer sides anf inradius r=2016.Then N is divisible by which of the following:
2, 3 , 5 ,7

2. Jan 1, 2016

### collinsmark

Could you repeat that (perhaps with a touch more clarity)? I don't quite follow what the problem means.

Edit: I think I understand what you mean now. If I'm interpreting the problem correctly, you're asking,

Let N be the number of Pythagorean triples, corresponding to right triangles with integer length sides (all sides, including the hypotenuse, are integers), that have an inradius of 2016. This number N is evenly divisible by which of the following: 2, 3, 5, 7.​

Last edited: Jan 1, 2016
3. Jan 3, 2016

### collinsmark

I have a solution.

In this solution I am assuming that any congruent triangles are the same and only count once. Thus they can be represented in the form where $a \lt b \lt c$. In other words, a 3-4-5 triangle is the same as 4-3-5 is the the same as 5-4-3, etc, and they all can be called 3-4-5. Any set of such congruent triangles that have an inradius of 2016 only counts toward N once.

I wrote a computer program.

If you choose to solve this problem by writing a computer program, take special care that if you use the Pythagorean theorem, $c^2 = a^2 + b^2$, that you are using 64 bit integer types (e.g., type "long" in many programming languages) to represent the numbers. The squared terms, on occasion, will be too large to be represented by standard 32 bit integer types (e.g., type "int").

The least prime factor of N is 2. So "2" is the short answer.
(Assuming that my calculations are correct, of course.)

[Edit: although 3 and 5 are also factors of N. So you could also say that 2, 3 and 5 all evenly divide N. Not 7 though.]

The inradius, $r$, of a right triangle, where $c$ is the hypotenuse, can be found using $2r = a + b - c$.
The computer program uses this equation to "calculate" $a$.

The program then loops through all possible values of $c$ and $b$ (calculating $a$ each time) and checks whether the sides form a right triangle via the Pythagorean theorem. A couple of optimizations were also done:

Upper limit of the hypotenuse:
We can calculate the upper limit of the hypotenuse, $c$, by noting that as $c$ increases, $a$ approaches $2r$. If $a = 2r$, then $b$ and $c$ become infinite, parallel lines. So the smallest possible $a$ of interest in this problem is just a tad greater than that. Thus
$a_{min} = 2r + 1$.

An alternative approach is to note that as $c$ increases, $b$ approaches $c$. This leads us to say,
$b_{max} = c_{max} - 1$.

Using either of these approaches, and using a touch of algebra substitution together with $c^2 = a^2 + b^2$ and $2r = a + b - c$, we obtain,
$c_{max} = 2r^2 + 2r + 1$.

Important optimization:
For any given hypotenuse, $c$, we can reduce the number of $b$s searched through for Pythagorean testing. It might be tempting to just ignore all $b$s that are less than $c/2$. But using that optimization would still cause the program to run for many hours, even days perhaps, before completion, even on a fast computer (assuming single threaded operation). We can do better.

Note that as $c$ increases, $b$ also increases monotonically, but at a faster rate. From this we can conclude that the function "$c - b$" decreases monotonically with increasing $c$.

So the optimization that I used ignores all values of $b$ that are less than $c - (c_{previous} - b_{previous})$, where $c_{previous}$ and $b_{previous}$ are from the values of $c$ and $b$ of the last (i.e., previous) successfully found Pythagorean triple. This optimization reduces the running time of the program from hours or days to merely seconds.

So here is the program. It was written as a C# console application using Visual Studio 2015.
Code (Text):
using System;
using System.Collections.Generic;

namespace _2016_New_Years_Puzzle_2
{
class Program
{
static void Main(string[] args)
{
using (System.IO.StreamWriter oStream = new System.IO.StreamWriter("Output.txt"))
{
int r = 2016;
// The hypotensuse need not be any bigger than
// 2r^2 + 2r + 1.
int maxHypotenuse = (2 * r * r) + (2 * r) + 1;
List<ValidTriplet> TripletList = new List<ValidTriplet>();
long lastGoodCmB = (int)(1.4142 * r); // ~(r)sqrt(2) sounds like a good place to start to ensure a < b.
// c is surely greater than 2r, so let's start with c = 2r and go from there.
for (long c = 2 * r; c <= maxHypotenuse; c++)
{
// Our starting b for this loop need not be smaller than the previously successful c - b.
// This is because both c and b are monotonically increasing, but b increases faster.
for (long b = c - lastGoodCmB; b < c; b++)
{
// 2r = a + b - c
long a = (2 * r) + c - b;

if ((a * a) + (b * b) == (c * c))
{
// Got one! Add it to the list.
ValidTriplet triplet = new ValidTriplet(a, b, c);

// Make comment to screen and output file.
Console.WriteLine("a = " + a + ", b = " + b + ", c = " + c);
oStream.WriteLine("a = " + a + ", b = " + b + ", c = " + c);

// Update last good c - b.
lastGoodCmB = c - b;

// Since we've found a triple for this hypotenuse, let's
// move on to the next one.
break;
}
}
}

Console.WriteLine(Environment.NewLine + "Number of triples = " + TripletList.Count + Environment.NewLine + "Press a key to exit.");
oStream.WriteLine(Environment.NewLine + "Number of triples = " + TripletList.Count);
}
}
}

class ValidTriplet
{
// Assumed that a < b < c.
// Also assumed that a, b, and c are all positive.
public long a;
public long b;
public long c;

public ValidTriplet(long a, long b, long c)
{
this.a = a;
this.b = b;
this.c = c;
}
}
}
And the program's output:
a = 6720, b = 7056, c = 9744
a = 6678, b = 7104, c = 9750
a = 6624, b = 7168, c = 9760
a = 6384, b = 7488, c = 9840
a = 6336, b = 7560, c = 9864
a = 6300, b = 7616, c = 9884
a = 6080, b = 8001, c = 10049
a = 6048, b = 8064, c = 10080
a = 5824, b = 8568, c = 10360
a = 5796, b = 8640, c = 10404
a = 5760, b = 8736, c = 10464
a = 5600, b = 9216, c = 10784
a = 5568, b = 9324, c = 10860
a = 5544, b = 9408, c = 10920
a = 5376, b = 10080, c = 11424
a = 5355, b = 10176, c = 11499
a = 5328, b = 10304, c = 11600
a = 5208, b = 10944, c = 12120
a = 5184, b = 11088, c = 12240
a = 5166, b = 11200, c = 12334
a = 5056, b = 11970, c = 12994
a = 5040, b = 12096, c = 13104
a = 4928, b = 13104, c = 14000
a = 4914, b = 13248, c = 14130
a = 4896, b = 13440, c = 14304
a = 4816, b = 14400, c = 15184
a = 4800, b = 14616, c = 15384
a = 4788, b = 14784, c = 15540
a = 4704, b = 16128, c = 16800
a = 4680, b = 16576, c = 17224
a = 4620, b = 17856, c = 18444
a = 4608, b = 18144, c = 18720
a = 4599, b = 18368, c = 18935
a = 4544, b = 19908, c = 20420
a = 4536, b = 20160, c = 20664
a = 4480, b = 22176, c = 22624
a = 4473, b = 22464, c = 22905
a = 4464, b = 22848, c = 23280
a = 4424, b = 24768, c = 25160
a = 4416, b = 25200, c = 25584
a = 4410, b = 25536, c = 25914
a = 4368, b = 28224, c = 28560
a = 4356, b = 29120, c = 29444
a = 4326, b = 31680, c = 31974
a = 4320, b = 32256, c = 32544
a = 4288, b = 35784, c = 36040
a = 4284, b = 36288, c = 36540
a = 4256, b = 40320, c = 40544
a = 4248, b = 41664, c = 41880
a = 4228, b = 45504, c = 45700
a = 4224, b = 46368, c = 46560
a = 4221, b = 47040, c = 47229
a = 4200, b = 52416, c = 52584
a = 4194, b = 54208, c = 54370
a = 4179, b = 59328, c = 59475
a = 4176, b = 60480, c = 60624
a = 4160, b = 67536, c = 67664
a = 4158, b = 68544, c = 68670
a = 4144, b = 76608, c = 76720
a = 4140, b = 79296, c = 79404
a = 4130, b = 86976, c = 87074
a = 4128, b = 88704, c = 88800
a = 4116, b = 100800, c = 100884
a = 4113, b = 104384, c = 104465
a = 4104, b = 116928, c = 117000
a = 4096, b = 131040, c = 131104
a = 4095, b = 133056, c = 133119
a = 4088, b = 149184, c = 149240
a = 4086, b = 154560, c = 154614
a = 4081, b = 169920, c = 169969
a = 4080, b = 173376, c = 173424
a = 4074, b = 197568, c = 197610
a = 4068, b = 229824, c = 229860
a = 4064, b = 258048, c = 258080
a = 4060, b = 294336, c = 294364
a = 4059, b = 305088, c = 305115
a = 4056, b = 342720, c = 342744
a = 4053, b = 391104, c = 391125
a = 4050, b = 455616, c = 455634
a = 4048, b = 512064, c = 512080
a = 4046, b = 584640, c = 584654
a = 4044, b = 681408, c = 681420
a = 4041, b = 907200, c = 907209
a = 4040, b = 1020096, c = 1020104
a = 4039, b = 1165248, c = 1165255
a = 4038, b = 1358784, c = 1358790
a = 4036, b = 2036160, c = 2036164
a = 4035, b = 2713536, c = 2713539
a = 4034, b = 4068288, c = 4068290
a = 4033, b = 8132544, c = 8132545

Number of triples = 90

Since 90 is an even number, its least prime factor is "2". So the answer is "2".

[Edit: although 3 and 5 are also factors of 90. So you could also say that 2, 3 and 5 all evenly divide N. Not 7 though.]

Last edited: Jan 3, 2016
4. Jan 3, 2016

### collinsmark

I just came up with a completely different way to solve the puzzle.

This new way also involves writing a computer program, however.

But the algorithm is fundamentally different than the last one.

The new method generates a list of primitive Pythagorean triples using the B. Berggren tree method discussed here:
https://en.wikipedia.org/wiki/Tree_of_primitive_Pythagorean_triples
Primitive triples are generated until each branch in the tree corresponds to a primitive triple whose hypotenuse is greater than or equal to
$c_{max} = 2r^2 + 2r + 1$.

Then the inradius associated with each primitive triple is compared to the 2016 radius (technically I compare twice the primitive's inradius to 2 $\times$ 2016). if it divides evenly, the primitive Pythagorean triple is scaled by that ratio. If it doesn't divide evenly the primitive triple is thrown out.

So here's the program (again, written in C#)
Code (Text):
using System;
using System.Collections.Generic;

namespace _2016_12_03_thing
{
class Program
{
static void Main(string[] args)
{
int rMax = 2016;
int twoRmax = 2 * rMax;
int cMax = (2 * rMax * rMax) + (2 * rMax) + 1;
List<ValidTriple> primTriples = new List<ValidTriple>();

// Create matricies for triple generation
// See https://en.wikipedia.org/wiki/Tree_of_primitive_Pythagorean_triples
// for details.
Matrix A = new Matrix
(1, -2, 2,
2, -1, 2,
2, -2, 3);
Matrix B = new Matrix
(1, 2, 2,
2, 1, 2,
2, 2, 3);
Matrix C = new Matrix
(-1, 2, 2,
-2, 1, 2,
-2, 2, 3);

// Generate list of primitive Pythagorean triples.
List<ValidTriple> currentLevelTriples = new List<ValidTriple>();
currentLevelTriples.Add(new ValidTriple(3, 4, 5)); // seed triple.
bool finishedGenerating = false;
while (!finishedGenerating)
{
List<ValidTriple> nextLevelTriples = new List<ValidTriple>();
foreach (ValidTriple triple in currentLevelTriples)
{
// Add triple to total primitive triple list.

// Generate next level of triples based on this one.
}
// Clear out currentLevelTriples list.
currentLevelTriples.Clear();
// transfer nextLevelTriples to CurrentLevelTriples.
finishedGenerating = true; // Set this true for the moment.
foreach(ValidTriple triple in nextLevelTriples)
{
if (triple.c <= cMax)
{
finishedGenerating = false;
}
}
}

// if radius of a valid triple evenly divides rMax
// we scale it accordingly and put it in a new list.
List<ValidTriple> rMaxTriples = new List<ValidTriple>();
foreach(ValidTriple triple in primTriples)
{
if(twoRmax % triple.twoR == 0)
{
int scale = twoRmax / triple.twoR;
}
}

// Output list.
using (System.IO.StreamWriter oStream = new System.IO.StreamWriter("Output.txt"))
{
foreach(ValidTriple triple in rMaxTriples)
{
Console.WriteLine("a = " + triple.a + ", b = " + triple.b + ", c = " + triple.c);
oStream.WriteLine("a = " + triple.a + ", b = " + triple.b + ", c = " + triple.c);
}
Console.WriteLine(Environment.NewLine + "Triple count: " + rMaxTriples.Count + Environment.NewLine + "Press a key to exit.");
oStream.WriteLine(Environment.NewLine + "Triple count: " + rMaxTriples.Count);
}
}
}

class ValidTriple
{
public int a { get; private set; }
public int b { get; private set; }
public int c { get; private set; }
public int twoR { get; private set; }

public ValidTriple(int a, int b, int c)
{
this.a = a;
this.b = b;
this.c = c;
twoR = a + b - c;
}

public ValidTriple mult(int N)
{
return new ValidTriple(N * a, N * b, N * c);
}

public static ValidTriple operator *(int N, ValidTriple triple)
{
return triple.mult(N);
}
public static ValidTriple operator *(ValidTriple triple, int N)
{
return triple.mult(N);
}
}

class Matrix
{
int[,] elem;

public Matrix
(int a00, int a01, int a02,
int a10, int a11, int a12,
int a20, int a21, int a22)
{
elem = new int[3, 3];
elem[0, 0] = a00;
elem[0, 1] = a01;
elem[0, 2] = a02;
elem[1, 0] = a10;
elem[1, 1] = a11;
elem[1, 2] = a12;
elem[2, 0] = a20;
elem[2, 1] = a21;
elem[2, 2] = a22;
}

public ValidTriple mult(ValidTriple trip)
{
int a; int b; int c;

a = (elem[0, 0] * trip.a) + (elem[0, 1] * trip.b) + (elem[0, 2] * trip.c);
b = (elem[1, 0] * trip.a) + (elem[1, 1] * trip.b) + (elem[1, 2] * trip.c);
c = (elem[2, 0] * trip.a) + (elem[2, 1] * trip.b) + (elem[2, 2] * trip.c);

return new ValidTriple(a, b, c);
}

public static ValidTriple operator *(Matrix A, ValidTriple B)
{
return A.mult(B);
}
}
}

And here's the output:

a = 6048, b = 8064, c = 10080
a = 5040, b = 12096, c = 13104
a = 7056, b = 6720, c = 9744
a = 10080, b = 5376, c = 11424
a = 4704, b = 16128, c = 16800
a = 11088, b = 5184, c = 12240
a = 5544, b = 9408, c = 10920
a = 4536, b = 20160, c = 20664
a = 7560, b = 6336, c = 9864
a = 8736, b = 5760, c = 10464
a = 6624, b = 7168, c = 9760
a = 13104, b = 4928, c = 14000
a = 4788, b = 14784, c = 15540
a = 5208, b = 10944, c = 12120
a = 5600, b = 9216, c = 10784
a = 4896, b = 13440, c = 14304
a = 6384, b = 7488, c = 9840
a = 18144, b = 4608, c = 18720
a = 8568, b = 5824, c = 10360
a = 9324, b = 5568, c = 10860
a = 5166, b = 11200, c = 12334
a = 4914, b = 13248, c = 14130
a = 14616, b = 4800, c = 15384
a = 5796, b = 8640, c = 10404
a = 6300, b = 7616, c = 9884
a = 22176, b = 4480, c = 22624
a = 4368, b = 28224, c = 28560
a = 5328, b = 10304, c = 11600
a = 4816, b = 14400, c = 15184
a = 4464, b = 22848, c = 23280
a = 4320, b = 32256, c = 32544
a = 6678, b = 7104, c = 9750
a = 4620, b = 17856, c = 18444
a = 25200, b = 4416, c = 25584
a = 4410, b = 25536, c = 25914
a = 4284, b = 36288, c = 36540
a = 4599, b = 18368, c = 18935
a = 4473, b = 22464, c = 22905
a = 4680, b = 16576, c = 17224
a = 4256, b = 40320, c = 40544
a = 4424, b = 24768, c = 25160
a = 4356, b = 29120, c = 29444
a = 4326, b = 31680, c = 31974
a = 19908, b = 4544, c = 20420
a = 35784, b = 4288, c = 36040
a = 4248, b = 41664, c = 41880
a = 46368, b = 4224, c = 46560
a = 4200, b = 52416, c = 52584
a = 5355, b = 10176, c = 11499
a = 4221, b = 47040, c = 47229
a = 4228, b = 45504, c = 45700
a = 4176, b = 60480, c = 60624
a = 4158, b = 68544, c = 68670
a = 4179, b = 59328, c = 59475
a = 11970, b = 5056, c = 12994
a = 4194, b = 54208, c = 54370
a = 67536, b = 4160, c = 67664
a = 4144, b = 76608, c = 76720
a = 4140, b = 79296, c = 79404
a = 4128, b = 88704, c = 88800
a = 4116, b = 100800, c = 100884
a = 4130, b = 86976, c = 87074
a = 4104, b = 116928, c = 117000
a = 4113, b = 104384, c = 104465
a = 4095, b = 133056, c = 133119
a = 8001, b = 6080, c = 10049
a = 131040, b = 4096, c = 131104
a = 4088, b = 149184, c = 149240
a = 4086, b = 154560, c = 154614
a = 4080, b = 173376, c = 173424
a = 4081, b = 169920, c = 169969
a = 4074, b = 197568, c = 197610
a = 4068, b = 229824, c = 229860
a = 4064, b = 258048, c = 258080
a = 4060, b = 294336, c = 294364
a = 4059, b = 305088, c = 305115
a = 4056, b = 342720, c = 342744
a = 4053, b = 391104, c = 391125
a = 4050, b = 455616, c = 455634
a = 4048, b = 512064, c = 512080
a = 4046, b = 584640, c = 584654
a = 4044, b = 681408, c = 681420
a = 4041, b = 907200, c = 907209
a = 4040, b = 1020096, c = 1020104
a = 4039, b = 1165248, c = 1165255
a = 4038, b = 1358784, c = 1358790
a = 4036, b = 2036160, c = 2036164
a = 4035, b = 2713536, c = 2713539
a = 4034, b = 4068288, c = 4068290
a = 4033, b = 8132544, c = 8132545

Triple count: 90

...which is divisible by 2, 3, and 5; but not 7

Last edited: Jan 4, 2016
5. Jan 4, 2016

### PeroK

Solution without the aid of a computer programme!

Every pythagorean triple can be written (uniquely) as $2kst, k(t^2-s^2), k (t^2 + s^2)$ where $s, t$ are coprime and $t-s$ is odd.

$r = \frac{1}{2}(a + b - c) = ks(t-s)$

$ks(t-s) = 2016 = 2^5 3^2 7$

$k, s, t-s$ can be any factors of $2016$ as long as $s, t$ are coprime. Equivalent to $s, (t-s)$ are coprime. And $t-s$ is odd, which means $t-s$ cannot have any powers of $2$ as a factor.

Each choice of $k$ leaves $1, 2$ or $4$ choices for $t-s$, depending on whether $k$ takes out the factors of $7$ or $9$ or both I.e.

$k = 1, 2, 4, 8, 16, 32$ leaves $4$ choices for $t-s$. ($= 1, 7, 9$ or $63$). Hence $24$ triples.

$k = 3, 6, 12, 24, 48, 96$ likewise. Hence another $24$ triples.

$k = 9, 18, 36, 72, 144, 192$ gives $2$ choices for $t-s$. Another $12$ triples.

$k = 7, 14, 28, 56, 112, 224$ likewise. Another $12$ triples.

$k = 21, \dots$ likewise. Another $12$ triples.

$k = 63, \dots$ leaves $1$ choice for $t-s$. A final $6$ triples.

Grand total of $90$ triples.

6. Jan 4, 2016

### collinsmark

Nice! Thanks for that, @PeroK .

I had considered playing around with Euclid's formula, but at the time I figured it might have been over my head. A few days ago I wasn't even familar with the inradius of a triangle (or if I was at one time, I must have forgotten all about it). I only stumbled upon Euclid's formula in the process of researching this puzzle.

Your approach is very clear and concise to me, if I follow through it step by step. Thanks again!

I've learned a lot about Pythagorean triples in the last few days.

[Edit: I'm half-tempted to write another computer program using your method. Of course it would have to also work with other years, not just 2016. So I'd have to code up a prime factor generator. Even if I was sloppy with that, it would still be more computationally efficient than either of my two previous programs. We can assume that I don't have access to a quantum computer. (Prime factor generation is one of the few known strengths for possible applications of quantum computers [Shor's[/PLAIN] [Broken] algorithm] --but that's overkill for this puzzle). A brute force method should suit my needs adequately for this.]

Last edited by a moderator: May 7, 2017
7. Jan 5, 2016

### PeroK

You could write a program to calculate the number of triples for any given input. First, you'd have to find the prime factorisation of the input, and then codify the factor counting. You only need the number of prime factors and their powers. Everything follows from that.

8. Jan 6, 2016

### collinsmark

At long last, I wrote a computer program based on the advice from @PeroK . It's generalized and works for any "year" (where "year" is any positive integer greater than or equal to 2, up to some outlandish limit).

As you can tell from the code, this program is much more involved than my previous ones. It might be more efficient (it probably is), but I haven't checked.

I'm pretty proud of the prime number generator I wrote for it; it's fairly fast. It's based on the "Sieve of Sundaram."

The prime factoring algorithm is pretty straightforward. I didn't do much optimization there.

For the main algorithm of moving prime factors between $k$, $s$ and $(t - s)$, the only "real" optimization I made was keeping the even prime factors (2s) out of the $(t - s)$ bucket from the get-go. It turns out that other "true" optimizations such as keeping $s$ and $(t - s)$ coprime at all times was more trouble than it was worth for me. I haven't been exposed to enough combinatorics tools to tackle that. (Recall this needs to work for any "year" -- any combination of arbitrary prime factors.)

Instead, the program loops through all possible combinations of prime factors and variables (except the 2s are kept out of $(t - s)$ bucket, as mentioned earlier), and then tests whether any rules were violated, such as $s$ and $(t - s)$ not being coprime. If all the rules are satisfied, and it isn't a repeated combination, it makes a Pythagorean triple. Otherwise, it ignores the combination and moves onto the next.

Code (Text):
using System;
using System.Collections.Generic;

{
class Program
{
static void Main(string[] args)
{
ulong Number = 2016; // Radius of Pythagorean triples

using (System.IO.StreamWriter oStream
= new System.IO.StreamWriter("Output.txt"))
{
Console.WriteLine("Year: " + Number + Environment.NewLine);
oStream.WriteLine("Year: " + Number + Environment.NewLine);
}

// list of triples (ultimate goal of this program)
List<ValidTriple> tripleList = new List<ValidTriple>();

// Generate primes.
PrimeGenerator pGen = new PrimeGenerator(Number);

// Find prime factors of Number.
PrimeFactorGenerator facGen = new PrimeFactorGenerator(Number, pGen);

// Write prime factors to output file and screen.
using (System.IO.StreamWriter oStream
= new System.IO.StreamWriter("Output.txt", true))
{
Console.Write("Prime factors: ");
oStream.Write("Prime factors: ");
foreach (ulong prime in facGen.factors.Keys)
{
Console.Write(prime + "^" + facGen.factors[prime] + " ");
oStream.Write(prime + "^" + facGen.factors[prime] + " ");
}
Console.Write(Environment.NewLine + Environment.NewLine);
oStream.Write(Environment.NewLine + Environment.NewLine);
}

// Create an initialize the Combinator and LegalComboChecker
Combinator combinator = new Combinator(3, true);
LegalComboChecker legalChecker = new LegalComboChecker();
// Add prime factors into the combinator
foreach (ulong factor in facGen.factors.Keys)
for (int i = 0; i < facGen.factors[factor]; i++)

// The first combination is aldready present in the
// combinator, and is guaranteed to be valid.
// (We can check validity anyway though.) Create the
// first triple here.
if (legalChecker.IsLegal(combinator))
{
}

// Loop through all possible combinations. Check for
// valid combinations and create Phythorean triples
// of the valid ones.
while(!combinator.FinishedLooping)
{
// Increment state of combinator

// Check if combinator contains a valid triple.
// If so, generate triple and add it to the list.
if (legalChecker.IsLegal(combinator))
{
}
}

// Write the Pythagorean triples to the screen and file
using (System.IO.StreamWriter oStream
= new System.IO.StreamWriter("Output.txt", true))
{
foreach (ValidTriple triple in tripleList)
{
Console.WriteLine("a = " + triple.a
+ ", b = " + triple.b
+ ", c = " + triple.c);
oStream.WriteLine("a = " + triple.a
+ ", b = " + triple.b
+ ", c = " + triple.c);
}
Console.WriteLine(Environment.NewLine
+ "Total number of Pythagorean triples: "
+ tripleList.Count);
oStream.WriteLine(Environment.NewLine
+ "Total number of Pythagorean triples: "
+ tripleList.Count);
}
}
}

// Just a helper class with a static function
// that generates a ValidTriple from
// the current state of a Combinator. Assumes,
// Combinator.Buckets[0]: k
// Combinator.Buckets[1]: s
// Combinator.Buckets[2]: (t - s)
// Then generates (a, b, c) triple,
// a = k(t^2 - s^2)
// b = 2kst
// c = k(t^2 + s^2)
class GenerateTriple
{
public static ValidTriple GetTriple(Combinator combinator)
{
ulong k = combinator.Buckets[0].BallProduct;
ulong s = combinator.Buckets[1].BallProduct;
ulong tms = combinator.Buckets[2].BallProduct;
ulong t = tms + s;

ulong a = k * ((t * t) - (s * s));
ulong b = 2 * k * s * t;
ulong c = k * ((t * t) + (s * s));

return new ValidTriple(a, b, c);
}
}

// Checks for combinations in the associated
// Combinator instance whether it satisfies
// the requirements to form a Pythagorean
// triplet. It assumes that the associated
// Combinator.Buckets[0]: factors in k
// Combinator.Buckets[1]: factors in s
// Combinator.Buckets[2]: factors in tms
// The method IsLegal returns true if
// conditions are satisfied: s and tms
// must be coprime, and no factors of
// 2 can be in tms.
class LegalComboChecker
{
List<Tuple> previousGoodTuples;

// Constructor
public LegalComboChecker()
{
previousGoodTuples = new List<Tuple>();
}

public bool IsLegal(Combinator combinator)
{
// t must be greater than s for the combo
// to be valid. So make sure tms > 0
ulong tms = combinator.Buckets[2].BallProduct;
if (tms <= 0) return false;

// Check for '2' factors in the tms bucket.
if (combinator.Buckets[2].ContainsFactor(2))
return false;
List<ulong> sFactors = combinator.Buckets[1].GetFactors();
List<ulong> tmsFactors = combinator.Buckets[2].GetFactors();

// check if any sFactors are in tmsFactors
foreach (ulong factor in sFactors)
if (tmsFactors.Contains(factor)) return false;
// check if any tmsFactors are in sFactors
foreach (ulong factor in tmsFactors)
if (sFactors.Contains(factor)) return false;

// If we've gotten this far, it's a valid combination
// for a triple, unless it has been used before. So now
// we check for repeats
ulong k = combinator.Buckets[0].BallProduct;
ulong s = combinator.Buckets[1].BallProduct;
foreach (Tuple tuple in previousGoodTuples)
if ((tuple.k == k) && (tuple.s == s) && (tuple.tms == tms))
return false; // it's a repeat.

// And now we have a unique, valid combination.
return true;
}
}

// Just a set of three ulong numbers.
class Tuple
{
public ulong k;
public ulong s;
public ulong tms;

// Constructor
public Tuple(ulong k, ulong s, ulong tms)
{
this.k = k;
this.s = s;
this.tms = tms;
}
}

// Combinator creates a list of buckets. Each bucket
// can contain balls. Balls are added to the Combinator
// using the AddBall method, where a "factor" of a ball
// is specified; the ball's ID is assigned automatically.
// When a ball is added, Combinator
// automatically places it in a bucket of its own
// choosing. After all the balls
// have been added, Combinator can then, setp by step,
// loop through every possible combination of buckets
// and balls.
class Combinator
{
// List of buckets
public List<Bucket> Buckets { get; }
// Total number of balls in the combinator
public int BallCount { get { return _ballCount; } }
// Number of buckets in the combinator
public int BucketCount { get { return _bucketCount; } }
public bool FinishedLooping { get { return _loopingFinished; } }

private int _ballCount;
private int _bucketCount;
private int _nextID;
private bool _loopingStarted;
private bool _loopingFinished;
private List<Ball> _balls;
private bool _keep2factorFromLastBucket;

// Constructor
public Combinator(
int numberOfBuckets,
bool keep2factorFromLastBucket= false)
{
_bucketCount = numberOfBuckets;
Buckets = new List<Bucket>();
_balls = new List<Ball>();
for (int i = 0; i < _bucketCount; i++)
{
Bucket bucket = new Bucket();
bucket.ID = Buckets.IndexOf(bucket);
}
_ballCount = 0;
_nextID = 0;
_loopingStarted = false;
_loopingFinished = false;
_keep2factorFromLastBucket = keep2factorFromLastBucket;
}

// Takes the specified factor, creates a ball
// with that factor. Then assigns that ball
// a uniquie ID and puts it in the first bucket.
{
if (_loopingStarted)
{
throw new Exception("Too late to get new balls.");
}
else
{
// Make the ball and put it in the bucket
Ball ball = new Ball(factor, _nextID);
_nextID++;
_ballCount++;

// Also add it to a private list of balls.
}
}

// This is the main combinitorics (sort of) thing.
// With each call it advances the balls to a new
// unique combination of balls. All balls and
// buckets are considered distinguishable.
// Consideration of non-distigushable balls (i.e.,
// repeated prime factors) and legality of
// combinations (i.e., coprime or even/odd
// parameters) must be checked externally.
// Method returns 'true' when all combinations
// have been exhausted.
// A special optimization, specific for finding
// Pythagorean triples, was created to keep
// balls with a factor of '2' out of the last bin
// if _keep2factorFromLastBucket is set to 'true'.
{
// Set flag that looping has started.
_loopingStarted = true;
for(int i = 0; i < _ballCount; i++)
{
int modifiedBucketCount = _bucketCount;
if(_keep2factorFromLastBucket && (_balls[i].factor == 2))
{
// Special case where _balls[i] contains factor
// 2, and we want to keep it out of the tms
// bucket.
modifiedBucketCount = _bucketCount - 1;
}
if(_balls[i].BucketID < modifiedBucketCount - 1)
{
// Ball is not yet in last bucket.
// Place it in a higher bucket.
Buckets[_balls[i].BucketID].Remove(_balls[i].ID);
break; // We have a new combination. Quit loop.
}
else
{
// Ball is already in the last bucket. Move
// it to the first bucket, and let the loop
// handle the "carry" with the next ball.
Buckets[_balls[i].BucketID].Remove(_balls[i].ID);
}
}
// Check to see if the looping is finished. Return
foreach (Ball ball in _balls)
{
if (_keep2factorFromLastBucket && (ball.factor == 2))
{
if ((ball.BucketID < _bucketCount - 2))
return false;
}
else if (ball.BucketID < _bucketCount - 1)
return false;
}
// If we've got this far, the looping is finished.
_loopingFinished = true;
return true;
}
}

// Buckets hold balls.
class Bucket
{
// List of balls.
public List<Ball> Balls { get; set; }
// Returns the number of balls in the bucket.
public int BallCount { get { return Balls.Count; } }
// Returns the product of all Ball.factor of all
// balls in the bucket, multiplied together.
public ulong BallProduct
{
get
{
ulong tempProd = 1;
foreach (Ball ball in Balls)
tempProd *= ball.factor;
return tempProd;
}
}
public int ID { get; set; } // set externally

// Constructor
public Bucket()
{
Balls = new List<Ball>();
}

// Adds a ball to the bucket.
{
ball.BucketID = ID;
}

// Removes the first ball from the bucket with
// specified ID.
public void Remove (int ballID)
{
foreach (Ball ball in Balls)
if (ball.ID == ballID)
{
Balls.Remove(ball);
break;
}
}

// Returns true if a ball in the bucket
// contains a ball with the specified
// factor.
public bool ContainsFactor(uint factor)
{
foreach (Ball ball in Balls)
if (ball.factor == factor) return true;
return false;
}

// Returns a collection of unique factors in the
// bucket (factors are not repeated, even if
// multiple balls in the bucket contain that
// factor)
public List<ulong> GetFactors()
{
List<ulong> factors = new List<ulong>();
foreach (Ball ball in Balls)
if (!factors.Contains(ball.factor))
return factors;
}
}

// Simply stores a single prime factor (no repeats).
// Each ball can have its own ID. IDs are assigned
// externally, at the time of instantiation.
// BucketID can be changed externally at any time.
class Ball
{
public int ID { get; } // each ball can have it's own ID.
public ulong factor { get;} // prime factor associated with ball.
public int BucketID { get; set; } // changed externally.

public Ball (ulong factor, int ID)
{
this.ID = ID;
this.factor = factor;
}
}

// Generates a dictionary of prime factors.
//
// The TKey of the returned dicationary is the collection of prime numbers
// up to NumberToBeFactored.
// Each element of TValue corresponds to the
// number of particular primes that make up the factoring. For example,
// NumberToBeFactored = 20 will produce this result:
//   Key      Value
//    2        2
//    5        1
// The above meaning the prime factors contain two 2s and one 5.
// These get stored in the factors property.
class PrimeFactorGenerator
{
// Public properties
public Dictionary<ulong, int> factors { get; private set; }

// Private properties
private PrimeGenerator _primeGenerator;
private ulong _maxNumber;

// constructor
public PrimeFactorGenerator(ulong maxNumber, PrimeGenerator primeGenerator,
bool generatePrimeList = false)
{
_maxNumber = maxNumber;
_primeGenerator = primeGenerator;
if (generatePrimeList) primeGenerator.generatePrimes(maxNumber);

// Instantiate the dictionary of factors.
factors = new Dictionary<ulong, int>();

// Generate the factors.
generateFactors();
}

// generates the factors.
public void generateFactors()
{
ulong workingMaxNumber  = _maxNumber;
while (workingMaxNumber > 1)
{
foreach (uint prime in _primeGenerator.primes)
{
int repetition = 0;
while ((workingMaxNumber % prime) == 0)
{
repetition++;
workingMaxNumber = workingMaxNumber / prime;
}
if (repetition > 0) factors.Add(prime, repetition);
}
}
}
}

// Prime number generator that implements the Sieve of Sundaram.
// For that it generates a list of prime numbers up to the
// specified max integer of interest.
class PrimeGenerator
{
// Prime number list.
public List<ulong> primes { get; private set; }

// Constructor
public PrimeGenerator(ulong maxNumber)
{
generatePrimes(maxNumber);
}

// Main algorithm to generate primes
public void generatePrimes(ulong maxNumber)
{
primes = new List<ulong>();

// eventually tell us which numbers are prime,
// based on the individual bits in the array.
// The array comes in chucks of 32 bits (type int),
// which is why we divide by 32 and such. There's
// not much reason to alter the hardoded 32 to
// softcoded, since the size is fixed in the
// hardware.
uint[] boolArray = new uint[(maxNumber / 32) + 1];
ulong m = (maxNumber / 2) + 1;

// Fill boolArray with true.
for (ulong i = 0; i < (ulong)boolArray.Length; i++)
boolArray[i] = 0xblackff;

// Here's the first part of the Sieve of Sundaram.
for (ulong i = 1; i <= m; i++)
{
for (ulong j = i; j <= (m - i) / ((2 * i) + 1); j++)
{
ulong firstCalc = i + j + (2 * i * j);
ulong ByteIdx = fixedQuotient(firstCalc, 32);
uint BitIdx = (uint)fixedRemainder(firstCalc, 32);

// Here we clear that particular bit.
boolArray[ByteIdx]
= ~((~boolArray[ByteIdx]) | (uint)(1 << (int)BitIdx));
}
}

// Now let's get rid of that '1' in the prime list by zeroing
// out bit 0. (The Sieve of Sundaram seems to think 1 is prime for
// some reason.)
boolArray[0] = ~(~boolArray[0] | 1);

// And before we go on, let's add '2' to the actual prime list to
// start, since the Sieve of Sundaram only generates odd primes.

// The final part of the Sieve of Sundaram comes from the known prime
// at 2k + 1, where k represents each the remaining '1's.
// These are completely different 'i' and 'j' than before, by the way.
// Here, i and j are just looping over boolean array elements.
for (ulong i = 0; i < (ulong)boolArray.Length; i++)
{
for (int j = 0; j < 32; j++)
{
if (((boolArray[i] >> j) & (1)) == 1)
{
ulong k = (32 * i) + (ulong)j;
ulong foundPrime = (2 * k) + 1;
if (foundPrime > maxNumber) break;
}
}
}
}
// These methods calculate the quotient and remainder when doing
// fixed point division. They're important to figure out what
// uint block we need from the boolean array, and which bit from
// that block.
private ulong fixedQuotient(ulong numerator, ulong denominator)
{
return numerator / denominator;
}
private ulong fixedRemainder(ulong numerator, ulong denominator)
{
return numerator % denominator;
}
}

// Represents a valid Pythagorean triple. It's assumed to be valid.
class ValidTriple
{
public ulong a { get; private set; }
public ulong b { get; private set; }
public ulong c { get; private set; }

public ValidTriple(ulong a, ulong b, ulong c)
{
this.a = a;
this.b = b;
this.c = c;
}
}
}
And here's the program's output:

Year: 2016

Prime factors: 2^5 3^2 7^1

a = 6048, b = 8064, c = 10080
a = 5040, b = 12096, c = 13104
a = 4536, b = 20160, c = 20664
a = 4284, b = 36288, c = 36540
a = 4158, b = 68544, c = 68670
a = 4095, b = 133056, c = 133119
a = 4704, b = 16128, c = 16800
a = 4368, b = 28224, c = 28560
a = 4200, b = 52416, c = 52584
a = 4116, b = 100800, c = 100884
a = 4074, b = 197568, c = 197610
a = 4053, b = 391104, c = 391125
a = 10080, b = 5376, c = 11424
a = 7056, b = 6720, c = 9744
a = 5544, b = 9408, c = 10920
a = 4788, b = 14784, c = 15540
a = 4410, b = 25536, c = 25914
a = 4221, b = 47040, c = 47229
a = 4256, b = 40320, c = 40544
a = 4144, b = 76608, c = 76720
a = 4088, b = 149184, c = 149240
a = 4060, b = 294336, c = 294364
a = 4046, b = 584640, c = 584654
a = 4039, b = 1165248, c = 1165255
a = 22176, b = 4480, c = 22624
a = 13104, b = 4928, c = 14000
a = 8568, b = 5824, c = 10360
a = 6300, b = 7616, c = 9884
a = 5166, b = 11200, c = 12334
a = 4599, b = 18368, c = 18935
a = 4320, b = 32256, c = 32544
a = 4176, b = 60480, c = 60624
a = 4104, b = 116928, c = 117000
a = 4068, b = 229824, c = 229860
a = 4050, b = 455616, c = 455634
a = 4041, b = 907200, c = 907209
a = 4128, b = 88704, c = 88800
a = 4080, b = 173376, c = 173424
a = 4056, b = 342720, c = 342744
a = 4044, b = 681408, c = 681420
a = 4038, b = 1358784, c = 1358790
a = 4035, b = 2713536, c = 2713539
a = 4896, b = 13440, c = 14304
a = 4464, b = 22848, c = 23280
a = 4248, b = 41664, c = 41880
a = 4140, b = 79296, c = 79404
a = 4086, b = 154560, c = 154614
a = 4059, b = 305088, c = 305115
a = 4064, b = 258048, c = 258080
a = 4048, b = 512064, c = 512080
a = 4040, b = 1020096, c = 1020104
a = 4036, b = 2036160, c = 2036164
a = 4034, b = 4068288, c = 4068290
a = 4033, b = 8132544, c = 8132545
a = 6624, b = 7168, c = 9760
a = 5328, b = 10304, c = 11600
a = 4680, b = 16576, c = 17224
a = 4356, b = 29120, c = 29444
a = 4194, b = 54208, c = 54370
a = 4113, b = 104384, c = 104465
a = 18144, b = 4608, c = 18720
a = 11088, b = 5184, c = 12240
a = 7560, b = 6336, c = 9864
a = 5796, b = 8640, c = 10404
a = 4914, b = 13248, c = 14130
a = 4473, b = 22464, c = 22905
a = 8736, b = 5760, c = 10464
a = 6384, b = 7488, c = 9840
a = 5208, b = 10944, c = 12120
a = 4620, b = 17856, c = 18444
a = 4326, b = 31680, c = 31974
a = 4179, b = 59328, c = 59475
a = 46368, b = 4224, c = 46560
a = 25200, b = 4416, c = 25584
a = 14616, b = 4800, c = 15384
a = 9324, b = 5568, c = 10860
a = 6678, b = 7104, c = 9750
a = 5355, b = 10176, c = 11499
a = 5600, b = 9216, c = 10784
a = 4816, b = 14400, c = 15184
a = 4424, b = 24768, c = 25160
a = 4228, b = 45504, c = 45700
a = 4130, b = 86976, c = 87074
a = 4081, b = 169920, c = 169969
a = 131040, b = 4096, c = 131104
a = 67536, b = 4160, c = 67664
a = 35784, b = 4288, c = 36040
a = 19908, b = 4544, c = 20420
a = 11970, b = 5056, c = 12994
a = 8001, b = 6080, c = 10049

Total number of Pythagorean triples: 90

Last edited: Jan 7, 2016