14
\$\begingroup\$

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.

triangle_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:

enter image description here

enter image description here

enter image description here

Current top score: Sp3000 - 7000 points - C
\$\endgroup\$
1
  • \$\begingroup\$ Comments are not for extended discussion; this conversation has been moved to chat. \$\endgroup\$
    – Doorknob
    Jul 15, 2015 at 20:22

9 Answers 9

7
\$\begingroup\$

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[1]);
    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;
}
\$\endgroup\$
2
  • \$\begingroup\$ 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. \$\endgroup\$ Dec 17, 2016 at 19:35
  • \$\begingroup\$ 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/… \$\endgroup\$ Dec 17, 2016 at 21:27
6
\$\begingroup\$

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[1]);
    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.)

\$\endgroup\$
2
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ Jul 15, 2015 at 13:30
  • \$\begingroup\$ @RetoKoradi Ah thanks - that was one surprisingly not-minor f, actually. \$\endgroup\$
    – Sp3000
    Jul 15, 2015 at 13:46
3
\$\begingroup\$

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"
FLX_DEF_THREAD_FRAME
//Thread Frame Constructor
thread_frame_t::thread_frame_t(
) :
  gcp(0),
  shape_list_head(&thread_frame_t_ptr_map)
{}

//-----------------------------------------
//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[1])).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
\$\endgroup\$
2
\$\begingroup\$

C# (about 11000?)

using System;
using System.Collections.Generic;

namespace PPCG
{
    class PPCG53100
    {
        static void Main(string[] args)
        {
            int n = int.Parse(args[0]);
            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.

\$\endgroup\$
1
  • \$\begingroup\$ 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. \$\endgroup\$ Mar 30, 2017 at 18:30
0
\$\begingroup\$

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[724]={0};
 if(barr[0]==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[1]==0){P("uso: questo_programma.exe  arg1\n ove arg1 e\' un numero positivo\n");R 0;}
 ni=atoi(a[1]);
 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?

\$\endgroup\$
1
0
\$\begingroup\$

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 
  
\$\endgroup\$
0
\$\begingroup\$

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)  [[14951]]
                                                  Type: List List Integer
        Time: 19.22 (IN) + 36.95 (EV) + 0.05 (OT) + 3.62 (GC) = 59.83 sec
\$\endgroup\$
0
\$\begingroup\$

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
\$\endgroup\$
1
  • 1
    \$\begingroup\$ Welcome to PPCG! \$\endgroup\$ Mar 18, 2018 at 13:37
0
\$\begingroup\$

Rust, brute force, 4500

On my computer, the program calculated that the answer at n=4500 is 7844386, within 60 seconds.

use std::sync::mpsc;
use std::thread;
use std::time::Duration;
use rayon::prelude::*;

#[inline]
fn is_square(v: i32) -> bool {
    let s = (f32::sqrt(v as f32) + 0.5).floor() as i32;
    s * s == v
}

#[inline]
fn is_median(v: i32) -> bool {
    v % 4 == 0 && is_square(v / 4)
}

#[inline]
fn calculate_triangles(n: i32) -> i32 {
    (1..=n)
        .into_par_iter()
        .map(|c| {
            let mut count = 0;
            for b in (c + 1) / 2..=c {
                for a in (c - b + 1)..=b {
                    if is_median(2 * (b * b + c * c) - a * a)
                        || is_median(2 * (a * a + c * c) - b * b)
                        || is_median(2 * (a * a + b * b) - c * c)
                    {
                        count += 1;
                    }
                }
            }
            count
        })
        .sum()
}

fn main() {
    let n: i32 = 4500;
    let (tx, rx) = mpsc::channel();

    let _computation_thread = thread::spawn(move || {
        let n_tri = calculate_triangles(n);
        tx.send(n_tri).unwrap();
    });

    match rx.recv_timeout(Duration::from_secs(60)) {
        Ok(n_tri) => println!("Result: {}", n_tri),
        Err(_) => println!("Timed out after 60 seconds."),
    }

    // The thread will be terminated when the main function ends, regardless of whether it has finished its computation.
}

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.