1 /*
   2  * Ovo je moj pokušaj implementacije Routh-Hurwitzovog kriterija stabilnosti
   3  * linearnih vremenski nepromjenjivih sustava u programskom jeziku AEC.
   4  * U biti, to je algoritam da se pomoću matrica provjeri koliko korijena neki
   5  * polinom ima desno od imaginarne osi (dakle, s pozitivnim realnim dijelom),
   6  * bez potrebe da se oni računaju (što je u principu nemoguće za polinome
   7  * višeg reda).
   8  */
   9 
  10  #target WASI // Nećemo pozivati nikakve JavaScriptine funkcije, pa će
  11               // potprogram pisan u JavaScriptu biti jednostavniji ako AEC
  12               // compiler misli da ciljamo WASI, a ne browser.
  13 
  14 Integer16 stupanj_polinoma, jesmo_li_obrnuli_polinom;
  15 Decimal64 polinom[20], matrica[20 * 20];
  16 
  17 Function postavi_stupanj_polinoma(Integer16 n) Which Returns Nothing Does
  18   jesmo_li_obrnuli_polinom := 0;
  19   stupanj_polinoma := n;
  20 EndFunction
  21 
  22 Function postavi_polinom(Integer16 indeks, Decimal64 broj) Which Returns Nothing
  23   Does
  24     polinom[indeks] := broj;
  25 EndFunction
  26 
  27 Function signum(Decimal64 x) Which Returns Integer16 Does
  28   Return x < 0 ? -1 :
  29          x = 0 ? 0 :
  30          1;
  31 EndFunction
  32 
  33 Function provjeri_nuzne_uvjete() Which Returns Integer16 Does
  34   Integer16 i := 1;
  35   While i < stupanj_polinoma + 1 Loop
  36     If not(signum(polinom[i - 1]) = signum(polinom[i])) Then
  37       Return 0;
  38     EndIf
  39     i += 1;
  40   EndWhile
  41   Return 1;
  42 EndFunction
  43 
  44 Function obrni_polinom() Which Returns Nothing Is Declared;
  45 
  46 Function f(Integer16 i, Integer16 j) Which Returns Integer16 Does
  47   // Za pretvaranje indeksa dvodimenzionalnog polja u indeks jednodimenzionalnog
  48   // polja. Kada u svoj AEC compiler još nisam implementirao dvodimenzionalna
  49   // polja...
  50   Return 20 * i + j;
  51 EndFunction
  52 
  53 Function popuni_matricu() Which Returns Integer16 Does
  54   Integer16 broj_stupaca :=
  55               (stupanj_polinoma + 1) / 2 + mod(stupanj_polinoma + 1, 2),
  56             broj_redaka := stupanj_polinoma + 1;
  57   Integer16 i := 0;
  58   //Popunimo matricu prvo not-a-numbersima...
  59   While i < broj_redaka Loop
  60     Integer16 j := 0;
  61     While j < broj_stupaca Loop
  62       matrica[f(i, j)] := 0. / 0.;
  63       j += 1;
  64     EndWhile
  65     i += 1;
  66   EndWhile
  67   //Zatim idemo primjeniti Hurwitzov algoritam...
  68   i := 0;
  69   While i < broj_redaka Loop
  70     Integer16 j := 0;
  71     While j < broj_stupaca Loop
  72       If i = 0 Then // Prvi redak
  73         matrica[f(i, j)] := polinom[j * 2];
  74       ElseIf i = 1 Then // Drugi redak
  75         matrica[f(i, j)] := (j * 2 + 1 < stupanj_polinoma + 1) ?
  76                             polinom[j * 2 + 1] :
  77                             0;
  78       Else // Ostali reci...
  79         If matrica[f(i - 1, 0)] = 0 Then // Posebni slučajevi, kada se u prvom
  80                                          // stupcu matrice pojavi nula.
  81           If jesmo_li_obrnuli_polinom Then // Obrtanje polinoma nije "upalilo".
  82             Return 0;
  83           Else // Možda obrtanje polinoma "upali"...
  84             obrni_polinom(); // https://www.forum.hr/showpost.php?p=97955497&postcount=16
  85             jesmo_li_obrnuli_polinom := 1;
  86             Return popuni_matricu();
  87           EndIf
  88         EndIf
  89         matrica[f(i, j)] := (matrica[f(i - 1, 0)] *
  90                             (j + 1 < broj_stupaca ?
  91                              matrica[f(i - 2, j + 1)] : 0) -
  92                             (matrica[f(i - 2, 0)] *
  93                             (j + 1 < broj_stupaca ?
  94                              matrica[f(i - 1, j + 1)] : 0))) /
  95                             matrica[f(i - 1 , 0)];
  96       EndIf
  97       j += 1;
  98     EndWhile
  99     i += 1;
 100   EndWhile
 101   If matrica[f(broj_redaka - 1, 0)] = polinom[stupanj_polinoma] Then
 102     Return 1;
 103   EndIf
 104   Return 0;
 105 EndFunction
 106 
 107 Function broj_korijena_u_desnoj_poluravnini() Which Returns Integer16 Does
 108   Integer16 i := 1, brojac := 0;
 109   While i < stupanj_polinoma + 1 Loop
 110     brojac += not(signum(matrica[f(i, 0)]) = signum(matrica[f(i - 1, 0)]));
 111     i += 1;
 112   EndWhile
 113   Return brojac;
 114 EndFunction
 115 
 116 Function obrni_polinom() Which Returns Nothing Does
 117   // Ako vas zanima kako taj trik sa obrtanjem polinoma funkcionira, evo
 118   // jednostavnog objašnjenja: https://math.stackexchange.com/a/4689107/791819
 119   Decimal64 pomocni_polinom[20];
 120   Integer16 i := 0, j := stupanj_polinoma;
 121   While i < stupanj_polinoma + 1 Loop
 122     pomocni_polinom[i] := polinom[j];
 123     i += 1;
 124     j -= 1;
 125   EndWhile
 126   i := 0;
 127   While i < stupanj_polinoma + 1 Loop
 128     polinom[i] := pomocni_polinom[i];
 129     i += 1;
 130   EndWhile
 131 EndFunction
 132 
 133 Function broj_iz_matrice(Integer16 i, Integer16 j) Which Returns Decimal64 Does
 134   //Ova se funkcija poziva iz JavaScripta kad se crta tablica...
 135   Return matrica[f(i, j)];
 136 EndFunction
 137