La sfida “pilota” (si lo ammetto non ero sicuro avrebbe potuto interessare), imperniata sulla ricerca delle Triplette Pitagoriche Primitive, è stata ampiamente raccolta e retro-elaborata attraverso l’uso di vari linguaggi e su vari computer. Personalmente sono molto contento dal numero di riscontri ottenuti e ancor di più dall’altissima qualità e varietà dei lavori pervenuti.

Elaborati

La sfida (https://sys64738.org/wp-content/uploads/2019/06/00_TriplettePitagoriche.pdf) chiedeva di trovare il modo “migliore” per cercare alcune triplette pitagoriche, la richiesta era di visualizzare quelle dalla numero 981 alla 1000, secondo l’enumerazione indicata nell’articolo, comunque ognuno poi sarebbe stato libero di cercare un’altra enumerazione valida e utilizzare quella. Avendo aperto la “sfida” a qualunque retro computer e qualunque linguaggio, risulta impossibile fare una classifica che metta in fila gli elaborati in una sorta di valutazione, soprattutto se si volesse utilizzare il mero numero di secondi utilizzati per arrivare al risultato richiesto. Ho visto l’impegno profuso dagli “sfidanti” e solo per quello meritano tutti, nessuno escluso, la medaglia d’oro, il solo vedere le persone alla continua ricerca di superare se stesse rosicando piccole porzioni di secondo è stato bellissimo e appagante.

Cercherò di elencare i vari elaborati arrivati in redazione evidenziandone le peculiarità, non senza l’aiuto degli autori stessi, i quali sono i veri protagonisti di queste sfide.

Cominciamo con Antonio Porcino, il quale sin da subito ha cercato di abbattere i tempi di calcolo cercando di evitare calcoli inutili e spulciando particolarità del linguaggio di programmazione da lui scelto, il C, grazie al quale ha creato una versione per 6502 su C64 e una per Z80 su Laser500.

Inizialmente Antonio aveva creato un programma che eseguiva il lavoro nel giro di 4 secondi, ma non contento, ottimizzazione su ottimizzazione ha abbassato il tutto di un ordine di grandezza sino ad arrivare a meno di mezzo secondo nella versione per Z80.

Riporto le note di implementazione scritte direttamente da Antonio:

Sono partito da un prototipo in JavaScript per prendere confidenza col problema ed avere una versione funzionante con cui confrontarmi (triplette.js) e da cui ricavare un sorgente in C.
In JavaScript non ho fatto nessuna ottimizzazione, ho solo evitato di utilizzare la ricorsione anticipando che sarebbe stato un problema nell’ambiente a 8 bit.
Analizzando il problema ho notato che non è necessario memorizzare tutto l’albero delle soluzioni, ma solo due livelli: il livello “n” contenete le triplette di partenza, e “n+1” le triplette generate.
La mia soluzione consiste quindi nell’utilizzo di due buffer/vettori: un primo buffer di lettura contenente le triplette “origine” e un secondo buffer di scrittura per memorizzare le triplette generate.
Ad ogni iterazione, il buffer di scrittura diventa il buffer di lettura per la prossima iterazione, operazione molto veloce in C poichè basta scambiare i puntatori.
Ho utilizzato interi a 16 bit (unsigned int in C) poichè sono sufficienti a memorizzare i numeri delle triplette almeno fino al punto richiesto dal problema.
In C ho fatto una serie di piccole ottimizzazioni che hanno via via diminuito il tempo di esecuzione, rendendo però il codice sempre meno leggibile.

Alcune tecniche utilizzate:

    • variabili tutte globali

    • niente passaggio di parametri alle funzioni

 

  • la moltiplicazione delle matrici è stata velocizzata precalcolando gli elementi che sono comuni a M1, M2 e M3

  • ho scoperto che è più veloce scrivere:

  • buf[i] = valore; ++i;
    al posto di
    buf[i++] = valore;
    probabilmente ciò è dovuto al fatto che il 6502 non ha abbastanza registri per memorizzare il valore temporaneo di i prima di incrementarlo, per cui non si ha nessun guadagno al contrario di quanto si potrebbe intuire.

A questo indirizzo potrete guardare ed apprezzare il suo magnifico lavoro: https://github.com/nippur72/8-bit-projects/tree/master/c64-triplette-pitagoriche

Passiamo a Fabrizio Caruso.

Anche Fabrizio è un appassionato di C in tutte le sue forme. Ha prodotto i sorgenti per C64 e ZX Spectrum e durante questo piccolo viaggio Pitagorico è riuscito a scoprire e far correggere un bug del compilatore ZSCDD https://github.com/z88dk/z88dk/issues/1207
Oltre a estremizzare le ottimizzazioni del codice, Fabrizio si propone come “must” la portabilità del codice su vari sistemi, cosa che con il C viene permesso abbastanza agevolmente, a patto che il codice venga scritto in ANSI C, cosa che Fabrizio si ostina a fare.

 

Riporto la traduzione (vi chiedo di perdonarmi eventuali strafalcioni) della spiegazione delle ottimizzazioni da lui utilizzate (la versione originale la troverete all’interno de link in calce):

Triplette Pitagoriche Primitive

Useremo il teorema di Barning per calcolare una sequenza di Triplette Pitagoriche Primitive (TPP)
Ad ogni iterazione verranno calcolate e scritte 

nell’array di output tre nuove TPP.

Data la tripletta a,b,c in input, si vogliono calcolare le TPP a1,b1,c1; a2,b2,c2; a3,b3,c3.

Questa è una implementazione in ANSIC (senza iniezione di codice assembly). Il codice è parametrizzato opportunamente e può essere compilato con qualsiasi ANSI C compiler per qualsiasi sistema.

Ci sono 2 tipi di ottimizzazioni generiche:

  1. I prodotti matriciali sono ottimizzati sfruttando 2 strategie principali.

    1. pre-calcolo delle somme e le fattorizzazioni ricorrenti:da = a<<1;
      db = b<<1;
      dc = c<<1;
      tc = dc+c;
      s = db+tc;

    2.  
    3. Sfrutteremo le relazioni lineari tra i componenti delle nuove TPP:

      • Esempio 1:
        Possiamo scrivere:
        a2 = a1 + (db<<1);
        al posto di:
        a2 = a+db+dc;
        in quanto noi abbiamo:
        a1 = a+dc-db;

      • Esempio 2:

        Sfruttiamo il fatto che b2=b1+db.

  2. Le operazione di lettura/scrittura su array vengono ottimizzate utilizzando due puntatori (testa e coda) al posto di semplici indici. Questa tecnica risulta più veloce in quanto l’utilizzo di un indice causa un risultato scadente in compilazione in quanto produce calcoli di offset (posizionamento) a 16 bit (raddoppiando inutilmente, in questo caso, il lavoro essendo processori ad 8 bit [ndr]).

Per quanto riguarda la cpu 6502 sfrutteremo due ottimizzazioni specifiche per tale microprocessore:

  1. Imporremo al CC65 (il compilatore C utilizzato [ndr]) di memorizzare tutte le variabili globali nella pagina zero
    ciò si ottiene attraverso l’uso di:

    • direttive #pragma (vedi extern_vars.h)

    • variabili esterne (vedi extern_vars.h)

    • un file con estensione .s che permette di mappare le variabili nella pagina zero (vedi extern_c64_vars.s)

  1. Un trucco per evitare che il CC65, in alcuni casi, generi codice poco ottimizzato.
    Questa è una sorta di “Magia nera” in quanto dipende molto dalle regole di ottimizzazione interna del CC65. In questo specifico caso abbiamo la possibilità di usare un generico workaround riguardante un problema riconosciuto nel CC65:
    con i tipi dato 8 e 16 bit, CC65 potrebbe generare del codice migliore utilizzando assegnando un valore ad una variabile di comodo prima di depositarla nella posizione finale attraverso un ountatore:
    Ad Esempio, è meglio scrivere:
    c1 = da+tc-db;
    *(tail++)=c1;
    piuttosto di:
    *(tail++) = da+tc-db;

Il lavoro di Fabrizio lo si può trovare all’indirizzo: https://github.com/Fabrizio-Caruso/Primitive-Pythagorean-Triples

 

Tra i lavori pervenuti vi è anche l’ottimo lavoro di Francesco Sblendorio, il quale non contento di lavorare con qualche linguaggio ben noto quale C, Basic, Assembly ha voluto cimentarsi con il poco noto Modula-2 (nella versione Turbo) su CP/M, il successore del Pascal.

Sebbene i tempi ottenuti da Francesco siano più allungati rispetto a quelli visti fin’ora, anche il suo lavoro è un capolavoro. Nel suo algoritmo ha voluto affrontare l’albero delle soluzioni in maniera molto “ingegneristica”, ma lascio la parola allo stesso Francesco:

Perché Turbo Modula-2?

La scelta è dettata principalmente da ragioni di interesse storico: questo compilatore (mai distribuito ufficialmente da Borland) è stato scritto da Martin Odersky (l’inventore di Scala language) e non è mai stato usato da una massa critica di utenti. Pertanto è sia una scelta “sfidante” che un motivo per una ricerca filologica.

Tecnica utilizzata: Breadth Search

La soluzione consiste nell’applicazione del teorema di Barning, che consente di ottenere tutte le terne pitagoriche avendo come input una di partenza, che tipicamente è la minimale (3, 4, 5). Ciò implica l’esplorazione di un albero infinito, pertanto diventa impossibile una sua esplorazione in termini di “Depth-Search”, dato che questa tecnica presuppone un numero di nodi finito:

Viceversa, l’esplorazione è da effettuare in larghezza, tutti i nodi di un livello, prima di passare al successivo:

Ciò è stato implementato tramite una struttura a coda (queue):

PROCEDURE DoBreadthSearch();
VAR value: VECTOR;
    r: RESULTS;
    i: CARDINAL;
BEGIN
    WHILE count < 1000 DO
        INC(count);
        Pull(value);
        IF count >= 981 THEN
           INC(lines);
           Assign(output[lines], value);
        END;
        GetNext3(value, m, r);
        FOR i := 1 TO 3 DO
            Push(r[i])
        END
    END
END DoBreadthSearch;

Da notare il limite (1000) imposto nel ciclo WHILE: senza questo limite non ci sarebbe (teoricamente) termine all’algoritmo, dato che l’albero da esplorare è infinito.

NOTA: Per ottenere una Depth Search (ricerca in profondità), ammesso che l’albero da esplorare sia finito, è sufficiente utilizzare una struttura dati a pila (push-pop) anziché a coda (push-pull).

CP/M Plus (3.0) e gestione del tempo: interrupt BDOS 105

Per la rilevazione del tempo utilizzato è stato necessario (non avendo Turbo Modula-2 funzioni native in tal senso) ricorrere agli interrupt di CP/M, in particolare all’interrupt 105 della componente BDOS:

PROCEDURE GetTime(): LONGINT;
CONST GetDT = 105; (* BDOS Function *)
VAR dat: ARRAY[0..1] OF CARDINAL;
    hour, minute, second: INTEGER;
BEGIN
    BDOS(GetDT,ADR(dat));
    second := BcdToDec(IORESULT);
    minute := BcdToDec(dat[1] DIV 256);
    hour   := BcdToDec(INTEGER(BITSET(dat[1])*BITSET(65535)));
    RETURN LONG(second) + (LONG(minute)*60L) + (LONG(hour)*3600L)
END GetTime;

Tramite la chiamata BDOS viene richiamata la routine di sistema (disponibile solo nel CP/M 3.0, o “plus”, e non nelle release precedenti come la 2.2) che restituisce:

  • IORESULT, che contiene i secondi in formato BCD

  • la variabile dat, che contiene nei 16 bit meno significativi, ore e minuti, concatenati e in formato BCD (quindi “packed BCD”)

Tramite una funzione che fa uso di operatori bitwise i valori BCD vengono convertiti in decimale:

PROCEDURE BcdToDec(n: INTEGER): INTEGER;
BEGIN
    RETURN (INTEGER(BITSET(n) * BITSET(240)) DIV 16) * 10
          + INTEGER(BITSET(n) * BITSET(15))
END BcdToDec;
Risultati e prestazioni

Le prestazioni sono state misurate eseguendo il programma, sia in versione “m-code” che in codice nativo, su un emulatore di C128 in modalità CP/M
Eseguibile nativo: 114 secondi
Eseguibile m-code: 161 secondi

Questo il link per vedere il lavoro di Francesco: https://github.com/sblendorio/pythagorean-triples

 

 

Proseguiamo con Daniele Verzotto, il quale ha voluto esagerare producendo codice sia in C con iniezioni di assembly che in assembly puro, utilizzando diversi tipi di dati interi (16, 24, 32 bit), diversi microprocessori (6502, Z80) e diversi computer (C64, C128 e ZX Spectrum).

Con i suoi elaborati Daniele è riuscito ad ottenere tempi veramente mirabili, stiamo parlando di meno di mezzo secondo, stampa dei risultati compresa.

Principalmente sia nelle versioni Assembly che in C viene usato lo stesso principio: computare le nuove terne in una zona fissa della memoria (in pagina zero [notoriamente più veloce ndr]) pe poi ricopiarle nel punto di destinazione.

Il sorgente in C vengono ridotte tutte le operazioni necessarie a delle mere somme, in modo da far usare variabili locali che risultano essere più rapide. Ho scritto del codice in assembly posizionato inline, in quanto la funzione di libreria memcpy (funzione di copia della memoria da una zona ad un altra [ndr]) è clamorosamente lenta.

Nei sorgenti in assembly invece ho usato una successione di calcolo degli elementi delle terne, in modo da ottimizzare il tempo macchina impiegato, e sfruttare quanto più possibile i registri della cpu. Ad esempio, il termine 3 del primo vettore viene calcolato in mezzo al termine 2 assieme ad altri stravolgimenti dell’ordine logico che normalmente verrebbe tenuto.
Con tutta probabilità ho lasciato delle istruzioni più lente rispetto ad altre che potevano essere utilizzate.

Tutte le versioni mi sono “costate” parecchio nella parte di output a video, perché ho dovuto elaborare soluzioni e ottimizzarle, e nel caso della versione Sinclair ho dovuto anche ri-studiarmi l’organizzazione del video, che è abbastanza “originale”.

L’assembly della versione Z80 è abbastanza elementare, non conoscendone approfonditamente l’architettura, ciò non mi ha permesso di spingere a fondo le ottimizzazioni. Nella versione per Z80 inoltre non ho usato la tecnica di copiare/ricopiare dalla matrice finale per elaborare in locale; ho letto e memorizzato i dati direttamente nella loro posizione finale.”

I lavori di Daniele sono a disposizione al link che segue:

https://sys64738.org/wp-content/uploads/2019/07/TPP-Daniele_Verzotto.zip

 

Altro lavoro eccellente è quello che ci ha fatto pervenire Felice Nardella, il quale si è cimentato con un codice scritto in Basic e quindi riutilizzato con molti emulatori (C64-Sinclair QL-Atari-BBC Micro ….).

Felice ha fatto un ottimo lavoro di codifica e studio del problema ma soprattutto ha prodotto una documentazione stupenda, della quale riporterò una parte moto interessante (la documentazione e sorgente completa la potrete raggiungere tramite un link in calce alle parole di Felice) che centra in pieno lo spirito con cui vanno affrontate queste sfide:

Elaborazione di un algoritmo di ricerca delle terne pitagoriche

Come suggerito nel testo della sfida, per la ricerca delle terne, ci si è avvalsi del teorema di Barning, il quale utilizza delle moltiplicazioni matriciali, per trovare le terne pitagoriche successive, moltiplicando una terna iniziale (racchiusa in un vettore) per altre tre matrici, che sono sempre le stesse, e sono le seguenti:

Le operazioni da svolgere sono, dunque, abbastanza semplici; in definitiva, se si vuole calcolare le tre terne, che chiameremo [X2, Y2, Z2]; [X3, Y3, Z3]; [X4, Y4, Z4], successive ad una data terna, che chiameremo [A, B, C], basta eseguire la moltiplicazione matriciale di quest’ultima, con le tre matrici date M1, M2 e M3.

In particolare, schematizzando, possiamo per semplicità inserire i valori di M1 in una matrice di dimensioni 3×3, in questo modo:
A(1,1) = 1 ; A(1,2) = -2 ; A(1,3) = 2
A(2,1) = 2 ; A(2,2) = -1 ; A(2,3) = 2
A(3,1) = 2 ; A(3,2) = -2 ; A(3,3) = 3

Eseguendo il prodotto matriciale, otterremo i tre valori, che rappresentano la prima terna successiva (terna 2), che chiameremo [X2, Y2, Z2], rispetto a quella data [A, B, C] (terna 1):
X2 = A * A(1,1) + B * A(1,2) + C * A(1,3)
Y2 = A * A(2,1) + B * A(2,2) + C * A(2,3)
Z2 = A * A(3,1) + B * A(3,2) + C * A(3,3)

In effetti la matrice 3×3 non serve affatto, ma servono soltanto i valori che vi si trovano all’interno, per cui, andando a sostituire tali valori e svolgendo i calcoli, otterremo:

X2 = A * 1 + B * (-2) + C * 2 = A – 2B + C = A – B – B + C + C
Y2 = A * 2 + B * (-1) + C * 2 = 2A – B + 2C = A + A – B + C + C
Z2 = A * 2 + B * (-2) + C * 3 = 2A – 2B + 3C = A + A – B – B + C + C + C

Come si può notare, i calcoli da eseguire sono riconducibili a delle semplici addizioni e sottrazioni.
Racchiudendo il valore di X2 così ottenuto, in una variabile K, ovvero ponendo K = A – B – B + C + C otterremo:

X2 = A – B – B + C + C = K
Y2 = A + A – B + C + C = K + A + B
Z2 = A + A – B – B + C + C + C = K + A + C

Nella matrice M2, i valori sono identici alla M1 con la sola differenza che nella seconda colonna, i valori sono positivi; quindi per trovare la terna successiva [X3, Y3, Z3], svolgeremo i seguenti calcoli:
X3 = A * 1 + B * 2 + C * 2 = A + 2B + 2C = A + B + B + C + C
Y3 = A * 2 + B * 1 + C * 2 = 2A + B + 2C = A + A + B + C + C
Z3 = A * 2 + B * 2 + C * 3 = 2A + 2B + 3C = A + A + B + B + C + C + C

Ora ponendo K = K + 4 * B si avrà che:
X3 = K
Y3 = K + A – B
Z3 = K + A + C

Per l’ultima terna successiva [X4, Y4, Z4], ponendo K = K – A – A e svolgendo i calcoli, si avrà senz’altro che:
X4 = K
Y4 = K – A – B
Z4 = K – A + C

Abbiamo così trovato le tre terne successive alla prima terna data.
Volendo continuare a trovare le successive tre terne, non occorre far altro che considerare stavolta come terna data, la terna [X2, Y2, Z2] appena trovata, e rifare gli stessi calcoli sin qui esposti.
In definitiva basterà porre A = X(N), B = Y(N) e C = Z(N), ed incrementare di volta in volta la variabile N = 2 di una unità e svolgere sempre gli stessi calcoli.
A questo punto ci si potrebbe chiedere quante volte occorrerà ripetere tali operazioni, per trovare tutte le prime 1000 terne: ebbene la risposta è 333 volte! Infatti, escludendo la prima terna data e cioè A = 3, B = 4 e C = 5, le terne successive da trovare saranno 1000 – 1 = 999; quindi 999/3 = 333 volte.

I sorgenti e la documentazione completa del bellissimo lavoro di Felice potranno essere reperiti al link:

https://sys64738.org/wp-content/uploads/2019/07/TPP_Felice_Nardella.pdf

 

 

Altro sfidante altro computer, Francesco Ugga, si è cimentato in basic per MSX e NEC-PC8801.

Il listato basic di Francesco è contenuto “agiatemente” in sole 22 righe, ma potrebbe benissimo essere “compresso” in molto meno.

Se con MSX Basic Francesco ha registrato 95 secondi, è anche vero che utilizzando il NEC-PC8801 è riuscito ad arrivare sui 25.

100 time=0: rem per Nec INIT$=TIME$
105 DIM A(1000) : DIM B(1000) : DIM C(1000)
300 REM tripletta iniziale
310 I=0 :  A(I)=3 : B(I)=4 : C(I)=5
400 REM calcolo TPP
450 FOR V=1 TO 333
500 A(V+S)= A(V-1)*(1)+B(V-1)*(-2)+C(V-1)*(2)
510 B(V+S)= A(V-1)*(2)+B(V-1)*(-1)+C(V-1)*(2)
520 C(V+S)= A(V-1)*(2)+B(V-1)*(-2)+C(V-1)*(3)
530 A(1+V+S)= A(V-1)*(1)+B(V-1)*(2)+C(V-1)*(2)
540 B(1+V+S)= A(V-1)*(2)+B(V-1)*(1)+C(V-1)*(2)
550 C(1+V+S)= A(V-1)*(2)+B(V-1)*(2)+C(V-1)*(3)
560 A(2+V+S)= A(V-1)*(-1)+B(V-1)*(2)+C(V-1)*(2)
570 B(2+V+S)= A(V-1)*(-2)+B(V-1)*(1)+C(V-1)*(2)
580 C(2+V+S)= A(V-1)*(-2)+B(V-1)*(2)+C(V-1)*(3)
585 S=S+2 : NEXT V : ti=time/50: rem Per nec GO$=TIME$
590 FOR I=979 TO 997 STEP 3
600 PRINT i+1;":";A(I);B(I);C(I)
610 PRINT i+2;":";A(I+1);B(I+1);C(I+1)
620 PRINT i+3;":";A(I+2);B(I+2);C(I+2)
630 NEXT I
640 print "time: ";ti: rem Per Nec Print TIME$, GO$

Have your say