Routh-Hurwitz Stability Criterion

[Ravenshire]

With the help from people from forum.hr, I have modified the code of my program like this:
[hide]

Code:

Function obrni_polinom() Which Returns Nothing Is Declared;

Function f(Integer16 i, Integer16 j) Which Returns Integer16 Does
 // Za pretvaranje indeksa dvodimenzionalnog polja u indeks jednodimenzionalnog
 // polja. Kada u svoj AEC compiler još nisam implementirao dvodimenzionalna
 // polja...
 Return 20 * i + j;
EndFunction

Function popuni_matricu() Which Returns Integer16 Does
 Integer16 broj_stupaca :=
             (stupanj_polinoma + 1) / 2 + mod(stupanj_polinoma + 1, 2),
           broj_redaka := stupanj_polinoma + 1;
 Integer16 i := 0;
 //Popunimo matricu prvo not-a-numbersima...
 While i < broj_redaka Loop
   Integer16 j := 0;
   While j < broj_stupaca Loop
     matrica[f(i, j)] := 0. / 0.;
     j += 1;
   EndWhile
   i += 1;
 EndWhile
 //Zatim idemo primjeniti Hurwitzov algoritam...
 i := 0;
 While i < broj_redaka Loop
   Integer16 j := 0;
   While j < broj_stupaca Loop
     If i = 0 Then // Prvi redak
       matrica[f(i, j)] := polinom[j * 2];
     ElseIf i = 1 Then // Drugi redak
       matrica[f(i, j)] := (j * 2 + 1 < stupanj_polinoma + 1) ?
                           polinom[j * 2 + 1] :
                           0;
     Else // Ostali reci...
       If matrica[f(i - 1, 0)] = 0 Then // Posebni slučajevi, kada se u prvom
                                        // stupcu matrice pojavi nula.
         If jesmo_li_obrnuli_polinom Then // Obrtanje polinoma nije "upalilo".
           Return 0;
         Else // Možda obrtanje polinoma "upali"...
           obrni_polinom(); // https://www.forum.hr/showpost.php?p=97955497&postcount=16
           jesmo_li_obrnuli_polinom := 1;
           Return popuni_matricu();
         EndIf
       EndIf
       matrica[f(i, j)] := (matrica[f(i - 1, 0)] *
                           (j + 1 < broj_stupaca ?
                            matrica[f(i - 2, j + 1)] : 0) -
                           (matrica[f(i - 2, 0)] *
                           (j + 1 < broj_stupaca ?
                            matrica[f(i - 1, j + 1)] : 0))) /
                           matrica[f(i - 1 , 0)];
     EndIf
     j += 1;
   EndWhile
   i += 1;
 EndWhile
 If matrica[f(broj_redaka - 1, 0)] = polinom[stupanj_polinoma] Then
   Return 1;
 EndIf
 Return 0;
EndFunction

Function broj_korijena_u_desnoj_poluravnini() Which Returns Integer16 Does
 Integer16 i := 1, brojac := 0;
 While i < stupanj_polinoma + 1 Loop
   brojac += not(signum(matrica[f(i, 0)]) = signum(matrica[f(i - 1, 0)]));
   i += 1;
 EndWhile
 Return brojac;
EndFunction

Function obrni_polinom() Which Returns Nothing Does
 Decimal64 pomocni_polinom[20];
 Integer16 i := 0, j := stupanj_polinoma;
 While i < stupanj_polinoma + 1 Loop
   pomocni_polinom[i] := polinom[j];
   i += 1;
   j -= 1;
 EndWhile
 i := 0;
 While i < stupanj_polinoma + 1 Loop
   polinom[i] := pomocni_polinom[i];
   i += 1;
 EndWhile
EndFunction