momente şi schiţe de informatică şi matematică
To attain knowledge, write. To attain wisdom, rewrite.

Introducere practică în Prolog (I)

Prolog | backtracking | partiţii | triplete pitagoreice primitive
2011 jun

În programarea logică, ideea de bază este exprimată prin Algorithm = Logic + Control. Componenta Logic - asimilată cu "program Prolog", sau "bază de date" - vizează descrierea problemei: ce fapte se cunosc (aserţiuni iniţiale asupra obiectelor problemei) şi ce reguli generale ("parametrizate" faţă de obiectele problemei) trebuie folosite pentru a deduce ulterior alte fapte (noi).

Componenta Control reprezintă însuşi interpretorul (ori fiinţa raţională) care va "executa" programele; interpretorul de Prolog "aplică" metoda deducerii (plecând de la baza de date furnizată de program) şi un mecanism intern de backtracking recursiv - căutând sistematic valori cu care să "înlocuiască" necunoscutele implicate în scopul ("goal") satisfacerii clauzelor curente.

Sunt în uz diverse implementări ale componentei "Control" (vezi Prolog Systems) şi în mare ele respectă un acelaşi standard (ISO Prolog). Aici avem în vedere şi utilizăm SWI-Prolog.

O exemplificare clasică: relaţii de rudenie

În privinţa ordinii exemplificărilor cuprinse în lucrări introductive de Prolog - prezentarea unei baze de date pe o temă larg cunoscută - cum ar fi, relaţiile de rudenie - este clasică (la fel cu programele "Hello world!" cu care debutează prezentarea în alte limbaje).

Părinţii lui Fred sunt Ed şi Edna (vezi The Flintstones); Pebbles este fiica lui Fred; etc.

man('Fred').    % exprimă faptul că "Fred este bărbat" 
man('Ed').
woman('Edna').  % Edna este femeie
woman('Wilma').
woman('Pebbles').
parent('Ed', 'Fred').      % Fred este un copil al lui Ed
parent('Edna', 'Fred').
parent('Fred', 'Pebbles').
parent('Wilma', 'Pebbles').
grandparent(X, Z) :-       % X este bunic/bunică a lui Z dacă (notat :-):
    parent(X, Y)           % există Y încât este valabil parent(X, Y)
    ,                      % şi (virgula conjugă clauzele)
    parent(Y, Z).          % este satisfăcut parent(Y, Z).

Am declarat nişte fapte - man('Fred', woman('Wilma') etc. - nişte relaţii - parent(Ed,Fred) etc. - şi o regulă (pentru a fi "grandparent"). Caracterul . (denumit "full-point") încheie declaraţia.

'Fred' - cu apostrof - este o constantă (alternativa era fred, cu litere mici dar fără încadrare cu apostrof); Fred - cu majusculă iniţială şi neîncadrat de apostrof - ar fi interpretat ca o variabilă (cum sunt mai sus X, Y, Z). Interpretorul va substitui o variabilă cu un obiect sau altul, după caz.

Să lansăm interpretorul SWI-Prolog, swipl:

vb@vb:~/docere/doc/Prolog$   swipl
% library(swi_hooks) compiled into pce_swi_hooks 0.00 sec, 2,224 bytes
Welcome to SWI-Prolog (Multi-threaded, 32 bits, Version 5.10.1)
Copyright (c) 1990-2010 University of Amsterdam, VU Amsterdam
SWI-Prolog comes with ABSOLUTELY NO WARRANTY. This is free software,
and you are welcome to redistribute it under certain conditions.
Please visit http://www.swi-prolog.org for details.

For help, use ?- help(Topic). or ?- apropos(Word).

?- 

Fie "flintstones.pl" fişierul în care am salvat declaraţiile de mai sus; la promptul ?- tastăm:

?- [flintstones].
% flintstones compiled 0.00 sec, 1,604 bytes
true.

?- 

Interpretorul a încărcat şi a compilat fişierul respectiv, afişând din nou promptul. Putem pune întrebări (interogări, "queries") şi putem formula diverse scopuri (similar cu "program main()").

?- man('Fred').
true.  % acest fapt există în baza de date

?- man('Mister').
false.  % Nu există (şi nici nu poate fi dedus)

?- man(Fred).  % aici Fred este o variabilă!
Fred = 'Fred' ;  % se tastează ; pentru a obţine răspunsuri alternative
Fred = 'Ed'.

?- 

Interogarea man(Fred) - unde Fred este o variabilă (puteam folosi la fel de bine man(X), de exemplu) - are drept răspunsuri toate valorile cu care se poate înlocui "necunoscuta" Fred astfel încât să fie satisfăcută clauza "man". Interpretorul caută o potrivire în baza de date, afişează răspunsul găsit şi apoi aşteaptă: dacă se tastează ; (pentru "sau" soluţie alternativă) atunci reia căutarea.

În momentul încărcării şi compilării programului, interpretorul indexează intern termenii respectivi (după nume şi după ordinea apariţiei lor în program) - astfel că acum găseşte foarte uşor termenul "man" (împreună cu "ramurile" respective: 'Fred' şi 'Ed', în această ordine); găsind prima potrivire posibilă - anume man(Fred) = man('Fred') - interpretorul marchează intern "ramura" respectivă - astfel că (după ce se tastează ;) căutarea unei alternative va fi reluată începând din acest loc.

?- parent(X, 'Pebbles').  % Cine sunt părinţii lui 'Pebbles'?
X = 'Fred' ;  % sau 
X = 'Wilma'.

?- parent(X, 'Pebbles'), woman(X).  % Dar cine este mama lui 'Pebbles'?
X = 'Wilma'.

Pentru ultima interogare de mai sus, interpretorul găseşte întâi "potrivirea" parent('Fred', 'Pebbles') şi consideră mai departe valoarea 'Fred' pentru toate apariţiile lui X în cadrul interogării; dar cu această instanţiere a variabilei X, clauza woman('Fred') este nesatisfăcută - "Fail" - şi atunci revine - "Redo" - la cealaltă "ramură", instanţiind X cu 'Wilma':

?- trace.  % pentru a trasa pas cu pas execuţia următoarei interogări
true.

[trace]  ?- parent(X, 'Pebbles'), woman(X).
   Call: (7) parent(_G793, 'Pebbles') ? creep
   Exit: (7) parent('Fred', 'Pebbles') ? creep
   Call: (7) woman('Fred') ? creep
   Fail: (7) woman('Fred') ? creep
   Redo: (7) parent(_G793, 'Pebbles') ? creep
   Exit: (7) parent('Wilma', 'Pebbles') ? creep
   Call: (7) woman('Wilma') ? creep
   Exit: (7) woman('Wilma') ? creep
X = 'Wilma'.

Dacă valorile unuia dintre argumentele prevăzute în clauze sunt indiferente, atunci putem folosi o variabilă anonimă - reprezentată prin _:

[trace]  ?- notrace.  % Anulează trasarea pas cu pas
true.

[debug]  ?- parent(X, _). % Care sunt părinţii?
X = 'Ed' ;    % sau (se tastează ;) 
X = 'Edna' ;  % sau 
X = 'Fred' ;  % sau 
X = 'Wilma'.  

[debug]  ?- grandparent(X, 'Pebbles').  % Cine sunt bunicii lui 'Pebbles'?
X = 'Ed' ;
X = 'Edna' ;
false.

[debug]  ?- findall(X, parent(X, _), Result_list).  % Listă cu părinţii
Result_list = ['Ed', 'Edna', 'Fred', 'Wilma'].

[debug]  ?- findall(X, grandparent(X, _), L).  % Listă cu bunicii
L = ['Ed', 'Edna'].

?- explain(findall) va expune o documentare pentru predicatul predefinit respectiv (aici, findall); pentru documentare se poate folosi şi help(findall), dar "explain" indică şi fişierul-sursă aferent (uneori, aici găseşti cea mai bună lămurire a lucrurilor!).

Am observat mai sus, dar este de subliniat: căutarea în baza de date se face "de sus în jos" - ceea ce înseamnă că ordinea clauzelor este importantă - cu eventuală revenire la precedentul punct în care există alternative neexplorate; în cazul "execuţiei" unei clauze (sau interogări) compuse din mai multe subclauze conjugate prin , - ordinea operării este "de la stânga la dreapta" (iar aici "operare" nu înseamnă "evaluare şi atribuire" ca în alte limbaje, ci substituţie şi unificare, urmărind "potrivirea" finală a termenilor).

Exemplificări aritmetice

Prolog a apărut în 1972 şi are tangenţe conceptuale mai degrabă cu SQL (unde deasemenea folosim interogări şi baze de date relaţionale), decât cu limbajele "clasice" (Pascal a apărut în 1970, iar C în 1978). Pentru tratarea serioasă a unor aspecte matematice, limbajele "clasice" (începând cu limbajul de asamblare) sunt desigur cele mai potrivite; altfel însă, diverse aspecte aritmetice pot fi tratate în orice limbaj şi este interesant de văzut asemenea exemple în Prolog şi în SQL.

Cel mai mare divizor comun

Amintim cum s-ar scrie o funcţie Pascal care furnizează "Gcd":

function GCD(a, b: integer): integer;
    begin
        if b = 0 then GCD := a
        else GCD := GCD(b, a mod b);
    end;

Se recunoaşte desigur, algoritmul lui Euclid şi fiind formulat recursiv, poate fi reflectat uşor în Prolog (unde recursivitatea este "limba maternă"):

gcd(A, B, G) :-            % gcd(A,B) este G dacă:    
    (B =:= 0 -> G is A     %   în cazul când B are valoarea 0, G = A
    ;                      %   altfel,
     R is A mod B,         %     (R fiind restul împărţirii lui A la B),
     gcd(B, R, G)          %   avem G = gcd(B,R).
    ).

Am definit un predicat cu o singură clauză, gcd/3 (de aritate 3, adică având 3 argumente). Construcţia ( clauză -> clauză_1 ; clauză_2 ). corespunde cu "IF clauză THEN clauză_1 ELSE clauză_2".

Pentru cazul valorilor numerice, operatorii prevăzuţi pentru egalitate, inegalitate şi "mai mic sau egal" diferă ca simbol de alte limbaje: =:= (egalitate), =\= (inegalitate), =< (mai mic sau egal).

R is Expresie asigură evaluarea expresiei aritmetice din partea dreaptă şi eventual, instanţierea variabilei R cu valoarea rezultată; forma echivalentă este is(R, Expresie).

Lansăm swipl şi tastăm la promptul interpretorului:

?- [pitagora].    % consultă fişierul pitagora.pl 
% pitagora compiled 0.00 sec, 592 bytes
true.

?- listing.    % listează definiţiile

gcd(C, A, B) :-
	(   A=:=0
	->  B is C
	;   D is C mod A,
	    gcd(A, D, B)
	).
true.

?- gcd(12, 21, G).    % o interogare
G = 3.

listing înlocuieşte identificatorii de variabile din program cu câte o anumită majusculă.

Dar putem face şi probe mai interesante decât am făcut mai sus (pe numerele 12 şi 21):

?- X is 3**20 * 5**10 * 7**12, Y is 3**25 * 7**8 * 11**11, gcd(X, Y, G).
X = 471304534201247574228515625,      % 320 * 510 * 712 
Y = 1393590653142003761801219090073,  % 325 * 78 * 1111 
G = 20100618201669201.  % cel mai mare divizor comun 

?- G is 3**20 * 7**8.   % verificare: cel mai mare divizor comun este 320 * 78 
G = 20100618201669201.

Pentru aritmetica numerelor întregi sau raţionale, SWI-Prolog - ca şi alte interpretoare, pentru diverse limbaje - foloseşte biblioteca GMP (scrisă în C şi limbaj de asamblare) - încât putem folosi "gcd" şi pentru numere mari, cum am ilustrat mai sus.

Însă din punctul de vedere al definiţiei "gcd", proba cea mai interesantă este aceasta:

?- gcd(0, 0, G).
G = 0.

Deci se "încalcă" definiţia uzuală, în care cel mai mare divizor comun este definit cu excepţia cazului când ambii operanzi sunt nuli. Corect ar fi fost ca gcd(0, 0, G) să fie "false" (nu "G = 0"), la fel cum:

?- gcd(12, 21, 5).  % "este adevărat că 5 este cel mai mare divizor pentru 12 şi 21"?
false.              % Nu. 

?- gcd(12, 21, 3).  % "este adevărat că 3 este cel mai mare divizor pentru 12 şi 21"?
true.               % Da.

Tratarea acestui caz de excepţie (omis şi în funcţia Pascal de mai sus) ar necesita de obicei instrucţiuni în plus în cadrul funcţiei, constituţia "clasică" a unei funcţii recursive fiind: verifică "condiţia de oprire"; dacă este satisfăcută, atunci returnează un anumit rezultat şi "exit" - altfel, determină noii parametri şi reapelează funcţia cu aceştia.

Şi pe de o parte, verificarea condiţiei de oprire va decurge la fiecare nouă reapelare, iar pe de alta - "exit" înseamnă "curăţarea" succesivă a cadrelor-stivă create pe parcursul reapelărilor, revenind "înnapoi" la contextul iniţial de apel al funcţiei recursive respective.

Putem imita şi în Prolog această manieră "clasică" de tratare, dar dispunem şi de o metodă specifică, în fond mai "directă" (şi posibil, mai eficientă): anume, putem specifica diverse cazuri particulare în clauze separate (dar "în cadrul" aceluiaşi predicat):

gcd(A, 0, A) :- A > 0.     % dacă A > 0, atunci gcd(A, 0) este A

gcd(A, B, G) :-            % gcd(A,B) este G dacă:    
    B > 0,                 % B > 0 (n-a ajuns zero) şi 
    R is A mod B,          %   (R fiind restul împărţirii lui A la B),
    gcd(B, R, G).          % avem G = gcd(B,R).

Acum predicatul "gcd" este format din două clauze; condiţiile A > 0 şi B > 0 exclud cazul gcd(0, 0, G) (acum obţinem corect "false" şi nu valoarea 0). În plus, am reuşit astfel (separând clauzele) să evităm construcţia IF-Then-Else (care este în general costisitoare, ca timp de execuţie).

Lucrurile nu decurg totuşi după cum ar fi de dorit:

[trace]  ?- gcd(2, 0, 2).
   Call: (6) gcd(2, 0, 2) ? creep  % găseşte clauza gcd(A, 0, A), "înlocuind" A cu 2
^  Call: (7) 2>0 ? creep  % "A > 0" ? (din prima clauză)
^  Exit: (7) 2>0 ? creep
   Exit: (6) gcd(2, 0, 2) ? creep
true ;  % prima clauză este satisfăcută, dar există alternative (tastăm ;)
   Redo: (6) gcd(2, 0, 2) ? creep  % găseşte şi a doua clauză, cu A <-- 2 şi B <-- 0
^  Call: (7) 0>0 ? creep  % "B > 0" ? (din a doua clauză)
^  Fail: (7) 0>0 ? creep
   Fail: (6) gcd(2, 0, 2) ? creep
false.  % a doua clauză nu este satisfăcută

Având scopul gcd(2,0,2), interpretorul a căutat în baza de date curentă o declaraţie potrivită scopului respectiv şi a găsit întâi regula cu antetul gcd(A,0,A) (prin înlocuirea "necunoscutei" A cu valoarea 2, capul acestei reguli devine identic scopului de satisfăcut).

Găsirea unei "potriviri" are de obicei, două urmări: interpretorul reţine intern locul în care a găsit-o (dacă ulterior va fi cazul, va relua căutarea din acel loc) şi respectiv, se trece la "executarea" corpului constituent al regulii găsite.

În cazul redat mai sus, se execută (6) şi apoi (7), încheind cu rezultatul "true" verificarea primei reguli găsite. Dar după ce afişează acest răspuns, interpretorul trece în starea de aşteptare: dacă se tastează ; atunci interpretorul execută "Redo" - adică reia căutarea din precedentul punct în care a reţinut că există alternative (în cazul de faţă, alternativa următoare în baza de date este gcd(A,B,G) care se potriveşte scopului prin substituţia lui A cu 2, a lui B cu 0 şi a lui G cu 2).

Este clar că execuţia ar fi trebuit încheiată imediat după verificarea primei reguli. O metodă uzuală pentru a impune aceasta (influenţând din afară, mecanismul intern de backtracking) constă în folosirea operatorului denumit cut şi notat prin ! (acesta elimină - sau "taie" - punctele de revenire precedente). Prima clauză se rescrie atunci astfel:

gcd(A, 0, A) :- A > 0, !.   % ("cut") ignoră alternativele precedentului "Call"

Cu această mică modificare (am adăugat ! în prima clauză "gcd"), dacă reluăm trasarea pentru gcd(2,0,2) (redată mai sus pentru cazul fără "cut") - constatăm că acum execuţia se încheie corect, adică imediat după constatarea verificării primei clauze.

Triplete pitagoreice primitive

Să găsim triunghiurile dreptunghice distincte care au catetele numere naturale coprime; acestea sunt numite primitive Pythagorean triple (PPT).

Vorbim deci de tripletele (a,b,c) de numere întregi pozitive care satisfac relaţiile: a < b, gcd(a, b) = 1 şi c2 = b2 + a2 (unde gcd este predicatul definit mai sus, pentru cel mai mare divizor comun).

"a < b" exclude triunghiurile congruente (triplete ca (3,4,5) şi (4,3,5)), iar "a şi b sunt coprime" exclude asemănarea (de exemplu (3k, 4k, 5k) în care k este un factor arbitrar).

O descriere brută a proprietăţii PPT poate fi următoarea:

ppt(A, B, C) :-     % (A,B,C) este PPT (dar cu C <= 100) dacă:
    between(5, 100, C),  % C are o valoare 5..100
    between(4, C, B),    % B referă un întreg 4..C (valoarea curentă a lui C)
    between(3, B, A),    % A referă un întreg 3..B_curent 
    gcd(A, B, 1),        % A şi B nu au divizori comuni (afară de 1)
    C^2 =:= A^2 + B^2.   % şi în final, este verificată teorema lui Pitagora

Privind between/3: ?- help(between)., sau ?- explain(between).

Adăugăm această definiţie în fişierul "pitagora.pl" (unde există deja predicatul gcd/3) şi reîncărcăm fişierul (la promptul interpretorului tastăm "scurtătura" [pitagora].).
Formulăm întâi interogări de forma "este adevărat că tripletul cutare este PPT ?":

?- ppt(3, 4, 5).
true.  % se îndeplinesc toate clauzele PPT
?- ppt(4, 3, 5).
false.  % Nu satisface a treia clauză ("not" between(3,B,A), unde B=3 şi A=4)
?- ppt(30, 40, 50).
false.  % Nu satisface a patra clauză (avem "not" gcd(30,40, 1))

Pentru un alt exemplu, (20, 99, 101) este PPT - dar ppt(20,99,101) va răspunde "false", pentru că nu este satisfăcută prima clauză şi ca urmare, nu se mai ajunge la verificarea celorlalte clauze.

Dar definiţia de mai sus permite şi întrebări de genul "care sunt valorile care satisfac" PPT:

?- ppt(A, 20, C).
false.  % Nu există A, C încât A este prim cu 20 şi A2 + 202 = C2
?- ppt(A, 40, C).
A = 9,
C = 41 ;  % tripletul (9, 40, 41) este PPT  (şi tastăm ;)
false.    % Nu există o altă soluţie (A, 40, C)
?- ppt(X, Y, 65).
X = 33,
Y = 56 ;  % tripletul (33, 56, 65) este PPT
X = 16,
Y = 63 ;  % (16, 63, 65) este PPT
false.    % Nu mai sunt soluţii (pentru C=65)

Dacă vrem o afişare mai convenabilă, putem folosi predicatul predefinit writeln; iar dacă vrem să forţăm "Redo" (reluarea căutării pentru o altă soluţie, fără să mai tastăm ;) atunci putem adăuga fail:

... % operăm la sfârşitul definiţiei lui ppt:
    C^2 =:= A^2 + B^2,   % am şters . final, înlocuind cu ,
    writeln((A, B, C)), fail.

Reîncărcând apoi "pitagora.pl", vom obţine:

?- ppt(X, Y, 65).
33,56,65
16,63,65
false.

fail forţează backtracking-ul, încât primim deodată toate răspunsurile.

De obicei însă, alta ar fi perspectiva care ne-ar interesa: vrem tripletele "libere" (fără să fixăm vreo valoare A, B, sau C) până la o anumită limită (de exemplu, toate sistemele PPT cu C ≤ 1000; sau eventual, numai pe acelea pentru care 1000 ≤ C ≤ 2000).

Deci, în primul rând, trebuie să considerăm două argumente Min şi Max, între valorile cărora să varieze "necunoscuta" C. Putem considera pentru Min valoarea implicită 5:

ppt(N) :- ppt(_, _, N, 5). % toate PPT (c,b,a) cu 5 ≤ c ≤ N  (c > b > a) 

ppt(N, Min) :- ppt(_, _, N, Min). % toate PPT (c,b,a) cu Min ≤ c ≤ N  (c > b > a)

ppt(A, B, Max, Min) :-  
    between(Min, Max, C), % generează o valoare C încât Min ≤ C ≤ Max
    between(4, C, B),     % 4 ≤ B ≤ C_curent
    between(3, B, A),     % 3 ≤ A ≤ B_curent
    A^2 =:= (C - B) * (C + B), % ceva mai convenabil, decât C^2 =:= A^2 + B^2 
    gcd(A, B, 1),              % A şi B sunt coprime
    writeln((C, B, A)), fail.  % afişează soluţia şi reia ("Redo")

Dacă vrem toate PPT până la 100: ppt(100) (vezi imaginea alăturată). Renunţând la clauza gcd(A,B,1) (o putem "comenta", prefixând cu %) şi repetând - reîncărcăm fişierul şi lansăm ppt(100) - vom obţine toate cele 52 de triplete pitagoreice până la 100 (nu numai pe cele 16 care sunt PPT).

Însă de exemplu pentru ppt(600, 500), rezultatul (lista PPT, între 500 şi 600) se obţine mult mai încet. Aceasta, pentru că - faţă de versiunea "brută" iniţială - nu am prevăzut decât "optimizări" minore (nelegate de fondul problemei): am înlocuit "C^2 =:= A^2 + B^2" cu o expresie care se evaluează ceva mai repede şi pe de altă parte, am inversat locurile pentru această clauză şi respectiv, clauza gcd(A,B,1) (ceea ce este totuşi important: pentru exemplul din imagine, gcd(A,B,1) se va apela acum numai pentru cele 52 de triplete pitagoreice).

Un PPT (a,b,c) are câteva proprietăţi implicite foarte simple, de care am putea ţine seama în vederea eficientizării programului.

a şi b fiind coprime, rezultă că a, b, c sunt prime două câte două şi c este impar; deci between(Min, Max, C) ar trebui să genereze doar valorile C impare.

Pătratele dau sau restul 0, sau restul 1 la împărţirea prin 3; deci, reducând modulo 3 egalitatea a2 + b2 = c2, avem sau 0+1=1, sau 1+0=1 (cazul 0+0=0 trebuie respins fiindcă a,b,c sunt coprime; cazul 1+1=1 trebuie evident, respins). Rezultă că sau a, sau b este divizibil cu 3 (iar c nu poate fi multiplu de 3).

O analiză similară pentru resturile împărţirii prin 4 arată că valorile C trebuie să fie nu numai impare (cum am văzut mai sus), dar chiar congruente cu 1 modulo 4.

Analog avem: modulo 16, pătratele sunt 0, 1, 4, sau 9; analizând cazurile posibile pentru a2 + b2 = c2 (modulo 16) rezultă că fie a, fie b este multiplu de 4.

Putem evita deocamdată, modificarea "radicală" a programului de mai sus (ar trebui rescris between, încât să genereze numai anumite valori - dar în feluri diferite pentru C, B şi respectiv pentru A). Putem sintetiza cele stabilite mai sus în două clauze: produsul a*b se divide cu 12 şi respectiv, C ia numai valori congruente modulo 4 cu 1:

ppt(A, B, Max, Min) :-  
    between(Min, Max, C),  % generează o valoare C încât Min <= C <= Max
    1 =:= C mod 4,         % Dar respinge imediat ("Redo") dacă C nu este 1 modulo 4
    between(4, C, B), 
    between(3, B, A), 
    0 =:= A*B mod 12,      % "Redo" dacă A*B nu se divide cu 12
    A*A =:= (C - B) * (C + B),      % A*A este mai eficient ca A^2 (sau A**2) 
    gcd(A, B, 1),         
    writeln((C, B, A)), fail.  % afişează soluţia şi reia ("Redo")

Reîncărcând "pitagora.pl" putem constata că viteza de execuţie se măreşte de câteva ori, faţă de prima versiune. De exemplu, lista PPT până la 500 rezultă cam în 8 secunde; dar mai departe până la 1000, se ajunge totuşi la 1 minut:

?- time(ppt(1000)).    % pentru documentare, ?- help(time).
5,4,3
13,12,5
. . . . . . 
985,864,473
997,925,372

/* în cazul A^2 =:= (C - B) * (C + B)
      52,959,287 inferences, 64.073 CPU in 64.121 seconds (100% CPU, 826543 Lips) */
   
/* în cazul A*A =:= (C - B) * (C + B)   
      52,959,130 inferences, 55.258 CPU in 55.295 seconds (100% CPU, 958397 Lips) */

Iniţial, am folosit în program operatorul de ridicare la putere, A^2; dar - cum se poate deduce din listingul de mai sus - acesta necesită inferenţe în plus şi timp mai mare decât înmulţirea de întregi (fiindcă angajează reprezentări şi calcul în virgulă mobilă).

Dacă îndrăznim mai departe… ppt(2000, 1000) rezultă cam în 6.5 minute (folosind A*A; respectiv, peste 7 minute când s-ar folosi A^2).

Dar programul poate servi pentru diverse alte obiective (nu numai pentru a lista tripletele PPT între două limite fixate, ca mai sus). De exemplu: care PPT (A,B,C) cu C < 1000 au A divizibil cu 100:

?- L=[100,200,300,400,500,600,700,800,900],  % Lista multiplilor de 100 (< 1000)
|    member(A, L),          % selectează în A un membru din lista L
|    ppt(A, Y, 1000, 100).  % PPT (A, Y, C) cu C < 1000
629,621,100
641,609,200
661,589,300
689,561,400
false.

Care triplete PPT între 1000 şi 2000, au A impar între 49 şi 63:

?- between(49, 63, A),  1 =:= A mod 2,  ppt(A, Y, 2000, 1000).
1201,1200,49
1301,1300,51
1405,1404,53
1513,1512,55
1625,1624,57
1741,1740,59
1861,1860,61
1985,1984,63
false.

Pentru tripletele obţinute mai sus avem C - B = 1. Există şi alte asemenea triplete? Putem folosi "ppt" inclusiv în scopul verificării satisfacerii proprietăţii PPT prin ?- ppt(45, 1012, 1013, 5). (sau, un exemplu cu A=900: ?- ppt(900,2419,2581, 2581).).

Numărarea soluţiilor, în Prolog

Dacă ne interesează câte soluţii sunt şi nu lista acestora, în SWI-Prolog putem folosi meta-predicatul aggregate(count, Goal, Result): apelează clauzele sau predicatul Goal şi contorizează soluţiile, returnând numărul acestora în variabila "Result".

Să eliminăm în prealabil, ultima linie din programul de mai sus (nu ne mai interesează "writeln" - deci comentăm întreaga linie şi punem . după clauza gcd(A, B, 1)); reîncărcând apoi "pitagora.pl", putem obţine numărul de triplete PPT între diverse limite, astfel:

?- aggregate(count, ppt(1000), S).
S = 158.

?- aggregate(count, ppt(2000,1000), S).
S = 161.

Deci 158 de triplete PPT au ipotenuza mai mică decât 1000 - ceea ce este corect (vezi A101931) şi sunt 161 de triplete PPT cu ipotenuza cuprinsă între 1000 şi 2000.

Obs. Numărul de PPT cu ipotenuza < N este aproximativ N / (2*PI) ≊ 0.1591549*N (D. N. Lehmer, 1900) şi avem prilejul să verificăm calitatea acestei aproximări: pentru N=1000 ea dă 159, iar pentru N=2000 dă 318 (faţă de valorile exacte 158 şi respectiv 158+161 = 319).

Putem formula şi o investigaţie mai "amănunţită":

?- time((member(N, [100,200,300,400,500,600,700,800,900,1000]),
         Min is N-100,
         aggregate(count, ppt(N, Min), S),  % PPT cu ipotenuza între Min şi N
         writeln((Min-N:S)), fail)).
0-100:16    % 16 PPT cu ipotenuza < 100
100-200:16  % 16 PPT cu ipotenuza între 100 şi 200
200-300:15
300-400:16
400-500:17
500-600:15
600-700:17
700-800:16
800-900:12
900-1000:18
% % 52,959,340 inferences, 55.136 CPU in 55.187 seconds (100% CPU, 960527 Lips)

time/1 cere un singur argument şi de aceea, am parantezat setul celor patru clauze care formează interogarea (regulă uzuală, când mai multe clauze ale unui aceluiaşi scop trebuie văzute ca un singur argument). Dacă adunăm valorile obţinute aici pe intervale de lungime 100 - obţinem într-adevăr, totalul de 158 găsit anterior pentru numărul de PPT până la 1000.

Dar aggregate face parte dintr-un modul specific interpretorului SWI-Prolog. O soluţie mai portabilă constă în folosirea unui predicat ca findall (sau "find_all" în alte interpretoare):

?- findall( (A, B), ppt(A, B, 500, 5), L ), length(L, S).
L = [(3,4),(5,12),(8,15),(7,24),(20,21),(12,35),(9,40),(28,45), (..., ...)|...],
S = 80.

Orice triplet PPT este unic determinat de catetele sale; findall/3 colectează soluţiile (A,B) ale ppt(A,B,500,5) în lista L; apoi, length/2 furnizează în S lungimea listei L, deci numărul de PPT.

Desigur, am parantezat (A, B) încât să constituie un singur argument în "findall"; dar în acest caz, "parantezarea" este în acelaşi timp un şablon pentru elementele listei (în L avem: (3, 5), etc.) şi se putea folosi şi vreun alt format (de exemplu, findall(A-B, ...) ducea la L = [ 3-4, 5-12, ...]).
De observat în plus, verificarea estimării lui Lehmer: 0.1591549 * 500 ≊ 79.58 ≊ 80.

Mai dăm un exemplu (în care avem iarăşi o "parantezare", în scopul unificării într-un singur argument): să găsim PPT în care catetele sunt numere consecutive:

?- findall( A-B, (ppt(A, B, 1000, 5), B-A =:= 1), L ), length(L, S).
L = [3-4, 20-21, 119-120, 696-697],
S = 4.

Această metodă serveşte situaţia în care este necesară lista soluţiilor (ca obiect în contextul execuţiei programului, nu doar ca afişare), în scopul prelucrării în program. Dacă vrem doar numărul de soluţii, atunci metoda bazată pe "findall" devine artificială (nu este necesară lista soluţiilor).

Prin urmare, ajungem la această problemă: cum se contorizează soluţiile unui predicat, în Prolog? Într-un limbaj ca C++ instituim o variabilă globală iniţializată cu 0 şi o incrementăm după fiecare soluţie obţinută. Dar în Prolog variabilele sunt totdeauna "locale" predicatului în care apar şi în plus, nu există "atribuire" (memorare a unei valori într-o variabilă) ci "unificare"; intern, "ppt(A,B,100,5)" (şi orice alt termen) este reprezentat printr-o structură de pointeri (unul pentru functorul "ppt" şi câte unul pentru fiecare argument) şi "unificarea" revine în fond la validarea faptului că doi pointeri (sau copii locale ale lor) referă o aceeaşi zonă (deci variabilele corespondente "au aceeaşi valoare").

În Prolog, "program" înseamnă declaraţii într-o "bază de date" şi interogări - încât "global" revine la un fapt/regulă/predicat care, fie există din start (predicate predefinite) fie este definit în baza de date indicată interpretorului de către utilizator, fie este adăugat bazei de date existente, prin efectul curent al execuţiei interogărilor programului-utilizator. (este de observat că este vorba de "global", dar nici într-un caz de "variabilă globală").

Se pot adăuga declaraţii în baza de date existentă folosind assertz (sinonim cu assert) şi asserta ("la sfârşit", respectiv "la început"); iar pentru "extragere" şi eliminare - retract:

?-     assert(elev(['Popa', 'Ion'])).
true.
?-     elev(Nume).
Nume = ['Popa', 'Ion'].

?-     retract(elev(X)).
X = ['Popa', 'Ion'].
?-     elev(Nume).
false.

În acest exemplu de lucru, prin prima interogare s-a adăugat elev(['Popa', 'Ion']) în baza de date; apoi s-a chestionat: există o valoare pentru "necunoscuta" Nume care să satisfacă termenul "elev"?.
Într-o interogare ulterioară, am folosit retract pentru a extrage într-o variabilă X argumentul termenului "elev" existent şi a elimina apoi acest termen din baza de date (interogarea ulterioară "elev(Nume)" nu mai poate fi satisfăcută, fiindcă faptul respectiv nu mai există).

Un meta-predicat count(Goal, N) care să lanseze repetat un predicat Goal indicat ca argument şi să contorizeze de câte ori este satisfăcut acesta (altfel spus, să dea numărul N de soluţii ale predicatului indicat ca argument) se poate defini folosind assert şi retract astfel:

% count(Goal, N) furnizează în N numărul soluţiilor predicatului Goal
count(Goal, N):-
    assertz(contor(0)),  % iniţializează contorul de soluţii 
    count_(Goal),        % lansează repetat Goal şi actualizează contorul
    /* Următoarea clauză se va executa numai dacă count_(Goal) dă "true" */
    retract(contor(N)).  % obţine în N valoarea finală a contorului

count_(Goal):-  % Apelează repetat Goal, cât timp este satisfăcut ("true").
    call(Goal), % Când Goal dă "false", încheie count_(Goal) cu "false".
    retract(contor(N)), % Dacă Goal este "true": N = rangul soluţiei precedente,
    M is N + 1,         %   retractează contorul curent, incrementează N şi
    assertz(contor(M)), %   înscrie contorul corespunzător acestei ultime soluţii.
    fail.       % Forţează backtracking la call(Goal), pentru o nouă soluţie.        

count_(_).  % Încheie count_(Goal) cu "true" şi după ce s-au epuizat soluţiile.

Uneori e mai greu de explicat "ce face" programul, decât să-l scrii - dar ne-am străduit mai sus, să comentăm corespunzător clauzele respective. Dacă ultima clauză ar lipsi, atunci apelul "count_(Goal)" (a doua linie din "programul principal") s-ar încheia cu "false" şi nu s-ar mai ajunge la clauza "retract(contor(N))"; poate că formularea cu If-Then-Else este mai clară:

% formulare cu If-Then-Else (nu mai necesită clauza finală "count_(_).")
count(Goal, N):-
    assertz(contor(0)),  % iniţializează contorul de soluţii
    (  count_(Goal) % va apela recursiv Goal, până când rezultă "false"
    ;  % "ELSE" (count_(Goal) nu este satisfăcut) - execută "retract"-ul 
       retract(contor(N)) % 
    ).

Însă în program am preferat (am justificat undeva mai sus, de ce) varianta fără If-Then-Else, adăugând atunci clauza finală count_(_)..

Într-un exemplu redat mai sus am văzut cum putem folosi aggregate pentru a obţine numărul de soluţii: s-a generat lista soluţiilor şi s-a afişat lungimea acestei liste. Să lansăm acelaşi test, folosind de data aceasta count, pentru a număra direct soluţiile:

?- time((  % a doua paranteză este necesară fiindcă time/1 are un singur argument
          member(N,[100,200,300,400,500,600,700,800,900]),
          Min is N - 100,
          count(ppt(N, Min), S), 
          writeln(Min-N: S), fail
   )).
0-100:16
100-200:16
. . . . . . 
700-800:16
800-900:12
% 38,584,787 inferences, 39.969 CPU in 40.000 seconds (100% CPU, 965362 Lips)

Desigur, constatăm că - numărând direct soluţiile prin count, fără a genera lista acestora ca în aggregate - timpul este sensibil mai bun (40 secunde, faţă de 55 secunde).

În În jurul unei probleme de echilibru (II) ne-am ocupat de generarea unui anumit tip de partiţii întregi, constituind fişierul "espart2.pl" conţinând predicatele "espart2" şi "parts_2". Dacă vrem să contorizăm soluţiile folosind predicatul definit mai sus count, înlocuim linia

Suma_curenta =:= Numar_de_inclus -> scrie([Numar_de_inclus | Part])

din corpul predicatului parts_2" (linie prin care afişam soluţia curentă), cu:

Suma_curenta =:= Numar_de_inclus -> true  % Nu mai interesează soluţiile, ci numărul lor

şi adăugăm în fişierul "espart2.pl" următorul predicat (angajând şi count):

espart3(N, Count) :- 
    halfTrig(N, Suma_initiala),
    count(parts_2(N, [0], Suma_initiala), Count). 

Încărcând pe lângă "pitagora.pl" (în care avem definit count) şi "espart2.pl" (după salvarea modificărilor precizate mai sus), putem formula de exemplu, următoarea interogare:

?- between(10, 25, N), espart3(N, X), writeln(N: X), fail.

10:40       15:722       20:15272      25:353743
11:70       16:1314      21:28460      false.
12:124      17:2410      22:53222
13:221      18:4441      23:99820
14:397      19:8220      24:187692

Am obţinut astfel numărul de "s-partiţii" ale segmentului natural [n], pentru n=10..25 (dar vezi totuşi În jurul unei probleme de echilibru (I) o soluţie mai rapidă, angajând o funcţie generatoare pentru numărul de s-partiţii).

vezi Cărţile mele (de programare)

docerpro | Prev | Next