#include <stdio.h>

/* 
  durandkernerint.c (c) Tobias Nygren 2006
  
  Finds all (complex) roots of the quartic equation

  f(x)=x^4 + a*x^3 + b*x^2 + c*x + d = 0

  to an accuracy of 2^-8 (0.4%) Using only integer arithmetics.

  Does not use the square root function.
*/


void durand_kerner(int *roots, int a, int b, int c, int d, int runs);

#define FP 8

int main(void) {
    int roots[8];
    int i;

    /* Solve f(x)=x^4 - x^3 - 6*x^2 - x + 3 = 0 */

    durand_kerner(roots, (-1)<<FP, (-6)<<FP, (-1)<<FP, (3)<<FP, 10);

    /* print result */

    for(i=0; i<8; i+=2) {
	printf("root %i: (%f, %f)\n", i>>1, (float)roots[i]/(1<<FP), (float)roots[i+1]/(1<<FP));
    }
    return 0;
}


#define dk_eval(fr, fi, r, i) \
 t=(r*r-i*i)>>FP;\
 u=(2*r*i)>>FP;\
 fr=((t*t-u*u+(a*((t*r-u*i)>>FP))+b*t+c*r)>>FP)+d;\
 fi=(2*t*u+(a*((t*i+u*r)>>FP))+b*u+c*i)>>FP;

#define dk_mult3(rr, ri, ar, ai, br, bi, cr, ci) \
 tr=(ar*br-ai*bi)>>FP;\
 ti=(ai*br+ar*bi)>>FP;\
 rr=(tr*cr-ti*ci)>>FP;\
 ri=(ti*cr+tr*ci)>>FP;\

void durand_kerner(int *roots, int a, int b, int c, int d, int runs) {
    int pr,pi,qr,qi,rr,ri,sr,si;
    int er,ei,fr,fi,gr,gi,hr,hi,ir,ii,jr,ji;
    int tr, ti;
    int t,u;
    int k;
    int mr, mi;
    int nr, ni;
    int i;

    /* magic constants */
    pr= 256; pi=  0;
    qr=  96; qi=224;
    rr=-160; ri=168;
    sr=-207; si=-77;

    for(i=0;i<runs;i++) {
	er=pr-qr; ei=pi-qi;
	fr=pr-rr; fi=pi-ri;
	gr=pr-sr; gi=pi-si;
	hr=qr-sr; hi=qi-si;
	ir=rr-sr; ii=ri-si;
	jr=qr-rr; ji=qi-ri;

	dk_mult3(mr, mi, er, ei, fr, fi, gr, gi);
	k=(mr*mr+mi*mi)>>FP;
	dk_eval(nr, ni, pr, pi);
	pr-=(nr*mr+ni*mi)/k;
	pi-=(ni*mr-nr*mi)/k;

	dk_mult3(mr, mi, -er, -ei, jr, ji, hr, hi);
	k=(mr*mr+mi*mi)>>FP;
	dk_eval(nr, ni, qr, qi);
	qr-=(nr*mr+ni*mi)/k;
	qi-=(ni*mr-nr*mi)/k;

	dk_mult3(mr, mi, -fr, -fi, -jr, -ji, ir, ii);
	k=(mr*mr+mi*mi)>>FP;
	dk_eval(nr, ni, rr, ri);
	rr-=(nr*mr+ni*mi)/k;
	ri-=(ni*mr-nr*mi)/k;

	dk_mult3(mr, mi, -gr, -gi, -hr, -hi, -ir, -ii);
	k=(mr*mr+mi*mi)>>FP;
	dk_eval(nr, ni, sr, si);
	sr-=(nr*mr+ni*mi)/k;
	si-=(ni*mr-nr*mi)/k;

    }
    *(roots++)=pr; *(roots++)=pi;
    *(roots++)=qr; *(roots++)=qi;
    *(roots++)=rr; *(roots++)=ri;
    *(roots++)=sr; *(roots++)=si;
}