
/* riemann.c */

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

typedef struct {
	Zarray4 a,b,c,d;
	Complex p1, p2, p3, p4;
	Complex th1, th2;
	int iter;
} Param;

Complex ZUNITY = { 1.0, 0.0 };

main(int argc, char *argv[])
{
	FILE *fp;
	Param par;
	Zarray4 w;
	void read_par(Param *, FILE *);
	Complex zconj(Complex);
	Zarray4 temp, u;
	Zarray4 init_w(Param);
	Zarray4 zaddm(Zarray4, Zarray4);
	Zarray4 gena(Param, int, Zarray4);
	Zarray4 genb(Param, int, Zarray4);
	Zarray4 genc(Param, int, Zarray4);
	Zarray4 gend(Param, int, Zarray4);

	if (argc == 1)
		read_par(&par, stdin);
	else if ((fp = fopen(*++argv, "r")) == NULL) {
		fprintf(stderr, "can't open %s\n", *argv);
		exit(1);
	} else {
		read_par(&par, fp);
		fclose(fp);
	}

	par.th1 = zconj(par.th1);
	par.th2 = zconj(par.th2);

	u = par.a;
	printf("a:          %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n            %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n",
		u.e[0].r, u.e[0].i, u.e[1].r, u.e[1].i, 
		u.e[2].r, u.e[2].i, u.e[3].r, u.e[3].i);
	u = par.b;
	printf("b:          %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n            %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n",
		u.e[0].r, u.e[0].i, u.e[1].r, u.e[1].i, 
		u.e[2].r, u.e[2].i, u.e[3].r, u.e[3].i);
	u = par.c;
	printf("inv(a):     %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n            %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n",
		u.e[0].r, u.e[0].i, u.e[1].r, u.e[1].i, 
		u.e[2].r, u.e[2].i, u.e[3].r, u.e[3].i);
	u = par.d;
	printf("inv(b):     %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n            %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n",
		u.e[0].r, u.e[0].i, u.e[1].r, u.e[1].i, 
		u.e[2].r, u.e[2].i, u.e[3].r, u.e[3].i);
	u.e[0] = par.p1; u.e[1] = par.p2; u.e[2] = par.th1; u.e[3] = par.th2;
	printf("p1,p2,th1,th2: %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n            %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n",
		u.e[0].r, u.e[0].i, u.e[1].r, u.e[1].i, 
		u.e[2].r, u.e[2].i, u.e[3].r, u.e[3].i);

	w = init_w(par);
	printf("\nW at start: %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n            %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n",
		w.e[0].r, w.e[0].i, w.e[1].r, w.e[1].i, 
		w.e[2].r, w.e[2].i, w.e[3].r, w.e[3].i);

	if (par.iter >= 1) {
	  temp = gena(par,1,par.a); w = zaddm(w,temp);
	  temp = genb(par,1,par.b); w = zaddm(w,temp);
	  temp = genc(par,1,par.c); w = zaddm(w,temp);
	  temp = gend(par,1,par.d); w = zaddm(w,temp);
	}

	printf("\nW (%d iter): %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n            %15.7lfJ%15.7lf %15.7lfJ%15.7lf\n",
		par.iter,w.e[0].r, w.e[0].i, w.e[1].r, w.e[1].i, 
		w.e[2].r, w.e[2].i, w.e[3].r, w.e[3].i);

	exit(0);

}

void read_par(Param *par, FILE *fp)
{
	Zarray4 zinvm(Zarray4);

        printf("Reading param...\n");

	fscanf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf\n", 
	&par->a.e[0].r, &par->a.e[0].i, &par->a.e[1].r, &par->a.e[1].i,
	&par->a.e[2].r, &par->a.e[2].i, &par->a.e[3].r, &par->a.e[3].i);

	fscanf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf\n", 
	&par->b.e[0].r, &par->b.e[0].i, &par->b.e[1].r, &par->b.e[1].i,
	&par->b.e[2].r, &par->b.e[2].i, &par->b.e[3].r, &par->b.e[3].i);

	par->c = zinvm(par->a);
	par->d = zinvm(par->b);

/*	fscanf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf\n", 
	&par->c.e[0].r, &par->c.e[0].i, &par->c.e[1].r, &par->c.e[1].i,
	&par->c.e[2].r, &par->c.e[2].i, &par->c.e[3].r, &par->c.e[3].i);

	fscanf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf\n", 
	&par->d.e[0].r, &par->d.e[0].i, &par->d.e[1].r, &par->d.e[1].i,
	&par->d.e[2].r, &par->d.e[2].i, &par->d.e[3].r, &par->d.e[3].i); */

	fscanf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf\n", 
	&par->p1.r, &par->p1.i, &par->p2.r, &par->p2.i,
	&par->p3.r, &par->p3.i, &par->p4.r, &par->p4.i);

	fscanf(fp,"%lf %lf %lf %lf\n", 
	&par->th1.r, &par->th1.i, &par->th2.r, &par->th2.i);

	fscanf(fp,"%d\n", &par->iter);

}


Zarray4 init_w(Param par)
{
	Zarray4 w; 
        Complex u, v;
	Complex zlog(Complex);
	Complex zsub(Complex, Complex);
	Complex zmult(Complex, Complex);
	Complex zdiv(Complex, Complex);

	u = zmult(par.th1,par.p1); u = zsub(ZUNITY,u);
	v = zmult(par.th1,par.p2); v = zsub(ZUNITY,v);
	u = zdiv(v,u);
	w.e[0] = zlog(u);

	u = zmult(par.th1,par.p3); u = zsub(ZUNITY,u);
	v = zmult(par.th1,par.p4); v = zsub(ZUNITY,v);
	u = zdiv(v,u);
	w.e[1] = zlog(u);

	u = zmult(par.th2,par.p1); u = zsub(ZUNITY,u);
	v = zmult(par.th2,par.p2); v = zsub(ZUNITY,v);
	u = zdiv(v,u);
	w.e[2] = zlog(u);

	u = zmult(par.th2,par.p3); u = zsub(ZUNITY,u);
	v = zmult(par.th2,par.p4); v = zsub(ZUNITY,v);
	u = zdiv(v,u);
	w.e[3] = zlog(u);

	return w;
}


Zarray4 gena(Param par, int i, Zarray4 ga)
{
	Zarray4 w, temp;
	Complex romega(Complex, Complex, Complex, Zarray4);
	Zarray4 zaddm(Zarray4, Zarray4);
	Zarray4 zmultm(Zarray4, Zarray4);
	Zarray4 gena(Param, int, Zarray4);
	Zarray4 genb(Param, int, Zarray4);
	Zarray4 gend(Param, int, Zarray4);

	w.e[0] = romega(par.p1, par.p2, par.th1, ga);
	w.e[1] = romega(par.p3, par.p4, par.th1, ga);
	w.e[2] = romega(par.p1, par.p2, par.th2, ga);
	w.e[3] = romega(par.p3, par.p4, par.th2, ga);

	if (par.iter >= ++i) {
	        temp = zmultm(par.a,ga); temp = gena(par,i,temp); 
		w = zaddm(w,temp);
	        temp = zmultm(par.b,ga); temp = genb(par,i,temp); 
		w = zaddm(w,temp);
	        temp = zmultm(par.d,ga); temp = gend(par,i,temp); 
		w = zaddm(w,temp);
	}

/*	print_w(w,i); */
	return w;
}

Zarray4 genb(Param par, int i, Zarray4 gb)
{
	Zarray4 w, temp;
	Complex romega(Complex, Complex, Complex, Zarray4);
        Zarray4 zaddm(Zarray4, Zarray4);
        Zarray4 zmultm(Zarray4, Zarray4);
        Zarray4 gena(Param, int, Zarray4);
        Zarray4 genb(Param, int, Zarray4);
        Zarray4 genc(Param, int, Zarray4);

	w.e[0] = romega(par.p1, par.p2, par.th1, gb);
	w.e[1] = romega(par.p3, par.p4, par.th1, gb);
	w.e[2] = romega(par.p1, par.p2, par.th2, gb);
	w.e[3] = romega(par.p3, par.p4, par.th2, gb);

	if (par.iter >= ++i) {
	        temp = zmultm(par.a,gb); temp = gena(par,i,temp); 
		w = zaddm(w,temp);
	        temp = zmultm(par.b,gb); temp = genb(par,i,temp); 
		w = zaddm(w,temp);
	        temp = zmultm(par.c,gb); temp = genc(par,i,temp); 
		w = zaddm(w,temp);
	}

/*	print_w(w,i); */
	return w;
}

Zarray4 genc(Param par, int i, Zarray4 gc)
{
	Zarray4 w, temp;
	Complex romega(Complex, Complex, Complex, Zarray4);
        Zarray4 zaddm(Zarray4, Zarray4);
        Zarray4 zmultm(Zarray4, Zarray4);
        Zarray4 genb(Param, int, Zarray4);
        Zarray4 genc(Param, int, Zarray4);
        Zarray4 gend(Param, int, Zarray4);

	w.e[0] = romega(par.p1, par.p2, par.th1, gc);
	w.e[1] = romega(par.p3, par.p4, par.th1, gc);
	w.e[2] = romega(par.p1, par.p2, par.th2, gc);
	w.e[3] = romega(par.p3, par.p4, par.th2, gc);

	if (par.iter >= ++i) {
	        temp = zmultm(par.b,gc); temp = genb(par,i,temp); 
		w = zaddm(w,temp);
	        temp = zmultm(par.c,gc); temp = genc(par,i,temp); 
		w = zaddm(w,temp);
	        temp = zmultm(par.d,gc); temp = gend(par,i,temp); 
		w = zaddm(w,temp);
	}

/*	print_w(w,i); */
	return w;
}

Zarray4 gend(Param par, int i, Zarray4 gd)
{
	Zarray4 w, temp;
	Complex romega(Complex, Complex, Complex, Zarray4);
        Zarray4 zaddm(Zarray4, Zarray4);
        Zarray4 zmultm(Zarray4, Zarray4);
        Zarray4 gena(Param, int, Zarray4);
        Zarray4 genc(Param, int, Zarray4);
        Zarray4 gend(Param, int, Zarray4);

	w.e[0] = romega(par.p1, par.p2, par.th1, gd);
	w.e[1] = romega(par.p3, par.p4, par.th1, gd);
	w.e[2] = romega(par.p1, par.p2, par.th2, gd);
	w.e[3] = romega(par.p3, par.p4, par.th2, gd);

	if (par.iter >= ++i) {
	        temp = zmultm(par.a,gd); temp = gena(par,i,temp); 
		w = zaddm(w,temp);
	        temp = zmultm(par.c,gd); temp = genc(par,i,temp); 
		w = zaddm(w,temp);
	        temp = zmultm(par.d,gd); temp = gend(par,i,temp); 
		w = zaddm(w,temp);
	}

/*	print_w(w,i); */
	return w;
}

Complex romega(Complex p1, Complex p2, Complex th, Zarray4 elem)
{
	Zarray2 zv1, zv2;
	Complex z1, z2;
	Complex zconj(Complex);
	Zarray2 zmultv(Zarray4, Zarray2);
	Complex zdiv(Complex, Complex);
	Complex zlog(Complex);
	Complex zmult(Complex, Complex);
	Complex zsub(Complex, Complex);

	zv1.e[0] = p1;
	zv1.e[1] = ZUNITY;
	zv2.e[0] = p2;
	zv2.e[1] = ZUNITY;

	zv1 = zmultv(elem,zv1);	z1 = zdiv(zv1.e[0],zv1.e[1]); 
	z1 = zmult(th,z1); z1 = zsub(ZUNITY,z1);

	zv2 = zmultv(elem,zv2); z2 = zdiv(zv2.e[0],zv2.e[1]); 
	z2 = zmult(th,z2); z2 = zsub(ZUNITY,z2);

	z1 = zdiv(z2,z1);

	return zlog(z1);
}
