% Leta 1991 je Nick Gould objavil, da je nasel protiprimer za 30 let staro
% domnevo, da je pivotna rast pri kompletnem pivotiranju omejena z n, kjer 
% je n velikost matrike. Njegov protiprimer je bila matrika 13 x 13 s
% pivotno rastjo 13.02, ki jo je nasel s sistematicnim iskanjem takih
% matrik na superračunalniku Cray. 
%
% Ko je Alan Edelman s svojimi studenti za to matriko izracunal tocen LU
% razcep s kompletnim pivotiranjem v Mathematici, se je izkazalo, da je 
% pivotna rast le 7.35, pivotna rast 13.02 pa je posledica zaokrozitvenih 
% napak med racunanjem razcepa v dvojni natancnosti. Uspelo pa mu je
% Gouldov primer z malenkostno spremembo (en element je iz 1 spremenil v
% 1-10^(-7) tako popraviti, da tudi pri tocnem racunanju (kar je spodaj 
% simulirano z racunanjem v vecji natancnosti) pride res pivotna rast 13.02
%
% Prvi pravi protiprimer je tako objavil Alan Edelman leta 1992.

% Bor Plestenjak

% Na ta nacin racunamo v Matlabu z vecjo natancnostjo
digits(200)
L = vpa(1e18);

g = vpa([L, -L, -L, 660848918578853640, 350768677240296530, 139130936348087710, ...
    L, -L, 945463095088536990, -64358761317393848, -47259056539260776, ...
    981447528786957180, L, L, L, -L, -L, -882625441488454570, ...
    -793497892195840220, -L, -700496337540687080, L, L, -L, L, ...
    -651498589419302720, L, 493218479970826740, L, 523219868894640230, L, ...
    931478025815019150, -L, -L, -L, 906340171404097510, L, 196359942450215320, ...
    520200438016106050, -852377236166545040, L, -799595937286409320, L, ...
    -613950298735988050, -L, -L, L, L, L, L, -L, L, -641979766159483270, L, ...
    -823477739209516720, -L, L, -L, -L, L, -L, -L, -980475145622109130, ...
    L, L, -757461144210523130, 876253886818607830, -L, -L, ...
    -814104693902053870, L, L, -L, -L, L, ...
    -L, -L, L, L, L, L, 588225298469760790, L, -L, 117806934515049340, ...
    -L, L, -L, -L, -L, L, L, L, -123654398954411060, -L, -L, L, L, L, L, -L, L, ...
    L, L, -L, 167280198905618540, -L, -L, L, 670377079454039460, -L, ...
    -L, L, -L, L, -L, -L, -L, L, 734512344136362240, ...
    774209922789794840, L, L, L, L, L, L, -L, L, ...
    -L, -L, -322948030097235110, L, -L, 59471427088948606, -L, L, ...
    -773051215153670920, L, L, L, L, L, -L, -170078579523277070, L, L, ...
    -L, L, -L, -L, 918980310122519350, -L, -L, 250493402326499640, L, ...
    961431109359263460, -L, 724092990184259320, -L, L, L, -L, L, L, ...
    L, -L, -L, L]);

g = g/L;
G = reshape(g,13,13);
G = G';

% racun 1 - vzamemo Gouldovo matriko in racunamo v dvojni natancnosti
disp('Vzamemo originalno Gouldovo matriko in racunamo v dvojni natancnosti')
A1 = double(G);
[L1,U1,P1,Q1] = lucomplete(A1);
pivotna_rast1 = max(max(abs(double(U1))))

% racun 2 - vzamemo Gouldovo matriko in racunamo v visji natancnosti
disp('Vzamemo originalno Gouldovo matriko in racunamo v visji natancnosti')
A2 = G;
[L2,U2,P2,Q2] = lucomplete(A2);
pivotna_rast2 = max(max(abs(double(U2))))

% racun 3 - popravimo Gouldovo matriko in racunamo v visji natancnosti
disp('Vzamemo popravek Gouldove matrike in racunamo v visji natancnosti')
A3 = G;
A3(11,10) = vpa(1)-vpa(10^(-7));
[L3,U3,P3,Q3] = lucomplete(A3);
pivotna_rast3 = max(max(abs(double(U3))))
