Emanuele Bonin - Ricerca Mattoni diofantei

/*
    linea di comando per la compilazione con cc65
    cl65 -t c64 -O --codesize 800 -Cl diofanto.c diofantoExtern.s -o DioFanto.prg

    il programma per c64 è scritto in maniera da sfruttare la pagina 0 tramite l'inclusione
    del file diofantoExtern.s che riporto:
    Diofanto.s
    ## Inizio File ##
.export _m1;
_m1 = $40

.export _m2;
_m2 = $42

.export _m3;
_m3 = $44

.export _T1;
_T1 = $46

.export _T2;
_T2 = $48

.export _T3;
_T3 = $50

.export _SPIGA;
_SPIGA = $52

.export _SPIGB;
_SPIGB = $56

.export _SPIGC;
_SPIGC = $60

.export _MoltiplicatoreMax;
_MoltiplicatoreMax = $64
    ## Fine File ##


*/

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

#define MMAX        15
#define MAXLEVEL    5

#define POKE(addr,val)     (*(unsigned char*) (addr) = (val))


unsigned int TPP[1100][3];
unsigned int TPP_Index = 0;
unsigned int HH;
unsigned int MM;
unsigned int Sec;
unsigned int Milli;
unsigned int totale = 0;
unsigned long int start_t, end_t, total_t;

#ifdef __C64__
extern unsigned int m1;
#pragma zpsym("m1")
extern unsigned int m2;
#pragma zpsym("m2")
extern unsigned int m3;
#pragma zpsym("m3")
extern unsigned int T1;
#pragma zpsym("T1")
extern unsigned int T2;
#pragma zpsym("T2")
extern unsigned int T3;
#pragma zpsym("T3")
extern unsigned int MoltiplicatoreMax;
#pragma zpsym("MoltiplicatoreMax")
extern unsigned long int SPIGA;
#pragma zpsym("SPIGA")
extern unsigned long int SPIGB;
#pragma zpsym("SPIGB")
extern unsigned long int SPIGC;
#pragma zpsym("SPIGC")
#else

unsigned int m3;
unsigned int m2;
unsigned int m1;

unsigned long int SPIGA, SPIGB, SPIGC;
unsigned int T1, T2, T3;

unsigned int MoltiplicatoreMax;

#endif // C64


// Matrice generatrice delle triplette
typedef struct M {
    unsigned int Mx[3][3]; // = {{1,2,2},{2,1,2},{2,2,3}};
} tM;
tM MT;

void AssignNext(int TPP_Index_Padre, int nMode, int Livello);
void AddNextLevel(int Padre, int Livello) {

    AssignNext(Padre, 1, Livello);
    AssignNext(Padre, 2, Livello);
    AssignNext(Padre, 3, Livello);

}

// Funzione che calcola i tre figli generati dal nodo From
// e assegna i valori al nodo To
// Viene utilizzato il teorema di Barnings
// node indica quali modifiche vanno apportate alla matrice
// iniziale per generare uno dei tre "figli" della tripletta "genitore"
void AssignNext(int TPP_Index_Padre, int nMode, int Livello) {
    unsigned int a, b, c, sw;
    unsigned int M1=1, M2=1;
    if (nMode != 2) {
        if (nMode == 1)
            M2 = -1;
        else
            M1 = -1;
    }

    a = MT.Mx[0][0] * TPP[TPP_Index_Padre][0] * M1 +
        MT.Mx[0][1] * TPP[TPP_Index_Padre][1] * M2 +
        MT.Mx[0][2] * TPP[TPP_Index_Padre][2];

    b = MT.Mx[1][0] * TPP[TPP_Index_Padre][0] * M1 +
        MT.Mx[1][1] * TPP[TPP_Index_Padre][1] * M2 +
        MT.Mx[1][2] * TPP[TPP_Index_Padre][2];

    c = MT.Mx[2][0] * TPP[TPP_Index_Padre][0] * M1 +
        MT.Mx[2][1] * TPP[TPP_Index_Padre][1] * M2 +
        MT.Mx[2][2] * TPP[TPP_Index_Padre][2];

    // La generazione mantiene sempre la stessa parità
    // a; /* dispari */
    // b; /* pari */
    // c; /* dispari*/

    /* Ordino la terna in modo da avere sempre a < b
       a discapito dell'ordinamento per parità (a dispari e b pari)
       in modo da poter eseguire le ricerche successive
    */
    if (b < a) {
        sw = b;
        b=a;
        a=sw;
    }

    ++TPP_Index;

    TPP[TPP_Index][0] = a;
    TPP[TPP_Index][1] = b;
    TPP[TPP_Index][2] = c;

    if (Livello < MAXLEVEL) {
        AddNextLevel(TPP_Index, Livello+1);
    }
}



// Inizializzazione delle variabili comuni
void Initialize() {
    TPP_Index = 0;
    TPP[TPP_Index][0] = 3;
    TPP[TPP_Index][1] = 4;
    TPP[TPP_Index][2] = 5;

    AddNextLevel(TPP_Index, 1);
}

int MCD(unsigned int n1, unsigned int n2) {
    /* algoritmo di Euclide */
    unsigned int a, b, mcd = 1;
    if (n1 > n2) {
        a = n1;
        b = n2;
    } else {
        a = n2;
        b = n1;
    }
    if ((a % b) == 0)
        mcd = b;
    else
        mcd = MCD(b, a % b);
    return mcd;
}




void WalkOnList() {



    unsigned short ContaTerne= 0;
    MoltiplicatoreMax=MMAX;



    // MoltiplicatoreMax contiene il numero massimo che potrà assumere il valore di
    // "scalatura" della terna pitagorica B C D2 (quella più piccola)
    // che sarà la prima ad essere scelta.
    // Questo comporterà un ciclo ulteriore fino a tale valore

    // +T3---+---------------------------+
    // |   / |              F          T1|
    // | D/  A                           |
    // | /   |                           |
    // +--B--+----C----------------------+
    //       |                         T2|
    //       B             E             |
    //       +---------------------------+
    // Ip:F>E>D
    // F>E
    // A^2+C^2 > B^2+C^2 => A^2>B^2 => A>B
    // E>D
    // B^2+C^2 > B^2+A^2 => C^2>A^2 => C>A
    // F>D
    // A^2+C^2 > B^2+A^2 => C^2>B^2 => C>B
    //
    // Quindi: C>A>B
    //
    // Per ogni faccia, con le lettere minuscole indico i due cateti (a,b) e l'ipotenusa c
    // tenendo presente che a<b<c. l'indice dopo ma lettera minuscola
    // indica la tera di appartenenza.
    // Inoltre identifico ogni faccia con una terna e generalizzo pensando alla terna come
    // una terna che possa non essere primitiva, quindi ogniuno degli elementi della terna
    // è moltiplicato per una stessa costante.
    // chiamiamo m1, m2, m3 le tre costanti che moltiplicano le terne primitive di
    // ogn'una delle tre facce, in questo modo possiamo scrivere:

    // T1 = A,C,F = a1*m1,b1*m1,c1*m1
    // T2 = B,C,E = a2*m2,b2*m2,c2*m2
    // T3 = B,A,D = a3*m3,b3*m3,c3*m3
    //
    // Gli spigoli sono ordinati secondo la loro grandezza
    // (1) A = a1*m1 = b3*m3
    // (2) B = a2*m2 = a3*m3
    // (3) C = b1*m1 = b2*m2
    //
    // lavorando sulle prime due delle ultime tre equivalenze sopra esposte:
    // dalla (1) m1 = ((b3*m3)/a1) = A/a1 => che a1|A
    // dalla (2) m2 = ((a3*m3)/a2) = B/a2 => che a2|B
    // dalla (3) m2 = ((b1*m1)/b2) = C/b2 => che b2|C
    //

    // girando nella maniera opportuna le equivalenze (1) (2) (3)otteniamo:
    // m3 = a1*m1/b3, m2 = a3*m3/a2, m1=b2*m2/b1
    // m3 = a1*(b2*m2/b1)/b3 => m3 = a1*(b2*(a3*m3/a2)/b1)/b3
    // => m3 = (a1*b2*a3*m3)/(a2*b1*b3) => m3/m3 = (a1*b2*a3)/(b1*a2*b3)
    // quindi abbiamo che un mattone diofanteo rispetta l'equivalenza:
    // 1 = (a1*b2*a3)/(b1*a2*b3)
    // il che ci permette di escludere le triplette primitive che non rispettino tale regola
    // prima ancora di andare a cercare i moltiplicatori opportuni.

    // Se analizziamo la parità dei tre spigoli, partendo da:
    // A^2+B^2=D^2
    // A^2+C^2=F^2
    // B^2+C^2=E^2
    // Troviamo che:
    // D^2+E^2+F^2 = 2(A^2+B^2+C^2)
    // il che indica che la somma dei quadrati delle diagonali deve essere pari
    // il che vuol dire che o sono tutte e tre diagonali pari oppure
    // ve ne sono due dispari e una sola pari.
    // d'altronde non può essere nemmeno che siano tutte e tre pari altrimenti
    // vorrebbe dire che gli elementi di tutte le terne pitagoriche dovrebbero
    // essere pari, ma questo non è possibile.
    // quindi due diagonali devono essere pari e una dispari.
    // questo ci evita di cercare il moltiplicatore m3 pari che renderebbe due diagonali pari

#ifndef __C64__
#define FOUT 1
#endif // __C64__
#ifdef FOUT
    FILE *fp;
    fp=fopen("c:\\temp\\Mattoni.txt", "w");
#endif // FOUT

    // Ciclo sulla faccia più piccola
    // Primitiva è il puntatore alla tripletta che indica la faccia più piccola
    T3=0;
    while(T3<TPP_Index) {
        T2=0;
        while(T2<TPP_Index-1) {
            T1=0;
            while(T1<TPP_Index-2) {
                if ((unsigned long int)TPP[T1][1]*(unsigned long int)TPP[T2][0]*(unsigned long int)TPP[T3][1] == (unsigned long int)TPP[T1][0]*(unsigned long int)TPP[T2][1]*(unsigned long int)TPP[T3][0]) {
                    m3=1;
                    SPIGA = 0;
                    SPIGB = 0;
                    while(m3<=MoltiplicatoreMax) {
                        SPIGA += TPP[T3][1]; // TPP[T3][1] * m3
                        SPIGB += TPP[T3][0]; // TPP[T3][0] * m3
                        if (SPIGA % TPP[T1][0] ==0 && SPIGB % TPP[T2][0]==0) {
                            m2 = SPIGB / TPP[T2][0]; // m2 = B/a2
                            m1 = SPIGA / TPP[T1][0]; // m1 = A/a1
                            // Controllo che il parallelepipedo perfetto appena trovato non sia una mera scalatura
                            // se mcd(m0, m1, m2) == 1 allora è un parallelepipedo primo
                            if(MCD(m1, MCD(m2, m3))==1) {
                                ++ContaTerne;

                                SPIGC = TPP[T2][1]*m2;

                                end_t = clock()-start_t;
                                Sec = (unsigned short) (end_t / CLOCKS_PER_SEC);
                                HH = Sec/3600;
                                MM = (Sec%3600)/60;
                                Sec= (Sec%3600)%60;
                                Milli = ((end_t % CLOCKS_PER_SEC) * 1000) / CLOCKS_PER_SEC;
                                printf("M%2d - %03u:%02u:%u.%03u (%lu ticks)\n", ContaTerne, HH, MM, Sec, Milli, end_t);

                                printf("%8u %8u %8u %4u\n", TPP[T1][0], TPP[T1][1], TPP[T1][2], m1);
                                printf("%8u %8u %8u %4u\n", TPP[T2][0], TPP[T2][1], TPP[T2][2], m2);
                                printf("%8u %8u %8u %4u\n", TPP[T3][0], TPP[T3][1], TPP[T3][2], m3);
                                printf("%8lu %8lu %8lu\n", SPIGA, SPIGB, SPIGC);

#ifdef FOUT
                                fprintf(fp,"M%2d - %03u:%02u:%u.%03u (%lu ticks)\nPremere Invio per iniziare la ricerca\n", ContaTerne, HH, MM, Sec, Milli, end_t);
                                fprintf(fp,"%8u %8u %8u %4u\n", TPP[T1][0], TPP[T1][1], TPP[T1][2], m1);
                                fprintf(fp,"%8u %8u %8u %4u\n", TPP[T2][0], TPP[T2][1], TPP[T2][2], m2);
                                fprintf(fp,"%8u %8u %8u %4u\n", TPP[T3][0], TPP[T3][1], TPP[T3][2], m3);
                                fprintf(fp,"%8lu %8lu %8lu\n", SPIGA, SPIGB, SPIGC);
#endif // FOUT
                                printf("-------------------------------------\n");
                                m3 =MMAX; // forzo uscita ... forse col Break fa prima
                            }
                        }
                        ++m3;
                    }
                }
                ++T1;
            }
            ++T2;
        }
        ++T3;
    }
#ifdef FOUT

    fclose(fp);
#endif // FOUT
}



int main() {


    MT.Mx[0][0] = 1;
    MT.Mx[0][1] = 2;
    MT.Mx[0][2] = 2;
    MT.Mx[1][0] = 2;
    MT.Mx[1][1] = 1;
    MT.Mx[1][2] = 2;
    MT.Mx[2][0] = 2;
    MT.Mx[2][1] = 2;
    MT.Mx[2][2] = 3;

    start_t = clock();
    Initialize();

    end_t = clock()-start_t;
    Sec = (unsigned int) (end_t / CLOCKS_PER_SEC);
    HH = Sec/3600;
    MM = (Sec%3600)/60;
    Sec= (Sec%3600)%60;
    Milli = ((end_t % CLOCKS_PER_SEC) * 1000) / CLOCKS_PER_SEC;
    printf("\nCalcolo TPP\nTPP Totali:%d\n", TPP_Index + 1);
    printf("HHH:MM:ss.mm - %03u:%02u:%u.%03u (%lu ticks)\n", HH, MM, Sec, Milli, end_t);



    start_t = clock();
    WalkOnList();
    end_t = clock()-start_t;
    Sec = (unsigned int) (end_t / CLOCKS_PER_SEC);
    HH = Sec/3600;
    MM = (Sec%3600)/60;
    Sec= (Sec%3600)%60;
    Milli = ((end_t % CLOCKS_PER_SEC) * 1000) / CLOCKS_PER_SEC;
    printf("Fine HHH:MM:ss.mm - %03u:%02u:%u.%03u (%lu ticks)\n", HH, MM, Sec, Milli, end_t);

    return 0;
}