Keskustelu: Ohjelmointikysymykset: C++: Nopea tarkka neliöjuuri

Oma menetelmäni pohjautuu Metabolixin koodivinkkiin: https://www.ohjelmointiputka.net/koodivinkit/24733-c-neliöjuuriapproksimaatio-ieee-liukuluvuilla

Idea:

Alkuperäinen luku a = 100.0f.
arvataan neliöjuureksi SG = 20.0f. (Square root Guess)
tehdään uusi arvaus: SG = (SG + a/SG)/2.0f = (20 + 5)/2 = 12,5
toistetaan edellinen monta kertaa, niin päästään lähemmäksi oikeata lukua 10.

Toteutus funktiona:

float SQRT(float a)
{
   float sg = 20.0f;
   int i;

   for(i=0;i<12;i++)
   {
      sg = (sg + a/sg)/2.0f;
   }
   return sg;
}

Tämä on tarvittavan tarkka funktio mihin tahansa laskentaan, mutta se on armottoman hidas. On olemassa hieman hullunkurisen näköinen, mutta todella nopea ratkaisu.

#define SG 20.0f

//epätarkka. ei ole järkeä käyttää ( for(i=0;i<1;i++) )
#define SQRT(a) ((a)/SG + SG) * 0.5f

//tämäkin on liian epätarkka ( for(i=0;i<2;i++) )
#define SQRT(a) (((a)/SG + SG) * 0.5f + (a)/(((a)/SG + SG) * 0.5f)) * 0.5f

Koodi kasvaa huomattavasti mitä enemmän silmukoita tehdään. Siksi tarvitaan apuohjelma, joka generoi sqrt.h:n

#include <fstream.h>
#include <string>

using namespace std;

//SG = Square root guess.

int main()
{
    int i;
    ofstream o("sqrt.h");
    string str = "((a)/SG+SG)*0.5f";

    for(i=0;i<12;i++)
    {
        str = "(" + str + "+(a)/(" + str + "))*0.5f";
    }
    str = "#define SQRT(a) (" + str + ")";

    o << "#define SG 20.0f\n";
    for(i=0;i<str.length();i++)
    {
        o << str[i];
        if(i%100 == 0 && i>0) o << "\\\n";
    }

    o.close();

    return 0;
}

Tämä tekee saman, mitä yllä oleva funktio. Aikaa menee noin 0. Tarkkuus on loistava. Pyöritin tyhjää silmukkaa ja silmukkaa, joka laskee neliöjuuren, miljardi kertaa, eikä aikaeroa ollut havaittavissa.

Minun generoimani sqrt.h 12 silmukalla: http://zux.sjr.fi/touho/sqrt.h

Hauska laskutapa, mutta onkohan tuo sittenkään kovin nopea? Nimittäin yli sadan kilotavun pituisen makron laskeminen vie paljon aikaa.

Tämä testiohjelma laskee 10000 luvun neliöjuuren SQRT-makrolla:

#include <stdio.h>
#include "sqrt.h"

int main(void) {
    int i;
    double j;
    for (i = 0; i < 10000; i++) {
        j = SQRT(i);
    }
    return 0;
}

Aikaa kuluu testikoneella yli sekunti, mikä on melko paljon. Tärkeää on, että SQRT-makrolle annetaan parametrina muuttuja. Jos kutsuttaisiin esim. SQRT(64), kääntäjä voisi optimoida koko makron pois.

Mitenköhän tarkka ja nopea makro olisi, jos siinä olisi 4 - 5 tasoa laskentaa?

Aika jätti toi makro! :D

katos katos. En itse sijottanut sitä muuttujaan. Lisäksi tuon kääntämisessä menee jonkun verran aikaa...

Testauksien jälkeen voi päätellä, ettei ole mitään järkeä käyttää tätä, kun on Metabolixin keino olemassa. :)

EDIT: onko olemassa keinoa, jolla saisi tämän makron palauttamaan juuren?

#define pika_sqrtf(vastaus, luku) \
    *(unsigned long *)&(vastaus) = ((*(unsigned long *)&(luku)) >> 1) + (127UL << 22)

Tuo sinun SQRT

real    1m24.148s
user    1m21.277s
sys     0m1.648s

Math.h sqrt

real    0m0.026s
user    0m0.020s
sys     0m0.000s

Ja testikoodi oli tällainen. Käännetty "gcc sqrt.c -o sqrt -lm".

#include <stdio.h>
#include <math.h>
#include "sqrt.h"

int main(void) {
	int i;
	double j;
	for (i = 0; i < 1000000; i++) {
		//j = SQRT(i); // kommentti pois
		//j = sqrt(i); // aina toisesta
	}
	return 0;
}

Touho kirjoitti:

Pyöritin tyhjää silmukkaa ja silmukkaa, joka laskee neliöjuuren, miljardi kertaa, eikä aikaeroa ollut havaittavissa.

Luulitko, ettei kääntäjä optimoi tuota pois?

Muutenkin, tuo ylläoleva koodi käännettynä -O2 -flägillä antaa molemmissa tapauksissa (math.h ja tuo sinun) tulokseksi noin nolla sekuntia. Tällaisissa testeissä, missä tuolla neliöjuuren arvolla ei tehdä mitään, tulokset vääristyvät helposti juuri optimoinnin takia.

Tosiaan, optimointilippu tuhoaa koko silmukan, kun on selvää, ettei se tee oikeasti mitään. Usein kuitenkin noista tulee optimoinnin ansiosta paljon nopeampia kuin ilman, joten jälkimmäinenkään testi ei välttämättä ole aivan oikea. Yksi aika realistinen testaustapa on tehdä tyhjä funktio, jolle voi antaa liukuluvun parametriksi. Kun tämän sijoittaa eri kooditiedostoon, kääntäjä ei voi optimoida kutsua pois.

extern void funktio(double d);
for (i = MAKSIMI; i; --i) {
  double d = i + 0.5;
  // 1. silmukassa: funktio(d);
  // 2. silmukassa: funktio(sqrt(d));
  // 3. silmukassa: funktio(SQRT(d));
}

Otin aikaa clock-funktiolla, laitoin kaikki testit siis samaan ohjelmaan samalla kertaa.

$ gcc sqrt.c sqrt.2.c -lm -O2 -o sqrt.bin && ./sqrt.bin
määrä: 100000000
tyhjä:  0.47 s
 SQRT: 15.40 s
 sqrt:  1.78 s
(SQRT - tyhjä) / (sqrt - tyhjä) = 11.396947

Käännetyssä koodissa makron sisältö näyttää jokseenkin samalta kuin silmukan kirjoittaminen auki, nimittäin samat neljä komentoa toistuvat peräkkäin tuon 12 kertaa. Samaan lopputulokseen pääsee siis noin 200 tavun koodillakin.