std::complex lent ?

Version imprimable

std::complex lent ?

Bonjour.
Pour ecrire un bout de code fractal, j'ai voulue utiliser std::complex. Seulement c'est beaucoup plus lent que d'ecrire le calcul à la main :koi:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
double zx = xf + xRatio*j;
double zy = yf + yRatio*i;
double cx = zx;
double cy = zy;
 
while(buf[j]<itermax && zx*zx+zy*zy<4)
            {
		double xnew = zx * zx - zy * zy + cx;
		zy = (2 * zx * zy )+ cy;
                zx = xnew;
                ++(buf[j]);
 
            }

Code:

1
2
3
4
5
6
7
8
 
z = c = std::complex<double>( xf + xRatio*j, yf + yRatio*i) ;	
while(buf[j]<itermax && std::abs(z)<4)
            {
                z = z*z+c;
                ++(buf[j]);
 
            }

es ce normale? je développe sous visual c++ 2008 express

30/06/2008, 23h19
yan

J'ai trouvé cela :
http://www.ddj.com/linux-open-source/199702312?pgno=1

Citation:

Application of std::complex<T> could result in some shorter code, but the computational performance, which I try to maximize, would be lost.

:cry:
y as rien a faire??
01/07/2008, 01h12
droggo

Goe,

Citation:

Envoyé par Mongaulois

y as rien a faire??

Si, coder toi-même ce qu'il faut, c'est à dire pas grand chose au total, en fait.

En évitant de faire une classe si la vitesse est primordiale, ce sera toujours ça de gagné (structure + les fonctions ad-hoc est suffisant [Je vais me faire insulter par les puristes, mais il faut savoir ce qu'on veut]). :)
01/07/2008, 08h53
yan

En faite ce qui me choque c'est que c'est 5 à 10 fois plus lent avec les complex!!!!
je trouve cela enorme!!! non?
je m'attendais à quasiment le même temp
01/07/2008, 10h44
yan

y aurais peut être une début d'explication
http://readlist.com/lists/gcc.gnu.or...help/0/12.html

ce serait à cause du test

Citation:

std::abs(z)<4

Faudra que je teste

[edit]
http://shootout.alioth.debian.org/gp...ng=gcc%20&id=3
C'est bien sensé donner de trés bon resultat
http://shootout.alioth.debian.org/gp...&lang=gpp&id=0
01/07/2008, 10h51
Sylvain Togni

Cela n'a rien à voir avec les complexes. La différence de temps viens du fait que les deux codes ne font pas la même chose : le premier compare la norme au carré, le second la valeur absolue, c-à-d fait appel à sqrt().
01/07/2008, 10h58
yan

Citation:

Envoyé par Sylvain Togni

Cela n'a rien à voir avec les complexes. La différence de temps viens du fait que les deux codes ne font pas la même chose : le premier compare la norme au carré, le second la valeur absolue, c-à-d fait appel à sqrt().

Ha oui merci, je me suis embrouillé entre abs et norm :aie:

Normalement j'avais aussi essayé avec std::norm et c'était pareil...
Je vais vérifier ce soir

Pas tant que ça... c'est pas du tout le même code !

Regardons les codes:
Ton code:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 
double zx = xf + xRatio*j;
// 1 addition, 1 multiplication
double zy = yf + yRatio*i;
// 1 addition, 1 multiplication
double cx = zx;
// 0
double cy = zy;
// 0
 
while(buf[j]<itermax && zx*zx+zy*zy<4) 
{
       // 1 indirection (pas en release), 2 comparaisons, 2 multiplications, 3 additions
       double xnew = zx * zx - zy * zy + cx;
       // 2 multiplications, 1 addition  A noter qu'en mode release, le compilo ne fait pas les 2 multiplications (déjà effectuées pour le test de la boucle)
       zy = (2 * zx * zy )+ cy;
       // 2 multiplications, 1 addition
       zx = xnew;
       // 1 affectation
       ++(buf[j]);
       // 1 incrémentation (à noter qu'en mode release cette opération est fait sur un registre qui est affecté après la boucle)
}

Résultat, disons pour 10 itérations: 52 additions, 42 multiplications, 20 comparaisons, 1 affectation avec indirection

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
z = c = std::complex<double>( xf + xRatio*j, yf + yRatio*i) ;
// idem code précédant
while(buf[j]<itermax && std::abs(z)<4)
{
      // std::abs == sqrt(std::norm(z))
      // donc: 2 comparaisons, 2 multiplications, 3 additions, 1 square-root
      z = z*z+c;
      // le compilo n'a pas fait d'optimisation ici: 
      // donc: 4 multiplications (2 par opérande du z*z), 4 additions (2 par opérande du z*z et 2 pour le +c)
      // a rajouter 1 construction temporaire d'objet, et 2 affectations
      ++(buf[j]);
      // idem code précédant
}

D'une part la condition de sortie de la boucle n'est pas la même (equivalent à <16 dans le premier code), ce qui doit induire plus d'itérations.
Même à nombre d'itération identique (10) on aura:
72 additions, 62 multiplications, 20 comparaisons, 21 affectations, 10 créations d'objet temporaire, 10 square-root

01/07/2008, 11h06
yan

Citation:

Envoyé par nicroman

Pas tant que ça... c'est pas du tout le même code !

Merci, très intéressant.
Par contre

Citation:

double xnew = zx * zx - zy * zy + cx;
// 2 multiplications, 1 addition A noter qu'en mode release, le compilo ne fait pas les 2 multiplications (déjà effectuées pour le test de la boucle)

moi je voie 2 multiplications et 2 addition (et 1 changement de signe )
ce qui ferait avec la seconde instruction :
4 multiplication et 3 addition

Citation:

z = z*z+c;
// le compilo n'a pas fait d'optimisation ici:
// donc: 4 multiplications (2 par opérande du z*z), 4 additions (2 par opérande du z*z et 2 pour le +c)
// a rajouter 1 construction temporaire d'objet, et 2 affectations

pourquoi ne serait'il pas optimiser et comment tu trouve

Citation:

4 multiplications (2 par opérande du z*z), 4 additions (2 par opérande du z*z et 2 pour le +c)
01/07/2008, 16h00
nicroman

Citation:

Envoyé par Mongaulois

moi je voie 2 multiplications et 2 addition (et 1 changement de signe )

Oups oui... mais soustraction ou addition c'est le même cout.

Pour l'optimisation je sais pas... mais en tout cas, Visual a l'air de *refaire* les squares même avec optimisation !
Ce qui fait qu'on a 2 multiplications de plus à chaque itération (sans compter le sqrt)

Quant au nombre d'opérations:

z*z+c se décompose en (z*z)+c

z.a <= z.a * z.a - z.b * z.b
z.b <= 2 * z.a * z.b

puis

z.a <= z.a + c.a
z.b <= z.b + c.b

Soit... 3 additions, 4 multiplications (oups encore un fourvoiement :( )
01/07/2008, 17h01
yan

ok. Donc on est d'accord.
Je vais retester. J'ai dû faire une fausse manip.
Pour l'astuce du compilot qui ne refera pas les carré, j'y avais jamais pensé :aie:
merci des expliquations
01/07/2008, 20h42
yan

résultat des courses pour une image 1000*1000 sur l'intervalle x/y [-2,2] et une iteration max de 1000 (le calcule est parallélisé avec QtConcurrent :mouarf:) :
1- sans complex : moyenne de 78 ms
std::norm(z) et z.real()*z.real()+z.imag()*z.imag() ne change quasiment rien
2-avec z *= z; z+=c; moyenne de 200 ms
3-avec z = z*z+c; moyenne de 421 ms

donc std::complex est bien un peu plus lent. par contre la différence entre
z = z*z+c; et z *= z;z+=c; est je trouve énorme!!!!

merci de vos réponse
02/07/2008, 15h00
nicroman
C'est vraiment bizarre....
En tout cas avec Visual, j'ai beaucoup de mal à le faire optimiser les fonctions de std::complex
C'est d'autant plus bizarre, que j'ai un code de math dans mon moteur qui fait a peu pret la même chose (c'est juste une indirection sur des fonctions vectorielles, mais ca revient au même), et qui est particulièrement optimisé lui ! Est-ce du au templating ? Aucune idée...

Toujours est-il que:

z *= z; => appel de std::complex& std::complex::operator += (const std::complex& c);
z += c; => appel de std::complex& std::complex::operator + (const std::complex& c);

si on développe l'inline:
z *= z devient:

double t = z.a * z.a - z.b * z.b;
z.b = z.a * z.b + z.b * z.a;
z.a = t;

z += c devient:

z.a += c.a;
z.b += c.b;

Et après optimisation ca ressemble à:
double t = z.a * z.a - z.b * z.b;
z.b = 2 * z.a * z.b + c.b;
z.a = t + c.a;

Soit en tout... 4 multiplications, 3 additions, 2 écritures mémoire.
Reste la mise au carré non ré-utilisée à chaque itération...
Code:

1 2 3 4 5 while (std::norm(z) < 4) { z *= z; z += c; }
A donc l'air de couter, à chaque itération:
6 multiplications, 5 additions, 1 comparaison, 2 écritures mémoire.
Code:

1 2 3 4 while (std::norm(z) < 4) { z = z * z + c; }
A l'air de couter (après correction de mon post précedent), à chaque itération:
6 multiplications, 5 additions, 1 comparaison, comme précédemment, mais aussi 2 appels constructeur, 4 écritures mémoire (au lieu de 2).
02/07/2008, 15h20
yan

Citation:

Envoyé par nicroman

C'est vraiment bizarre....
En tout cas avec Visual, j'ai beaucoup de mal à le faire optimiser les fonctions de std::complex

je suis d'accord.
Je ferais un petit main pour tester la différence sous visual et g++.
J'ai aussi testé en utilisant des pow, log .. et ca semble plutôt optimisé...

voici un code de teste :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include <iostream>
#include <fstream>
#include <complex>
#include <ctime>
 
int main(int argc, char* argv[])
{
    const int w = 1000;
    const int h = 1000;
    const int iterMax = 1000;
    const double p[] ={-2.,-2.};
    const double size[] ={4.,4.};
 
    const double xRatio =size[0]/(w-1);
    const double yRatio =size[1]/(h-1);
 
    std::complex<double> z;
    std::complex<double> c;
 
 
    int * pBuffer = new int[w*h];
 
    int t1 = clock();
    for (int i=0;i<h;++i)
    {
        for (int j=0;j<w;++j)
        {
        int id = i*w+j;
        int & refRes = pBuffer[id];
        refRes =0;
 
        z = c = std::complex<double>( p[0] + xRatio*j, p[1] + yRatio*i) ;
        while(refRes<iterMax && std::norm(z)<4)
            {
            #if 1
                z = z*z+c;
            #else
                z*=z;
                z+=c;
            #endif
                ++refRes;
            }
        if (refRes==iterMax) refRes =0;
 
        }
    }
    std::cout<< (double)(clock() - t1)/CLOCKS_PER_SEC<<std::endl;
 
 
    std::ofstream fout("image.pgm");
    fout<<"P2"<<std::endl;
    fout<<w<<" "<<h<<std::endl;
    fout<<"256"<<std::endl;
    const double max =.1*iterMax;
    for (int i=0;i<w*h;++i)
    {
        fout<<(int)(255.*(pBuffer[i]>max ? 1. : (double)pBuffer[i]/max))<<" ";
    }
 
 
 
    system("pause");
 
	return 0;
}

methode 1 : z =z*z+c
methode 2 : z =*z; z+=c;

avec mingw :
methode 1 : a peu prés 7
methode 2 : a peu prés 6.5

avec visual 2005 :
methode 1 : a peu prés 12.5
methode 2 : a peu prés 3.5

02/07/2008, 16h25
Sylvain Togni

La pénalité restante avec les std::complex vient du fait que les compilateurs ne savent pas en général mettre les petit objets dans des registres. Ce qui fait qu'à chaque calcul il y a lecture mémoire->registre, calcul registre, écriture registre->mémoire.
02/07/2008, 16h30
yan

Ca expliquerai la différence entre l'utilisation de complex ou non.
Mais une différence de fois 4 entre les deux méthodes sous visual, je suis bluffé!!! :cry:
03/07/2008, 01h47
loufoque

Tests effectués avec un Core 2 Duo T7100.

GCC 4.2 : 6.2 dans les deux cas
GCC 4.3 : 8.5 dans les deux cas (hum, régression)
ICC 10.1 : 10.6 dans cas 1, 3.7 dans cas 2.

Options de GCC : -O3 -fomit-frame-pointer -march=native
Options de ICC : -fast -march=core2 -fomit-frame-pointer

Avec ajout de la methode sans complex

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#include <iostream>
#include <fstream>
#include <complex>
#include <ctime>
 
int main(int argc, char* argv[])
{
    const int w = 1000;
    const int h = 1000;
    const int iterMax = 1000;
    const double p[] ={-2.,-2.};
    const double size[] ={4.,4.};
 
    const double xRatio =size[0]/(w-1);
    const double yRatio =size[1]/(h-1);
 
    std::complex<double> z;
    std::complex<double> c;
 
 
    int * pBuffer = new int[w*h];
 
    int t1 = clock();
    for (int i=0;i<h;++i)
    {
        for (int j=0;j<w;++j)
        {
        int id = i*w+j;
        int & refRes = pBuffer[id];
        refRes =0;
 
        #if 0
            z = c = std::complex<double>( p[0] + xRatio*j, p[1] + yRatio*i) ;	
            while(refRes<iterMax && std::norm(z)<4)
                {
                #if 0
                    z = z*z+c;
                #else
                    z*=z;
                    z+=c;   
                #endif
                    ++refRes;
                }
        #else
            double zx = p[0] + xRatio*j;
            double zy = p[1] + yRatio*i;
            double cx = zx;
            double cy = zy;
            double xnew;
 
            while(refRes<iterMax && zx*zx+zy*zy<4)
            {
                xnew = zx * zx - zy * zy + cx;
		        zy = (2 * zx * zy )+ cy;
                zx = xnew;
                ++refRes;
            }
 
        #endif
 
 
        if (refRes==iterMax) refRes =0;
 
        }
    }
    std::cout<< (double)(clock() - t1)/CLOCKS_PER_SEC<<std::endl;;
 
 
    std::ofstream fout("image.pgm");
    fout<<"P2"<<std::endl;
    fout<<w<<" "<<h<<std::endl;
    fout<<"256"<<std::endl;
    const double max =.1*iterMax;
    for (int i=0;i<w*h;++i)
    {
        fout<<(int)(255.*(pBuffer[i]>max ? 1. : (double)pBuffer[i]/max))<<" ";
    }
 
 
 
    system("pause");
 
 
 
 
	return 0;
}

methode 1 : z =z*z+c
methode 2 : z =*z; z+=c;
methode 3 : sans complex

avec mingw :
methode 1 : a peu prés 7
methode 2 : a peu prés 6.5
methode 3 : a peu 1.2

avec visual 2005 :
methode 1 : a peu prés 12.5
methode 2 : a peu prés 3.5
methode 3 : a peu prés 1.15

Ca fait tout de même de sacrée différence.
Je suis trés étonné que icc à le même comportement de visual entre la methode 1 et 2.
Mais sur des calcul plus compliqué (en utilisant des log, pow, exp ...) le std::complex doit avoir de meilleur performance je pense (j'espère).