
/* Automorphisms of a K3 surface in P^1 x P^1 x P^1	*/
/* (1+x^2)(1+y^2)(1+z^2) + Axyz = B */

#include "k3.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

static double a=A, b=B, sep=SEP;

/* Very large arrays */
point qs[PSMAX];

/* Set p.z from {x,y} so p lies on K3 */

point lift(x,y)
	double x, y;
{
	double qa, qb, qc;
	point p;

	qa = (1+x*x)*(1+y*y);
	qb = a*x*y;
	qc = qa - b;
	p.x = x;
	p.y = y;
	p.z = (-qb + sqrt(qb*qb - 4*qa*qc))/(2*qa);
	return(p);
}

/* Check if (x,y) lifts to a real point */
int unreal(x,y)
	double x,y;
{
	double qa, qb, qc;

	qa = (1+x*x)*(1+y*y);
	qb = a*x*y;
	qc = qa - b;
	return(qb*qb - 4*qa*qc < 0);
}

/* Involutions on K3 */
point ix(p)
	point p;
{
	p.x = involute(p.z,p.y,p.x);
	return(p);
}

point iy(p)
	point p;
{
	p.y = involute(p.x,p.z,p.y);
	return(p);
}

point iz(p)
	point p;
{
	p.z = involute(p.x,p.y,p.z);
	return(p);
}

/* Involution of z, fixing x and y */
/* z -> (-axy)/((1+x^2)(1+y^2)) - z */
double involute(x,y,z)
	double x,y,z;
{
	z = (-a*x*y)/((1+x*x)*(1+y*y)) - z;
	return(z);
}

/* Symmetry (x,y,z) -> (x,-y,-z) */
point sym(p)
	point p;
{
	p.y = -p.y;
	p.z = -p.z;
	return(p);
}

/* Automorphism of K3 */
point f(p)
	point p;
{
	return( sym( ix( iy( iz(p)))) );
}
	

/* Change parameters for K3 surface */
void setk3(newa,newb)
	double newa, newb;
{
	a = newa;
	b = newb;
}

void setsep(s)
	double s;
{
	sep = s;
}

/* Test if real point is on upper sheet (max z) */
/* Just check if invz makes the z coordinate smaller */
int top(p)
	point p;
{
	point q;

	q = iz(p);
	return(q.z <= p.z);
}

/* Defining equation of surface */
double surf(p)
	point p;
{
	double s, x, y, z;

	x=p.x; y=p.y; z=p.z;
	s = (1+x*x)*(1+y*y)*(1+z*z)+a*x*y*z-b;
	return(s);
}


/* Gradient/|Grad|^2 of defining equation of surface */
point grad(p)
	point p;
{
	point v;
	double g,x,y,z;

	x=p.x; y=p.y; z=p.z;

	v.x = 2*x*(1+y*y)*(1+z*z)+a*y*z;
	v.y = 2*y*(1+z*z)*(1+x*x)+a*z*x;
	v.z = 2*z*(1+x*x)*(1+y*y)+a*x*y;
	g=norm(v);
	v.x = v.x/g;
	v.y = v.y/g;
	v.z = v.z/g;
	return(v);
}

/* Follow gradient from p to point on surface */
point newton(p)
	point p;
{
	point v;
	double s;
	int i;

	for(i=0; i<NEWTONMAX; i++)
	{	s = surf(p);
		if(s < NEWTONEPS && (-s) < NEWTONEPS) return(p);
		v = grad(p);
		p.x = p.x - s*v.x*NEWTONCAUTION;
		p.y = p.y - s*v.y*NEWTONCAUTION;
		p.z = p.z - s*v.z*NEWTONCAUTION;
	}
	fprintf(stderr,"Error:  Newton iteration failed after 100 steps\n");
	exit(1);
}

/* Subdivide a segment into n segments */
void subdivide(p,q,ps,n)
	point p, q, ps[];
	int n;
{
	int i;

	for(i=1; i<n; i++)
	{	ps[i].x = (i*q.x + (n-i)*p.x)/n;
	 	ps[i].y = (i*q.y + (n-i)*p.y)/n;
	 	ps[i].z = (i*q.z + (n-i)*p.z)/n;
		ps[i] = newton(ps[i]);
	}
}

/* Norm of a point */
double norm(p)
	point p;
{
	return(p.x*p.x + p.y*p.y + p.z*p.z);
}

/* Distance between points */
double dist(p,q)
	point p, q;
{
	p.x=p.x-q.x; 
	p.y=p.y-q.y; 
	p.z=p.z-q.z;
	return(sqrt(norm(p)));
}

/* Rotate point */
point rotate(p)
	point p;
{
	point q;
	double t;
	q.x = p.y;
	q.y = p.z;
	q.z = p.x;
	return(q);
}

/* Map segment */
/* Subdivides and returns new length */
int fseg(p,q,ps)
	point p, q, ps[];
{
	int i, n;

	ps[0]=f(p);
	ps[1]=f(q);
	n=dist(ps[0],ps[1])/sep;
	if(n <= 1) return(1);
	if(n >= SUBMAX) suberr(n);
	ps[n]=ps[1];
	subdivide(p,q,ps,n);
	for(i=1; i<n; i++) ps[i] = f(ps[i]);
	return(n);
}

/* Map path */
/* n = number of segments; returns new n */
int fpath(ps,np)
	point ps[];
{
	point rs[SUBMAX];
	int i, j, nq, nr;

	qs[0]=f(ps[0]);
	nq=0;
	for(i=0; i<np; i++)
	{	nr = fseg(ps[i], ps[i+1], rs);
		if(nq+nr >= PSMAX) 
		{	patherr();
			break;
		}
		for(j=1; j<=nr; j++)
		{	nq++;
			qs[nq]=rs[j];
		}
	}
	ps[0] = qs[0];
	np = 0;
	for(i=1; i<=nq; i++)
		if(dist(ps[np],qs[i]) >= sep/2.0)
		{	np++;
			ps[np]=qs[i];
		}
	if(np == 0)
	{	ps[1]=qs[nq]; 
		np=1;
	}
	return(np);
}

/* ERROR:  Too many subdivisions */
void suberr(n)
	int n;
{
	fprintf(stderr,
		"Error: %d subdivs requested, %d avail\n",
		n,SUBMAX);
	exit(1);
}

/* ERROR:  Path too long */
void patherr()
{
	fprintf(stderr,
		"Error: max path length of %d exceeded\n",
		PSMAX);
}
