Drsnej kód

[Blaazen]

Poradit nepotřebuju, je to taková zajímavost.

Našel jsem tenhle kód v knize Practical Analog and Digital Filter Design (PADFD.pdf), autor Les Thede, 2004 (i když on sám autorem kódu nebude, spíš nějaká třetí strana). V knize je to jako Appendix D.

Funkce sama počítá skvěle, pro vstupní parametry 0-0.9999 stačí 6 iterací.

Jen doplním, že MAX_TERMS (maximum iterací) je 100 a ERR_SMALL=1e-15.

Kód: [Vybrat]

/*====================================================
Ellip_Integral() - calcs complete elliptic integral
using arithmetic-geometric mean method
Prototype: void Ellip_Integral(double k);
Return:
complete elliptic integral value
Arguments: k - the modulus of the integral
====================================================*/
double Ellip_Integral(double k)
{ int
i;
/* Loop counter. */
double A[MAX_TERMS],B[MAX_TERMS],
C[MAX_TERMS]; /* Array storage values. */
/* Square the modulus as required by this method.*/
k = k * k;
/* Initialize the starting values. */
A[0] = 1;
B[0] = sqrt(1-k);
C[0] = sqrt(k);
/* Iterate until error is small enough. */
for(i = 1; i < MAX_TERMS ;i++)
{ A[i] = (A[i-1] + B[i-1])/2;
B[i] = sqrt(A[i-1]*B[i-1]);
C[i] = (A[i-1] - B[i-1])/2;
if(C[i] < ERR_SMALL)
{ break;}
}
return PI / (2 * A[i]);
}

Co vy na to?

[Reklama]

[Hynek Vychodil]

A co jako? Rychle konvergujících řad je spousta. Mimoto ten algoritmus používá sqrt což je iracionální funkce takže sama o sobě se musí počítat iterací. No a že ten algoritmus naprosto zbytečně alokuje 300 double hodnot, když mu ve skutečnosti stačí tři, o tom už bych vůbec pomlčel.

Zkus tuhle verzi:

Kód: [Vybrat]

double Ellip_Integral(double k)
{
  double a = 1.0;
  double b = sqrt(1-k*k);
  while ((a-b)/2 >= ERR_SMALL) {
    double tmp = a;
    a = (tmp+b)/2;
    b = sqrt(tmp*b);
  }
  return PI/(2.0*a);
}

Nebo ekvivalent v Erlangu, což je takový malý důkaz že imperativní programování způsobuje trvalé poškození mozku. To vysvětluje, že se někdo může rozplývat nad špatně implementovaným algoritmem.

Kód: [Vybrat]

-define(?ERR_SMALL, 1.0e-15).

ellip_integral(K) when K > = 0, K < 1 ->
   ellip(1.0, math:sqrt(1.0-K*K)).

ellip(A, B) when (A-B)/2 < ?ERR_SMALL -> math:pi()/(A+B);
ellip(A, B) -> ellip((A+B)/2, math:sqrt(A*B)).

No a nakonec bych asi skončil u něčeho takového i v C, protože každý trochu slušný kompilátor mi z toho vygeneruje stejně efektivní kód jako kdybych tam měl ten matoucí cyklus.

Kód: [Vybrat]

double ellip(double a, b) {
  if ((a-b)/2 < ERR_SMALL) {
    return PI/(a+b);
  } else {
    return ellip((a+b)/2, sqrt(a*b));
  }
}

double Ellip_Integral(double k) {
  return ellip(1.0, sqrt(1-k*k));
}


[Blaazen]

No já to sem dal právě pro pobavení, vždyť jsem napsal, že pomoc neptřebuju.

To, že si autor na začátku spočítá k^2 hned to zase odmocní ( C[0] = sqrt(k); ) a ve finále to vůbec nepoužije! 
 A tomu zbytečně dva a půl kilobajtu paměti.