Come ottimizzare il mio codice per trovare il conteggio di tutte le mediane integrali per tutti i possibili triangoli integrali con un <= b <= c <= 100000?

Question

sto facendo la soluzione per questo problema da Euler Project Problem 513, Integral median:

ABC è un triangolo con i lati lati integrale a≤b≤c. mc è la mediana che collega C e il punto medio di AB. F (n) è il numero di tali triangoli con c≤n per cui mc ha anche lunghezza integrale. F (10) = 3 e F (50) = 165.

Trova F (100000).

Analyse:

1. a <= b <= c <= n == 100000
2. ABC è un triangolo così dovrebbe: abs(a-b) < c < a+b
3. Mc = sqrt(2 * a^2+ 2 * b^2 - c^2)/2wikipedia
4. Mc è intero così 2 * a^2+ 2 * b^2 - c^2 dovrebbe essere un quadrato perfetto e divisibile per 4.

Codice:

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

#define N 100000 
#define MAX(a,b) (((a)>(b))?(a):(b)) 

void main(){ 
    unsigned long int count = 0; 
    unsigned long int a,b,c; 
    double mc; 

    for (a = 1; a <= N; a++) { 
     printf("%lu\n", a); 
     for (b = a; b <= N; b++) 
      for (c = MAX(b, abs(b-a)); c <=N && c < a+b; c++){ 
       mc = sqrt(2 *a *a + 2 * b * b - c * c)/2.0; 
       if (mc-(unsigned long)mc == 0) 
        count++; 
      } 
    } 
    printf("\ncpt == %lu\n", count); 

} 

Issues:

Funziona bene per le piccole n ma la complessità della soluzione è troppo alto, suppongo che sia O(n^3) (mi sbaglio?) che impiegherà giorni per n = 100000.

Come potrei migliorare questo con un metodo matematico o algoritmico?

Aggiornamenti

ho avuto questi suggerimenti:

• Calcolo potere dei a al di fuori dei b/c loop e la potenza di b fuori c ciclo. Questo ha migliorato leggermente le prestazioni.
• c non può essere dispari. quindi a e b devono avere la stessa parità. Questo ha migliorato le prestazioni 4 volte.
• Utilizzo di thread per dividere il lavoro su molti core. Può migliorare di un fattore vicino al numero di core.
• Una soluzione matematica pubblicata in math.stackexchange. Si richiede O(N^5/2) per una soluzione di base e può raggiungere O(N^2) utilizzando O(N^2) di memoria. Non l'ho ancora testato.

Answer 1

Poiché si tratta di un problema di Project Euler, si suppone di poterlo fare in circa un minuto del tempo di elaborazione su un computer moderno.Non si attengono sempre a questo, ma indica che un tempo di esecuzione di k*n^2 o k*n^2*log(n) probabilmente va bene se la costante non è troppo male, ma probabilmente non è k*n^2.5 o k*n^3.

Come ha commentato SleuthEye, il lato c deve essere pari, altrimenti l'interno della radice quadrata dovrebbe essere dispari, quindi prendere la radice quadrata e dividere per 2 non potrebbe fare un numero intero.

È possibile semplificare l'equazione su 4(mc^2+(c/2)^2) = 2(a^2+b^2).

Ecco un approccio: creare due dizionari, sinistra e destra. Per ciascuno, lasciare che le chiavi siano possibili valori di quel lato dell'equazione e lasciare che i valori siano un elenco delle coppie (mc,c/2) o (a,b) che producono la chiave. Per il dizionario giusto, dobbiamo solo considerare dove a e b hanno la stessa parità e dove 1<=a<=b<=n. Per la sinistra, abbiamo solo bisogno di prendere in considerazione 1<=c/2<=n/2 e 1<=mc<=sqrt(3)/2 n dal 4mc^2 = 2a^2+2b^2-c^2 <= 3b^2 <=3n^2.

Quindi passare attraverso le possibili chiavi e confrontare gli elementi dei valori di ciascun dizionario, trovando il numero di coppie compatibili ((mc,c/2),(a,b)) dove b <= c < a+b. Questo passaggio interno non è un tempo costante, ma le lunghezze massime e medie degli elenchi non sono troppo lunghe. I modi di scrivere n come somma di due quadrati corrispondono approssimativamente (fino a unità) ai modi per calcolare n negli interi gaussiani, e proprio come il più grande number of factors of an integer non cresce troppo rapidamente, lo stesso è vero negli interi gaussiani. Questo passaggio richiede il tempo per qualsiasi epsilon>0. Quindi, il tempo di esecuzione totale è O(n^(2+epsilon)) per qualsiasi epsilon>0.

In pratica, se si esaurisce la memoria, è possibile creare dizionari parziali in cui le chiavi sono limitate per essere in intervalli particolari. Anche questo è in parallelo.