#include <ctime></ctime>

#include <cstdlib></cstdlib>

#include <cassert></cassert>

#include <complex></complex>

#include <helianthus.h></helianthus.h>

typedef double Real;

typedef std::complex<real> Complex;</real>

const Real epsilon = 1e-4;

const Real epsilonSqr = epsilon*epsilon;

inline Real magSqr(Complex x)

  { return x.real()*x.real() + x.imag()*x.imag(); }

inline bool isZero(Real x)

  { return fabs(x) < epsilon; }

inline bool isZero(Complex x)

  { return magSqr(x) <= epsilonSqr; }

inline Complex verify2(Complex r, Real a, Real b, Real c) {

  assert(isZero(a*r*r + b*r + c));

  return r;

}

inline Complex verify3(Complex r, Real a, Real b, Real c, Real d) {

  assert(isZero(a*r*r*r + b*r*r + c*r + d));

  return r;

}

inline Complex verify4(Complex r, Real a, Real b, Real c, Real d, Real e) {

  assert(isZero(a*r*r*r*r + b*r*r*r + c*r*r + d*r + e));

  return r;

}

int solve2(Complex *roots, Real a, Real b, Real c) {

  if (isZero(a)) {

    if (isZero(b)) return 0;

    roots[0] = verify2(-c/b, a, b, c);

    return 1;

  }

  Real D = b*b - a*c*4;

  Real k = Real(0.5)/a;

  Complex d = D < 0 ? Complex(0, sqrt(-D)) : Complex(sqrt(D));

  roots[0] = verify2((-b - d)*k, a, b, c);

  roots[1] = verify2((-b + d)*k, a, b, c);

  return 2;

}

int solve3(Complex *roots, Real a, Real b, Real c, Real d) {

  if (isZero(a)) return solve2(roots, b, c, d);

  // x = y - b/(3*a)

  // y*y*y + p*y + q = 0

  Real p = (3*a*c - b*b)/(3*a*a);

  Real q = (27*a*a*d - 9*a*b*c + 2*b*b*b)/(27*a*a*a);

  Real Q = p*p*p/27 + q*q/4;

  Complex Qs = Q < 0 ? Complex(0, sqrt(-Q)) : Complex(sqrt(Q));

  Complex A = pow(-q/2 + Qs, Real(1)/3);

  Complex B = pow(-q/2 - Qs, Real(1)/3);

  // choose complimentary B for A (A*B must be equal -p/3)

  Complex rot(Real(-0.5), sqrt(Real(3))/2);

  if (!isZero(A*B + p/3)) B *= rot;

  if (!isZero(A*B + p/3)) B *= rot;

  Complex Y = (A - B)*Complex(0, sqrt(Real(3)));

  Complex y0 = A + B;

  Complex y1 = (-y0 - Y)/Real(2);

  Complex y2 = (-y0 + Y)/Real(2);

  Real dd = b/(3*a);

  roots[0] = verify3(y0 - dd, a, b, c, d);

  roots[1] = verify3(y1 - dd, a, b, c, d);

  roots[2] = verify3(y2 - dd, a, b, c, d);

  return 3;

}

int solve4(Complex *roots, Real a, Real b, Real c, Real d, Real e) {

  if (isZero(a)) return solve3(roots, b, c, d, e);

  //printf("%g*x^4 + %g*x^3 + %g*x^2 + %g*x + %g = 0\n", a, b, c, d, e);

  // x = y - b/(4*a)

  // y^4 + p*y^2 + q*y + r = 0

  Real dd = b/(4*a);

  Real p = (8*a*c - 3*b*b)/(8*a*a);

  Real q = (8*a*a*d - 4*a*b*c + b*b*b)/(8*a*a*a);

  Real r = (256*a*a*a*e - 64*a*a*b*d + 16*a*b*b*c - 3*b*b*b*b)/(256*a*a*a*a);

  //printf("y^4 + %g*y^2 + %g*y + %g = 0\n", p, q, r);

  if (isZero(q)) {

    // biquadratic equation

    // y^4 + p*y^2 + r = 0

    Complex y[2];

    solve2(y, 1, p, r);

    y[0] = sqrt(y[0]);

    y[1] = sqrt(y[1]);

    roots[0] = verify4(-y[0] - dd, a, b, c, d, e);

    roots[1] = verify4( y[0] - dd, a, b, c, d, e);

    roots[2] = verify4(-y[1] - dd, a, b, c, d, e);

    roots[3] = verify4( y[1] - dd, a, b, c, d, e);

    return 4;

  }

  // solve cubic equation

  // z*z*z + (p/2)*z*z + ((p*p - 4*r)/16)*z - q*q/64 = 0

  Real pp = p/2;

  Real qq = (p*p - 4*r)/16;

  Real rr = -q*q/64;

  //printf("z^3 + %g*z^2 + %g*z + %g = 0\n", pp, qq, rr);

  Complex z[3];

  solve3(z, 1, pp, qq, rr);

  z[0] = sqrt(z[0]);

  z[1] = sqrt(z[1]);

  z[2] = sqrt(z[2]);

  // we need to find signs combination where following is valid:

  // (+-z0)*(+-z1)*(+-z2) = -q/8

  Complex zzz = z[0]*z[1]*z[2];

  //printf(

  //  "sqrt(z): (%g, %g), (%g, %g), (%g, %g), zzz: (%g, %g)\n",

  //  z[0].real(), z[0].imag(),

  //  z[1].real(), z[1].imag(),

  //  z[2].real(), z[2].imag(),

  //  zzz.real(), zzz.imag() );

  assert(isZero(zzz.imag()));

  assert(isZero(fabs(zzz.real()) - fabs(q/8)));

  if ((zzz.real() > 0) == (q > 0))

    z[0] = -z[0];

  assert(isZero(z[0]*z[1]*z[2] + q/8));

  roots[0] = verify4( z[0] - z[1] - z[2] - dd, a, b, c, d, e);

  roots[1] = verify4(-z[0] + z[1] - z[2] - dd, a, b, c, d, e);

  roots[2] = verify4(-z[0] - z[1] + z[2] - dd, a, b, c, d, e);

  roots[3] = verify4( z[0] + z[1] + z[2] - dd, a, b, c, d, e);

  return 4;

}

double px, py;

double erx, ery, ex, ey;

void generate() {

  erx = randomNumber(0, 400);

  ery = randomNumber(0, 400);

}

void solve() {

  if (isZero(erx)) {

    ex = 0;

    ey = py < fabs(ery) ? -fabs(ery)

       : py > fabs(ery) ?  fabs(ery) : py;

    return;

  }

  if (isZero(ery)) {

    ey = 0;

    ex = px < fabs(erx) ? -fabs(erx)

       : px > fabs(erx) ?  fabs(erx) : px;

    return;

  }

  double k = 1/erx;

  double x0 = px*k;

  double y0 = py*k;

  k *= ery;

  k *= k;

  double l = k - 1;

  double a = l*l;

  double b = 2*l*x0;

  double c = x0*x0 + y0*y0*k - l*l;

  double d = -b;

  double e = -x0*x0;

  double dist = INFINITY;

  Complex roots[4];

  int cnt = solve4(roots, a, b, c, d, e);

  printf("%g*x^4 + %g*x^3 + %g*x^2 + %g*x + %g = 0", a, b, c, d, e);

  for(int i = 0; i < cnt; ++i) {

    printf(", (%g, %g)", roots[i].real(), roots[i].imag());

    if (!isZero(roots[i].imag())) continue;

    double x = roots[i].real();

    double y;

    if (isZero(fabs(x) - 1)) {

      y = 0;

    } else

    if (fabs(x) < 1) {

      y = sqrt(k*(1 - x*x));

      if (y0 < 0) y = -y;

    } else {

      continue;

    }

    double dd = (x0-x)*(x0-x) + (y0-y)*(y0-y);

    if (dd < dist) { ex = x*erx; ey = y*erx; dist = dd; }

    printf(", [%g, %g, %g]", x, y, dd);

  }

  printf("\n");

  printf("dist: %g\n", dist);

  assert(dist < INFINITY);

}

void init() {

  generate();

}

void draw() {

  double w = windowGetWidth();

  double h = windowGetHeight();

  saveState();

  translate(w/2, h/2);

  const double quant = 10;

  px = mouseTransformedX();

  py = mouseTransformedY();

  if (quant) {

    px = round(px/quant)*quant;

    py = round(py/quant)*quant;

  }

  if (keyWentDown("space")) generate();

  solve();

  noFill();

  line(ex, ey, px, py);

  ellipse(-erx, -ery, 2*erx, 2*ery);

  strokeWidth(3);

  point(ex, ey);

  point(px, py);

  restoreState();

}

int main() {

  windowSetVariableFrameRate();

  windowSetResizable(TRUE);

  windowSetInit(&init);

  windowSetDraw(&draw);

  windowRun();

  return 0;

}