Integralne trójkąty i integralne mediany

15

Rozważ trójkąt ABC, w którym każdy bok ma długość całkowitą ( integralny trójkąt ). Zdefiniuj medianę z ABC być odcinek od wierzchołka do punktu środkowego przeciwnej stronie. Na poniższym rysunku segmenty czerwonej linii przedstawiają mediany. Zauważ, że każdy trójkąt ma trzy mediany.

triangle_medians

Niech n będzie dodatnią liczbą całkowitą. Ile nieodegenerowanych trójkątów całkowitych o każdej długości boku mniejszej lub równej n ma co najmniej jedną integralną medianę?

Wyzwanie

Napisz program do obliczania liczby integralnych trójkątów z co najmniej jedną integralną medianą dla danej maksymalnej długości boku n . Kolejność długości boków nie ma znaczenia, tzn. <6,6,5> reprezentuje ten sam trójkąt co <5,6,6> i należy go policzyć tylko raz. Wyklucz zdegenerowane trójkąty, takie jak <1,2,3>.

Punktacja

Największy n, dla którego Twój program może wygenerować liczbę trójkątów w ciągu 60 sekund na mojej maszynie, to Twój wynik. Program z najwyższym wynikiem wygrywa. Moje urządzenie to Sony Vaio SVF14A16CLB, Intel Core i5, 8 GB pamięci RAM.

Przykłady

Niech T ( N ) będzie program z wejściem N .

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

Zauważ, że T (1) = T (2) = T (3) = T (4) = T (5) = 0, ponieważ żadna kombinacja całek boków nie da integralnej mediany. Gdy jednak dojdziemy do 6, widzimy, że jedna z median trójkąta <5,5,6> wynosi 4, więc T (6) = 1.

Zauważ również, że T (22) jest pierwszą wartością, przy której podwójne liczenie staje się problemem: trójkąt <16,18,22> ma mediany 13 i 17 (i 2sqrt (85)).

Obliczanie median

Mediany trójkąta można obliczyć za pomocą następujących wzorów:

wprowadź opis zdjęcia tutaj

wprowadź opis zdjęcia tutaj

wprowadź opis zdjęcia tutaj

Current top score: Sp3000 - 7000 points - C
Nicolás Siplis
źródło
Komentarze nie są przeznaczone do rozszerzonej dyskusji; ta rozmowa została przeniesiona do czatu .
Klamka

Odpowiedzi:

7

C, brutalna siła - n = 6080

Jest to bardziej punkt odniesienia niż poważny pretendent, ale przynajmniej powinien zacząć wszystko.

n = 6080 jest tak wysoki, jak w minutę w czasie wykonywania na własnym komputerze, jakim jest MacBook Pro z procesorem Intel Core i5. Wynik dla tej wartości to:

15041226

Kod jest czysto brutalną siłą. Wymienia wszystkie trójkąty w ramach limitu wielkości i testuje warunek:

#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;
}
Reto Koradi
źródło
W zależności od kompilatora można uzyskać szybsze i lepsze zaokrąglanie do najbliższego przy użyciu lrintf()lub (int)roundf()zamiast dodawania 0,5f i przy użyciu domyślnego obcięcia. Czasami jednak trzeba go użyć -ffast-math, aby skompilować go do pojedynczej cvtss2siinstrukcji. gcc inlines lrintf()i sqrtfonly with -fno-math-errno, dzięki czemu zyskujesz wydajność asm: godbolt.org/g/E3hncQ . (Użyłem, -march=ivybridgeponieważ to jest procesor OP). Z-ffast-math , clang zamienia sqrt w rsqrt + iteracja Newtona; IDK, jeśli to wygrana.
Peter Cordes,
Ups, zwykle nie roundf. Użyj (int)nearbyintf()if lrintf()nie inline, ponieważ używa bieżącego trybu zaokrąglania zamiast określonego dziwnego. stackoverflow.com/questions/37620659/…
Peter Cordes
6

C, około 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;
}

Tak naprawdę nie używam często C, ale przy dużej ilości arytmetyki wydawało się, że to dobry wybór języka. Podstawowym algorytmem jest brutalna siła jak odpowiedź @ RetoKoradi , ale z kilkoma prostymi optymalizacjami. Nie jestem jednak pewien, czy nasze wartości są porównywalne, ponieważ komputer @ RetoKoradi wydaje się być szybszy niż mój.

Główną optymalizacją jest % 4całkowite ominięcie kontroli. Kwadratem całkowitym n*njest albo 0, albo 1 modulo 4, w zależności od tego n, czy sam jest 0, czy 1 modulo 2. Zatem możemy przyjrzeć się wszystkim możliwościom (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

Dogodnie są tylko dwa przypadki do rozważenia: (0, 0, 0)i (1, 1, 0)które, biorąc pod uwagę dwie pierwsze strony a, b, równa się cparzystości trzeciej strony a^b:

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

a^bma taką samą parzystość jak a-b, więc zamiast przeszukiwać c = a-b+1i zwiększać się o 1s, pozwala nam to przeszukiwać c = a-b+2i zwiększać się o 1s .

Kolejna optymalizacja wynika z faktu, że w tym (1, 1, 0)przypadku wystarczy wywołać is_square tylko raz, ponieważ działa tylko jedna permutacja. Jest to szczególnie uwypuklone w kodzie przez rozwinięcie wyszukiwania.

Inna zawarta w tym optymalizacja to po prostu szybka is_squarefunkcja.

Kompilacja została zakończona -std=c99 -O3.

(Podziękowania dla @RetoKoradi za wskazanie, że 0.5in is_square musiało być, 0.5faby uniknąć podwójnej konwersji).

Sp3000
źródło
1
Bardzo niewielkie, ale możesz chcieć użyć 0.5fzamiast 0.5w is_square(). 0.5jest stałą typu double, więc wyrażenie doda podwójną wartość po dodaniu 0.5, w tym konwersję typu od floatdo doubledla drugiego terminu.
Reto Koradi,
@RetoKoradi Ach, dzięki - to był fnaprawdę zaskakująco mało znaczący .
Sp3000,
2

Felix, nieznany

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;

Zasadniczo port odpowiedzi C, ale jest szybszy od niego, przetestowany z clang -O3 i icc -O3. Felix i Nim są dosłownie jedynymi dwoma znanymi mi językami, które mogą pokonać C i C ++ w testach porównawczych. Pracuję nad wersją równoległą, ale trochę potrwa, dopóki się nie skończy, więc postanowiłem opublikować to wcześniej.

Umieściłem również „nieznane”, ponieważ mój komputer niekoniecznie jest najszybszy na świecie ...

Polecenie użyte do zbudowania:

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

Wygenerowane C ++ jest dość interesujące:

//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
kirbyfan64sos
źródło
2

C # (około 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 jest traktowany jako argument wiersza poleceń.

Wyjaśnienie

m=(2)za2)+2)b2)-do2))/42)za2)+2)b2)=4m2)+do2)do2)dodo=2)doza2)+b2)=2)(m2)+do2))za2)+b2)zab muszą mieć taką samą parzystości.

za2)+b2)=2)(m2)+do2)) stanowi podstawę zastosowanego tutaj algorytmu spełnienia w środku.

zab

Peter Taylor
źródło
Nie mogę zbudować Felixa na moim komputerze, ale moje czasy n=5000to 67 sekund na odpowiedź Reto Koradi, 48 sekund na odpowiedź Sp3000 i 13 sekund na moją odpowiedź.
Peter Taylor
0

C, n = 3030 tutaj

#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;
}

wyniki:

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

powyższy kod byłby tłumaczeniem w C odpowiedzi Axiom (jeśli nie policzymy funkcji isq ()).

Mój kompilator nie łączy funkcji, której inni używają sqrtf () ... tutaj nie ma funkcji sqrt dla float ... Czy są pewni, że sqrtf jest funkcją standardową w C?

RosLuP
źródło
(od C99)
Peter Taylor
0

APL NARS, n = 239 282 w 59 sekund

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
∇

(tłumaczę odpowiedź Axiom 1, w 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 
RosLuP
źródło
0

Aksjomat, n = 269 w 59 sek

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

Jeśli a, b, cx mają długość boków jednego trójkąta o maksymalnej długości boku n ...

Wiedzielibyśmy, że m: = sqrt ((2 * (a ^ 2 + b ^ 2) -cx ^ 2) / 4)

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

Jak powiedział Peter Taylor, 4 | (2 * (a ^ 2 + b ^ 2) -cx ^ 2) i ponieważ 2 | 2 * (a ^ 2 + b ^ 2) niż 2 | cx ^ 2 => cx = 2 * c. Więc od 1 będzie

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

a ib muszą mieć tę samą parzystość, abyśmy mogli napisać b w funkcji a

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

niż my to mamy

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

więc (1) można przepisać patrz (2) (3) (4) jako:

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

gdzie

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

wyniki

(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
RosLuP
źródło
0

15 000 VBA w ciągu dziesięciu sekund!

Spodziewałem się znacznie mniej po tych innych postach. Na procesorze Intel 7 z 16 GB RAM dostaję 13–15 000 w ciągu dziesięciu sekund. Na Pentium z 4 GB pamięci RAM dostaję 5-7 000 w ciągu dziesięciu sekund. Kod znajduje się poniżej. Oto najnowszy wynik 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

Dostał się do trójkąta o bokach 240, 239, 31 i medium 121. Liczba mediów to 7371.

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
George Let
źródło
1
Witamy w PPCG!
Martin Ender