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\$ Commented 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\$ Commented 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\$ Commented 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\$ Commented Jul 15, 2015 at 13:30
  • \$\begingroup\$ @RetoKoradi Ah thanks - that was one surprisingly not-minor f, actually. \$\endgroup\$ Commented 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\$ Commented 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
  • \$\begingroup\$ (since C99) \$\endgroup\$ Commented Aug 16, 2017 at 20:34
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\$ Commented 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\$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.