+
+
+
+
+//Normal version of the Ehrlich-Aberth method
+cuDoubleComplex FirstH_EA(int i, cuDoubleComplex *Z) {
+
+ cuDoubleComplex result;
+ cuDoubleComplex C,F,Fp;
+ int j;
+ cuDoubleComplex sum;
+ cuDoubleComplex un;
+
+ //evaluate the polynomial
+ F = Fonction(Z[i]);
+ //evaluate the derivative of the poly.
+ Fp=FonctionD(Z[i]);
+
+ double mod=Cmodule(F);
+ sum.x=0;sum.y=0;
+ un.x=1;un.y=0;
+ C=Cdiv(F,Fp); //P(z)/P'(z)
+
+ //for all roots, compute the sum
+ //for the Ehrlich-Aberth iteration
+ for ( j=0 ; j<P.PolyDegre ; j++ )
+ {
+ if ( i != j )
+ {
+ sum=Cadd(sum,Cdiv(un,Csub(Z[i],Z[j])));
+ }
+ }
+ sum=Cdiv(C,Csub(un,Cmul(C,sum))); //C/(1-Csum)
+ result=Csub(Z[i], sum);
+ return (result);
+}
+
+
+//Log Exp version of the Ehrlich-Aberth method
+cuDoubleComplex NewH_EA(int i, cuDoubleComplex *Z) {
+
+ cuDoubleComplex result;
+ cuDoubleComplex F,Fp;
+ cuDoubleComplex one,denominator,sum;
+ int j;
+ one.x=1;one.y=0;
+ sum.x=0;
+ sum.y=0;
+
+ //evaluate the polynomial with
+ //the LogExp version
+ Fp = LogFonctionD(Z[i]);
+ //evaluate the derivative of the polynomial
+ //with the LogExp version
+ F = LogFonction(Z[i]);
+
+ cuDoubleComplex FdivFp=Csub(F,Fp);
+
+ //for all roots, compute the sum
+ //for the Ehrlich-Aberth iteration
+ for ( j=0 ; j<P.degrePolynome ; j++ )
+ {
+ if ( i != j )
+ {
+ sum=Cadd(sum,Cdiv(un,Csub(Z[i],Z[j])));
+ }
+ }
+
+ //then terminate the computation
+ //of the Ehrlich-Aberth method
+ denominator=Cln(Csub(un,Cexp(Cadd(FdivFp,Cln(sum)))));
+ result=Csub(FdivFp,denominator);
+ result=Csub(Z[i],Cexp(res));
+
+ return result;
+
+}
+
+
+//kernels to update a root i
+cuDoubleComplex H_EA(int i, cuDoubleComplex *Z) {
+ cuDoubleComplex c;
+ //if the root needs to be updated
+ if(!finished[i]) {
+ //according to the module of the root
+ if (Cmodule(Z[i])<=maxRadius)
+ //selects the normal version
+ c = FirstH_EA(i,Z);
+ else
+ //of the Log Exp version
+ c = NewH_EA(i,Z);
+ return c;
+ }
+ else return Z[i];
+}
+