/********************************************************** * Solve complex equation of degree 3: * * a z^3 + b z^2 + c z + d = 0 * * ------------------------------------------------------- * * SAMPLE RUN: * * a = (1,0) b = (2,-10) c = (-4,1) d = (5,2) * * * * Equation az^3 + bz^2 + cz + d = 0 * * Solutions z: * * z = ( -0.4703, 0.6685) * * F(z) = ( -0.0000, 0.0000) * * z = ( 0.6018, -0.2911) * * F(z) = ( 0.0000, -0.0000) * * z = ( -2.1315, 9.6226) * * F(z) = ( -0.0000, 0.0000) * * * * J-P Moreau, Paris. * * (www.jpmoreau.fr) * *********************************************************** NOTE: z = Z - (b/3a) p = (c/a) - (b²/3a²) q = (2b^3/27a^3) + (d/a) - (bc)/3a²) Now the equation is reduced to Z^3 + pZ + q = 0 with the roots: Z[1..3] = cubroot(-q/2 + sqrt(q²/4 + p^3/27)) + cubroot(-q/2 - sqrt(q²/4 + p^3/27)) -------------------------------------------------------------------------------------*/ #include #include #define PI 4*atan(1.0) typedef struct { double x, y; } COMPLEX; COMPLEX a,b,c,d,u,v,w; // return ABS(z) double CABS(COMPLEX z) { return (sqrt(z.x*z.x + z.y*z.y)); } // z = z1 * z2 void CMUL(COMPLEX z1, COMPLEX z2, COMPLEX *z) { z->x = z1.x*z2.x - z1.y*z2.y; z->y = z1.x*z2.y + z1.y*z2.x; } // z = z1 / z2 void CDIV(COMPLEX z1, COMPLEX z2, COMPLEX *z) { double d; COMPLEX c; d = z2.x*z2.x + z2.y*z2.y; if (d<1E-15) printf(" Complex Divide by zero!\n"); else { c.x=z2.x; c.y=-z2.y; CMUL(z1,c,z); z->x=z->x/d; z->y=z->y/d; } } // z1 = SQRT(z) void CSQRT(COMPLEX z, COMPLEX *z1) { double r; r=sqrt(z.x*z.x+z.y*z.y); z1->x = sqrt((r+z.x)/2); z1->y = sqrt((r-z.x)/2); if (z.y<0) z1->y=-z1->y; } double ATAN(double numerator, double denominator) { // Return a phase between -PI and +PI double phase; if (fabs(denominator) < 1e-15) { if (numerator < 0) return (-PI/2); else return (PI/2); } else { phase=atan(numerator/denominator); if (denominator < 0) phase += PI; return phase; } } //ATAN() // z1, z2, z3 = z^(1/3) void CUBRT(COMPLEX z, COMPLEX *z1, COMPLEX *z2, COMPLEX *z3) { double r,t,tt; r=sqrt(z.x*z.x+z.y*z.y); r=pow(r,1.0/3.0); t = ATAN(z.y,z.x); t=t/3.0; z1->x=r*cos(t); z1->y=r*sin(t); tt=t-(2*PI/3); z2->x=r*cos(tt); z2->y=r*sin(tt); tt=t+(2*PI/3); z3->x=r*cos(tt); z3->y=r*sin(tt); } // z1 = a*z^3 + b*z2 + c*z + d void F(COMPLEX z, COMPLEX *z1) { CMUL(a,z,&u); CMUL(z,u,&v); CMUL(z,v,&w); z1->x=w.x; z1->y=w.y; CMUL(b,z,&u); CMUL(u,z,&v); z1->x=z1->x+v.x; z1->y=z1->y+v.y; CMUL(c,z,&u); z1->x=z1->x+u.x+d.x; z1->y=z1->y+u.y+d.y; } void main() { COMPLEX det, p, q, sdet, u1, u2, v1, v2, z1; COMPLEX z[10]; int i; // seek three complex roots of complex equation: // a z^3 + b z^2 + c z + d = 0 a.x = 1.0; a.y = 0.0; b.x = 2.0; b.y = -10.0; c.x = -4.0; c.y = 1.0; d.x = 5.0; d.y = 2.0; // calculate p = c/a - b²/3a² CDIV(c,a,&u); CMUL(b,b,&v); CMUL(a,a,&w); w.x=3.0*w.x; w.y=3.0*w.y; CDIV(v,w,&z1); p.x = u.x - z1.x; p.y = u.y - z1.y; // calculate q = 2b^3/27a^3 + d/a - bc/3a² CMUL(b,v,&w); w.x=2.0*w.x; w.y=2.0*w.y; CMUL(a,a,&u); CMUL(u,a,&v); v.x=27.0*v.x; v.y=27.0*v.y; CDIV(w,v,&q); CDIV(d,a,&u); q.x=q.x+u.x; q.y=q.y+u.y; CMUL(b,c,&u); CMUL(a,a,&v); v.x=3.0*v.x; v.y=3.0*v.y; CDIV(u,v,&w); q.x=q.x-w.x; q.y=q.y-w.y; // calculate det = q²/4 + p^3/27 CMUL(q,q,&u); u.x=u.x/4.0; u.y=u.y/4.0; CMUL(p,p,&v); CMUL(p,v,&w); w.x=w.x/27.0; w.y=w.y/27.0; det.x=u.x + w.x; det.y=u.y + w.y; CSQRT(det,&sdet); //now sdet contains sqrt(q²/4 + p^3/27) v.x = -q.x/2.0 + sdet.x; v.y = -q.y/2.0 + sdet.y; CUBRT(v,&u,&u1,&u2); // 3 cubic roots w.x = -q.x/2.0 - sdet.x; w.y = -q.y/2.0 - sdet.y; CUBRT(w,&v,&v1,&v2); // 3 cubic roots z[1].x = u.x + v.x; z[1].y = u.y + v.y; z[2].x = u.x + v1.x; z[2].y = u.y + v1.y; z[3].x = u.x + v2.x; z[3].y = u.y + v2.y; z[4].x = u1.x + v.x; z[4].y = u1.y + v.y; z[5].x = u1.x + v1.x; z[5].y = u1.y + v1.y; z[6].x = u1.x + v2.x; z[6].y = u1.y + v2.y; z[7].x = u2.x + v.x; z[7].y = u2.y + v.y; z[8].x = u2.x + v1.x; z[8].y = u2.y + v1.y; z[9].x = u2.x + v2.x; z[9].y = u2.y + v2.y; // Note: only 3 z[i] are correct solutions u.x=3.0*a.x; u.y=3.0*a.y; CDIV(b,u,&v); for (i=1; i<10; i++) { //substract b/3a z[i].x = z[i].x - v.x; z[i].y = z[i].y - v.y; } printf("\n Equation az^3 + bz^2 + cz + d = 0\n"); printf(" Solutions z:\n"); for (i=1; i<10; i++) { F(z[i],&z1); if (CABS(z1) < 0.001) { // detect if solution is correct printf(" z = (%8.4f,%8.4f)\n", z[i].x, z[i].y); printf(" F(z) = (%8.4f,%8.4f)\n", z1.x, z1.y); } } printf("\n"); } // end of file croot3.cpp