# Integral triangles and integral medians

Consider a triangle ABC where each side has integer length (an integral triangle). Define a median of ABC to be a line segment from a vertex to the midpoint of the opposing side. In the figure below, the red line segments represent the medians. Note that any given triangle has three medians. Let n be some positive integer. How many non-degenerate integral triangles with each side length less than or equal to n have at least one integral median?

### Challenge

Write a program to compute the number of integral triangles with at least one integral median for a given maximum side length n. The order of the side lengths does not matter, i.e. <6,6,5> represents the same triangle as <5,6,6> and should be counted only once. Exclude degenerate triangles such as <1,2,3>.

### Scoring

The largest n for which your program can generate the number of triangles in 60 seconds on my machine is your score. The program with the highest score wins. My machine is a Sony Vaio SVF14A16CLB, Intel Core i5, 8GB RAM.

### Examples

Let T(N) be the program with input N.

T(1) = 0
T(6) = 1
T(20) = 27
T(22) = 34


Note that T(1) = T(2) = T(3) = T(4) = T(5) = 0 because no combination of integral sides will yield an integral median. However, once we get to 6, we can see that one of the medians of the triangle <5,5,6> is 4, so T(6) = 1.

Note also that T(22) is the first value at which double-counting becomes an issue: the triangle <16,18,22> has medians 13 and 17 (and 2sqrt(85)).

### Computing the medians

The medians of a triangle can be calculated by the following formulas:   Current top score: Sp3000 - 7000 points - C

• Comments are not for extended discussion; this conversation has been moved to chat. – Doorknob Jul 15 '15 at 20:22

# C, brute force - n=6080

This is more a baseline than a serious contender, but at least it should get things started.

n=6080 is as high as I got in a minute of runtime on my own machine, which is a MacBook Pro with an Intel Core i5. The result I got for this value is:

15041226

The code is purely brute force. It enumerates all the triangles within the size limit, and tests for the condition:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

static inline int isSquare(int v) {
int s = (int)(sqrtf((float)v) + 0.5f);
return s * s == v;
}

static inline int isMedian(int v) {
return v % 4 == 0 && isSquare(v / 4);
}

int main(int argc, char* argv[]) {
int n = atoi(argv);
int nTri = 0;
int a, b, c;

for (c = 1; c <= n; ++c) {
for (b = (c + 1) / 2; b <= c; ++b) {
for (a = c - b + 1; a <= b; ++a) {
if (isMedian(2 * (b * b + c * c) - a * a) ||
isMedian(2 * (a * a + c * c) - b * b) ||
isMedian(2 * (a * a + b * b) - c * c)) {
++nTri;
}
}
}
}

printf("%d\n", nTri);

return 0;
}

• Depending on the compiler, you can get faster + better round-to-nearest from using lrintf() or (int)roundf() instead of adding 0.5f and using the default truncation. Sometimes you need to use -ffast-math to get it to compile to a single cvtss2si instruction, though. gcc inlines lrintf() and sqrtf with only -fno-math-errno, so you get efficient asm: godbolt.org/g/E3hncQ. (I used -march=ivybridge because that's the OP's CPU). With -ffast-math, clang turns the sqrt into a rsqrt + Newton iteration; IDK if that's a win. – Peter Cordes Dec 17 '16 at 19:35
• Oops, usually not roundf. Use (int)nearbyintf() if lrintf() doesn't inline, because it uses the current rounding mode instead of a specific weird one. stackoverflow.com/questions/37620659/… – Peter Cordes Dec 17 '16 at 21:27

# C, approx 6650 6900

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

static inline int is_square(int n) {
if ((n&2) != 0 || (n&7) == 5 || (n&11) == 8) {
return 0;
}

int s = (int) (sqrtf((float) n) + 0.5f);
return (s*s == n);
}

int main(int argc, char **argv) {
int n = atoi(argv);
int count = 0;

for (int a = 1; a <= n; ++a) {
if (a&1) {
for (int b = (a+1)/2; b <= a; ++b){
if (b&1) {
for (int c = a-b+2; c <= b; c += 2) {
if (is_square((a*a + b*b)/2 - (c*c)/4)) {
++count;
}
}
} else {
for (int c = a-b+2; c <= b; c += 2) {
if (is_square((a*a + c*c)/2 - (b*b)/4)) {
++count;
}
}
}
}
} else {
for (int b = (a+1)/2; b <= a; ++b){
if (b&1) {
for (int c = a-b+2; c <= b; c += 2) {
if (is_square((b*b + c*c)/2 - (a*a)/4)) {
++count;
}
}
} else {
for (int c = a-b+2; c <= b; c += 2) {
if (is_square((b*b + c*c)/2 - (a*a)/4) ||
is_square((c*c + a*a)/2 - (b*b)/4) ||
is_square((a*a + b*b)/2 - (c*c)/4)) {
++count;
}
}
}
}
}
}

printf("%d\n", count);
return 0;
}


I don't really use C often, but with the amount of arithmetic going on it seemed like a good choice of language. The core algorithm is brute force like @RetoKoradi's answer, but with a few simple optimisations. I'm not sure our values are comparable though, because @RetoKoradi's computer seems to be faster than mine.

The major optimisation is bypassing the % 4 check completely. An integer square n*n is either 0 or 1 modulo 4, depending on whether n itself is 0 or 1 modulo 2. Thus, we can take a look at all possibilities for (x, y, z) % 2:

x%2  y%2  z%2    (2*(x*x+y*y) - z*z) % 4
----------------------------------------
0    0    0              0
0    0    1              3
0    1    0              2
0    1    1              1
1    0    0              2
1    0    1              1
1    1    0              0
1    1    1              3


Conveniently, there are only two cases to consider: (0, 0, 0) and (1, 1, 0), which, given the first two sides a, b, equates to the third side c having parity a^b:

 a%2   b%2         c%2 must be
-----------------------------
0     0               0
0     1               1
1     0               1
1     1               0


a^b is the same parity as a-b, so rather than searching from c = a-b+1 and going up by 1s, this lets us search from c = a-b+2 and go up by 2s.

Another optimisation comes from the fact that, for the (1, 1, 0) case, we only need to call is_square once since only one permutation works. This is special cased in the code by unrolling the search.

The other optimisation included is simply a quickfail in the is_square function.

Compilation was done with -std=c99 -O3.

(Thanks to @RetoKoradi for pointing out that the 0.5 in is_square needed to be 0.5f to avoid a double conversion taking place.)

• Very minor, but you may want to use 0.5f instead of 0.5 in is_square(). 0.5 is a constant of type double, so the expression will produce a double value when you add 0.5, including type conversion from float to double for the other term. – Reto Koradi Jul 15 '15 at 13:30
• @RetoKoradi Ah thanks - that was one surprisingly not-minor f, actually. – Sp3000 Jul 15 '15 at 13:46

# Felix, unknown

fun is_square(v: int) => let s = int$sqrt$ v.float + 0.5f in s*s == v;
fun is_median(v: int) => v % 4 == 0 and (v/4).is_square;

proc main() {
n := int\$ System::argv 1;
var ntri = 0;

for var c in 1 upto n do
for var b in (c+1)/2 upto c do
for var a in c - b + 1 upto b do
if is_median(2*(b*b+c*c)-a*a) or
is_median(2*(a*a+c*c)-b*b) or
is_median(2*(a*a+b*b)-c*c) do ++ntri; done
done
done
done

ntri.println;
}

main;


Basically a port of the C answer, but it's faster than it, tested with clang -O3 and icc -O3. Felix and Nim are literally the only two languages I know of that can beat C and C++ at benchmarks. I'm working on a parallel version, but it'll be a bit till it's finished, so I decided to post this ahead.

I also put "unknown" because my computer isn't necessarily the fastest on earth...

Command used to build:

flx --usage=hyperlight -c --static -o sl0 sl0.flx


The generated C++ is pretty interesting to look at:

//Input file: /home/ryan/golf/itri/sl0/sl0.flx
//Generated by Felix Version 15.04.03
//Timestamp: 2015/7/16 20:59:42 UTC
//Timestamp: 2015/7/16 15:59:42 (local)
#define FLX_EXTERN_sl0 FLX_EXPORT
#include "sl0.hpp"
#include <stdio.h>
#define comma ,

//-----------------------------------------
//EMIT USER BODY CODE
using namespace ::flxusr::sl0;

//-----------------------------------------
namespace flxusr { namespace sl0 {

//-----------------------------------------
//DEFINE OFFSET tables for GC
#include "sl0.rtti"
) :
gcp(0),
{}

//-----------------------------------------
//DEFINE FUNCTION CLASS METHODS
#include "sl0.ctors_cpp"
//------------------------------
//C PROC <61624>: _init_
void _init_(FLX_APAR_DECL_ONLY){
int _i63436_v63436_s;
int _i63435_v63435_s;
int s;
int a;
int b;
int c;
int ntri;
int n;
n = static_cast<int>(::std::atoi((::std::string(1<0||1>=PTF argc?"":PTF argv)).c_str())); //assign simple
ntri = 0; //assign simple
c = 1; //assign simple
_63421:;
if(FLX_UNLIKELY((n < c))) goto _63428;
b = (c + 1 ) / 2 ; //assign simple
_63422:;
if(FLX_UNLIKELY((c < b))) goto _63427;
a = (c - b ) + 1 ; //assign simple
_63423:;
if(FLX_UNLIKELY((b < a))) goto _63426;
/*begin match*/
/*match case 1:s*/
s  = static_cast<int>((::std::sqrt(((static_cast<float>(((2 * (b * b  + (c * c ) )  - (a * a ) ) / 4 ))) + 0.5f ))))/*int.flx: ctor*/; //init
/*begin match*/
/*match case 1:s*/
_i63435_v63435_s  = static_cast<int>((::std::sqrt(((static_cast<float>(((2 * (a * a  + (c * c ) )  - (b * b ) ) / 4 ))) + 0.5f ))))/*int.flx: ctor*/; //init
/*begin match*/
/*match case 1:s*/
_i63436_v63436_s  = static_cast<int>((::std::sqrt(((static_cast<float>(((2 * (a * a  + (b * b ) )  - (c * c ) ) / 4 ))) + 0.5f ))))/*int.flx: ctor*/; //init
if(!((((2 * (b * b  + (c * c ) )  - (a * a ) ) % 4  == 0) && (s * s  == (2 * (b * b  + (c * c ) )  - (a * a ) ) / 4 )  || (((2 * (a * a  + (c * c ) )  - (b * b ) ) % 4  == 0) && (_i63435_v63435_s * _i63435_v63435_s  == (2 * (a * a  + (c * c ) )  - (b * b ) ) / 4 ) ) ) || (((2 * (a * a  + (b * b ) )  - (c * c ) ) % 4  == 0) && (_i63436_v63436_s * _i63436_v63436_s  == (2 * (a * a  + (b * b ) )  - (c * c ) ) / 4 ) ) )) goto _63425;
{
int* _tmp63490 = (int*)&ntri;
++*_tmp63490;
}
_63425:;
if(FLX_UNLIKELY((a == b))) goto _63426;
{
int* _tmp63491 = (int*)&a;
++*_tmp63491;
}
goto _63423;
_63426:;
if(FLX_UNLIKELY((b == c))) goto _63427;
{
int* _tmp63492 = (int*)&b;
++*_tmp63492;
}
goto _63422;
_63427:;
if(FLX_UNLIKELY((c == n))) goto _63428;
{
int* _tmp63493 = (int*)&c;
++*_tmp63493;
}
goto _63421;
_63428:;
{
_a12344t_63448 _tmp63494 = ::flx::rtl::strutil::str<int>(ntri) + ::std::string("\n") ;
::flx::rtl::ioutil::write(stdout,_tmp63494);
}
}

//-----------------------------------------
}} // namespace flxusr::sl0
//CREATE STANDARD EXTERNAL INTERFACE
FLX_FRAME_WRAPPERS(::flxusr::sl0,sl0)
FLX_C_START_WRAPPER_PTF(::flxusr::sl0,sl0,_init_)

//-----------------------------------------
//body complete


using System;
using System.Collections.Generic;

namespace PPCG
{
class PPCG53100
{
static void Main(string[] args)
{
int n = int.Parse(args);
Console.WriteLine(CountOOE(n) + CountEEE(n));
}

static int CountOOE(int n)
{
// Maps from a^2 + b^2 to (b - a, a + b), which are the exclusive bounds on c.
IDictionary<int, List<Tuple<int, int>>> pairs = new Dictionary<int, List<Tuple<int, int>>>();

for (int a = 1; a <= n; a += 2)
{
int k = 2 * a * a;
for (int b = a; b <= n; b += 2, k += 4 * (b - 1))
{
List<Tuple<int, int>> prev;
if (!pairs.TryGetValue(k, out prev)) pairs[k] = prev = new List<Tuple<int, int>>();
prev.Add(Tuple.Create(b - a, a + b));
}
}

int max = 2 * n * n;
int count = 0;
for (int x = 1; x <= n >> 1; x++)
{
int k = 4 * x * x;
for (int y = x; y <= n; y++, k += 4 * y - 2)
{
if (k > max) break;
List<Tuple<int, int>> ab;
if (pairs.TryGetValue(k, out ab))
{
foreach (var pair in ab)
{
// Double-counting isn't possible if a, b are odd.
if (pair.Item1 < x << 1 && x << 1 < pair.Item2)
{
count++;
}
if (x != y && y << 1 <= n && pair.Item1 < y << 1 && y << 1 < pair.Item2)
{
count++;
}
}
}
}
}

return count;
}

static int CountEEE(int n)
{
// Maps from a^2 + b^2 to (b - a, a + b), which are the exclusive bounds on c.
IDictionary<int, List<Tuple<int, int>>> pairs = new Dictionary<int, List<Tuple<int, int>>>();

for (int a = 2; a <= n; a += 2)
{
int k = 2 * a * a;
for (int b = a; b <= n; b += 2, k += 4 * (b - 1))
{
List<Tuple<int, int>> prev;
if (!pairs.TryGetValue(k, out prev)) pairs[k] = prev = new List<Tuple<int, int>>();
prev.Add(Tuple.Create(b - a, a + b));
}
}

// We want to consider m in the range [1, n] and c/2 in the range [1, n/2]
// But to save dictionary lookups we can scan x in [1, n/2], y in [x, n] and consider both ways round.
int max = 2 * n * n;
int count = 0;
for (int x = 1; x <= n >> 1; x++)
{
int k = 4 * x * x;
for (int y = x; y <= n; y++, k += 4 * y - 2)
{
if (k > max) break;
List<Tuple<int, int>> ab;
if (pairs.TryGetValue(k, out ab))
{
foreach (var pair in ab)
{
// (c1, m1) = (2x, y)
// (c2, m2) = (2y, x)

int a = (pair.Item2 - pair.Item1) / 2, b = (pair.Item2 + pair.Item1) / 2;
int c1 = 2 * x;

if (pair.Item1 < c1 && c1 < pair.Item2)
{
// To deduplicate: the possible sets of integer medians are:
//     m_c
//     m_a, m_c
//     m_b, m_c
//     m_a, m_b, m_c
// We only want to add if c is (wlog) the shortest edge whose median is integral (or joint integral in case of isosceles triangles).

if (c1 <= a) count++;
else if (!IsIntegerMedian(b, c1, a))
{
if (c1 <= b || !IsIntegerMedian(a, c1, b)) count++;
}
}

int c2 = 2 * y;
if (c1 != c2 && c2 <= n && pair.Item1 < c2 && c2 < pair.Item2)
{
if (c2 <= a) count++;
else if (!IsIntegerMedian(b, c2, a))
{
if (c2 <= b || !IsIntegerMedian(a, c2, b)) count++;
}
}
}
}
}
}

return count;
}

private static bool IsIntegerMedian(int a, int b, int c)
{
int m2 = 2 * (a * a + b * b) - c * c;
int s = (int)(0.5f + Math.Sqrt(m2));
return ((s & 1) == 0) && (m2 == s * s);
}
}
}


n is taken as a command-line argument.

### Explanation

We can rewrite $m = \sqrt{(2a^2 + 2b^2 - c^2) / 4}$ as $2a^2 + 2b^2 = 4m^2 + c^2$, whence it's obvious that $c^2$ must be even, and so $c$ is even. Let $c = 2C$ and we rewrite again as $a^2 + b^2 = 2(m^2 + C^2)$. Therefore $a^2 + b^2$ must be even, so $a$ and $b$ must have the same parity.

The equation $a^2 + b^2 = 2(m^2 + C^2)$ is the basis for the meet-in-the-middle algorithm employed here.

If $a$ and $b$ are odd then we have no risk of double-counting, because only one of the three medians can possibly be integral. If all three are even then we need to beware double-counting. Therefore I handle the two cases separately so that the odd-odd-even case can be processed faster than the even-even-even case.

• I can't build Felix on my machine, but my times for n=5000 are 67 seconds for Reto Koradi's answer, 48 seconds for Sp3000's answer, and 13 seconds for my answer. – Peter Taylor Mar 30 '17 at 18:30

# C, n=3030 here

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define R     return
#define u32 unsigned
#define F        for
#define P     printf

int isq(u32 a)
{u32 y,x,t,i;
static u32  arr720[]={0,1,4,9,16,25,36,49,64,81,100,121,144,169,196,225,256,289,324,361,400,441,484,529,576,625,676,180,241,304,369,436,505,649,160,409,496,585,340,544,145,601,244,580,481,640,385,265};
static char barr={0};
if(barr==0)F(i=0;i<(sizeof arr720)/sizeof(unsigned);++i)
if(arr720[i]<720) barr[arr720[i]]=1;
if(barr[a%720]==0) R 0;
y=sqrt(a);
R y*y==a;
}

int f(u32 a, u32 b, u32 c)
{u32 t,x;
if(c&1)R 0;
t= a*a+b*b;
if(t&1)R 0;
R isq((2*t-c*c)/4);
}

int h(u32 n)
{u32 cnt,a,c,k,ke,kc,d,v,l,aa,bb,cc;

cnt=0;
F(a=1;a<=n;++a)
{ke=(n-a)/2;
F(k=0;k<=ke;++k)
{v=a+k;
d=v*v+k*k;
l=sqrt(d);
v=n/2;
if(l>v)l=v;
v=a+k-1;
if(l>v)l=v;
F(c=k+1;c<=l;++c)
{if(isq(d-c*c))
{bb=a+2*k;cc=2*c;
if(bb>cc && f(a, cc,bb)) continue;
if( a>cc && f(cc,bb, a)) continue;
++cnt;
//P("|a=%u b=%u c=%u", a, bb, cc);
}
}
}
}
R cnt;
}

int main(int c, char** a)
{time_t  ti, tf;
double   d;
int     ni;
u32    n,i;

if(c!=2||a==0){P("uso: questo_programma.exe  arg1\n ove arg1 e\' un numero positivo\n");R 0;}
ni=atoi(a);
if(ni<=0){P("Parametro negativo o zero non permesso\n");R 0;}
n=ni;
if(n>0xFFFFF){P("Parametro troppo grande non permesso\n"); R 0;}
F(i=3;i<33;++i)if(i<10||i>21)P("T(%u)=%u|",i, h(i));
ti=time(0);
P("\nT(%u)=%u\n", n, h(n));
tf=time(0);
d=difftime(tf,ti);
P("Tempo trascorso = %.2f sec\n", d);
R 1;
}


results:

C:\Users\a\b>prog 3030
T(3)=0|T(4)=0|T(5)=0|T(6)=1|T(7)=1|T(8)=2|T(9)=3|T(22)=34|T(23)=37|T(24)=42|T(25)=
45|T(26)=56|T(27)=59|T(28)=65|T(29)=67|T(30)=74|T(31)=79|T(32)=91|
T(3030)=3321226
Tempo trascorso = 60.00 sec


the above code would be the traslation in C of the Axiom answer (if we not count the isq() function).

My compiler not link a function others use sqrtf()... here there is no sqrt function for float... Are they sure that sqrtf it is a C standard function?

# APL NARS, n=239 282 in 59 seconds

f←{(a b c)←⍵⋄1=2∣c:0⋄t←+/a b*2⋄1=2∣t:0⋄0=1∣√4÷⍨(2×t)-c*2}

∇r←g n;cnt;c;a;k;kc;ke;d;l;bb;cc
r←⍬⋄cnt←0
:for a :in 1..n
ke←⌊(n-a)÷2
:for k :in 0..ke
d←((a+k)*2)+k*2
kc←⌊⌊/(n÷2),(a+k-1),√d
→B×⍳kc<k+1
:for c :in (k+1)..kc
→C×⍳∼1e¯9>1∣√d-c*2
bb←a+2×k⋄cc←2×c
→C×⍳(bb>cc)∧f a  cc bb
→C×⍳( a>cc)∧f cc bb  a
cnt+←1
⍝r←r,⊂a bb cc
C:     :endfor
B:  :endfor
:endfor
r←r,cnt
∇


(i traslate the Axiom answer one, in APL) test:

  g 282
16712
v←5 6 10 20 30 41
v,¨g¨v
5 0  6 1  10 4  20 27  30 74  41 166



# Axiom, n=269 in 59 sec

isq?(x:PI):Boolean==perfectSquare?(x)

f(a:PI,b:PI,c:PI):Boolean==
c rem 2=1=>false
t:=a^2+b^2
t rem 2=1=>false
x:=(2*t-c^2)quo 4
isq?(x)

h(n)==
cnt:=0  -- a:=a   b:=(a+2*k)  c:=
r:List List INT:=[]
for a in 1..n repeat
ke:=(n-a)quo 2
for k in 0..ke repeat
d:=(a+k)^2+k^2 -- (a^2+b^2)/2=(a+k)^2+k^2   m^2+c^2=d
l:=reduce(min,[sqrt(d*1.), n/2.,a+k-1])
kc:=floor(l)::INT
for c in k+1..kc repeat
if isq?(d-c^2) then
bb:=a+2*k; cc:=2*c
if bb>cc and f(a,cc,bb) then iterate   -- 2<->3
if  a>cc and f(cc,bb,a) then iterate   -- 1<->3
cnt:=cnt+1
--r:=cons([a,a+2*k,2*c],r)
r:=cons([cnt],r)
r


If a,b,cx are length of the sides of one triangle of max lenght side n...

We would know that m:=sqrt((2*(a^2+b^2)-cx^2)/4)

(1) m^2=(2*(a^2+b^2)-cx^2)/4


As Peter Taylor had said, 4|(2*(a^2+b^2)-cx^2) and because 2|2*(a^2+b^2) than 2|cx^2 => cx=2*c. So from 1 will be

(2) m^2=(a^2+b^2)/2-c^2


a, and b has to have the same parity, so we could write b in function of a

(3) a:=a   b:=(a+2*k)


than we have that

(4)(a^2+b^2)/2=(a^2+(a+2*k)^2)/2=(a+k)^2+k^2


so the (1) can be rewritten see (2)(3)(4) as:

m^2+c^2=(a+k)^2 + k^2=d         a:=a  b:=(a+2*k)  cx:=2*c


where

a in 1..n
k in 0..(n-a)/2
c in k+1..min([sqrt(d*1.), n/2.,a+k-1])


results

(16) -> h 269
(16)  []
Type: List List Integer
Time: 19.22 (IN) + 36.95 (EV) + 0.05 (OT) + 3.62 (GC) = 59.83 sec


VBA 15,000 in TEN seconds!

I expected much less after these other posts. On an Intel 7 with 16 GB RAM I get 13-15,000 in TEN seconds. On a Pentium with 4 GB RAM, I get 5-7,000 in TEN seconds. The code is below. Here is the latest result on the Pentium

abci= 240, 234, 114, 7367, 147
abci= 240, 235, 125, 7368, 145
abci= 240, 236, 164, 7369, 164
abci= 240, 238, 182, 7370, 221
abci= 240, 239, 31, 7371, 121


It got up to a triangle with sides 240, 239, 31 and a medium of 121. The count of mediums is 7,371.

Sub tria()
On Error Resume Next
Dim i As Long, a As Integer, b As Integer, c As Integer, ma As Double, mb As Double, mc As Double, ni As Long, mpr As Long
Dim dtime As Date
dtime = Now
Do While Now < DateAdd("s", 10, dtime)  '100 > DateDiff("ms", dtime, Now) '
a = a + 1
' Debug.Assert a < 23
b = 1: c = 1
Do
ma = 0
If a < b + c And b < a + c And c < a + b Then
ma = ((2 * b ^ 2 + 2 * c ^ 2 - a ^ 2) / 4) ^ 0.5
If ma <> 0 Then ni = i + 1 * -1 * (0 = ma - Fix(ma))
If ni > i Then
If ma <> mpr Then
i = ni
mpr = ma
Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & ma
GoTo NextTri  'TO AVOID DOUBLE COUNTING
End If
End If
'End If

mb = 0
'If b < a + c Then
mb = ((2 * a ^ 2 + 2 * c ^ 2 - b ^ 2) / 4) ^ 0.5
If mb <> 0 Then ni = i + 1 * -1 * (0 = mb - Fix(mb))
If ni > i Then
If mb <> mpr Then
i = ni
mpr = mb
Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & mb
GoTo NextTri  'TO AVOID DOUBLE COUNTING
End If
End If
'End If

mc = 0
'IfThen
mc = ((2 * b ^ 2 + 2 * a ^ 2 - c ^ 2) / 4) ^ 0.5
If mc <> 0 Then ni = i + 1 * -1 * (0 = mc - Fix(mc))
If ni > i Then
If mc <> mpr Then
i = ni
mpr = mc
Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & mc
End If
End If
End If
NextTri:
Do While c <= b
'c = c + 1
ma = 0
If a < b + c And b < a + c And c < a + b Then

ma = ((2 * b ^ 2 + 2 * c ^ 2 - a ^ 2) / 4) ^ 0.5
If ma <> 0 Then ni = i + 1 * -1 * (0 = ma - Fix(ma))
If ni > i Then
If ma <> mpr Then
mpr = ma
i = ni
End If
Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & ma
GoTo NextTri2  'TO AVOID DOUBLE COUNTING
End If
'End If

mb = 0
'If b < a + c Then
mb = ((2 * a ^ 2 + 2 * c ^ 2 - b ^ 2) / 4) ^ 0.5
If mb <> 0 Then ni = i + 1 * -1 * (0 = mb - Fix(mb))
If ni > i Then
If mb <> mpr Then
mpr = mb
i = ni
Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & mb
GoTo NextTri2  'TO AVOID DOUBLE COUNTING
End If
End If
'End If

mc = 0
'If c < b + a Then
mc = ((2 * b ^ 2 + 2 * a ^ 2 - c ^ 2) / 4) ^ 0.5
If mc <> 0 Then ni = i + 1 * -1 * (0 = mc - Fix(mc))
If ni > i Then
If mc <> mpr Then
mpr = mc
i = ni
Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & mc
End If
End If
End If
' Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i
c = c + 1
Loop 'While c <= a
NextTri2:
b = b + 1
c = 1
Loop While b <= a
Loop
Debug.Print i

End Sub

• Welcome to PPCG! – Martin Ender Mar 18 '18 at 13:37