               -x,-y,-z, t );
 }
 
-// Use vec2 return?
-bool quadratic0(float A, float B, float C, out float x0, out float x1) {
+int quadratic(float A, float B, float C, out vec3 x) {
    float D = B*B - 4.0*A*C;
-   if (D < 0.0) return false;
+   if (D < 0.0) return 0;
    D = sqrt(D);
    if (B < 0.0) D = -D;
-   x0 = (-B-D)/(2.0*A);
-   //x1 = (-B+D)/(2.0*A);
-   x1 = C/(A*x0);
-   if (x1 < x0) {
-     float tmp = x0; x0 = x1; x1 = tmp;
-   }
-   return true;
-}
-
-bool quadratic(float A, float B, float C, float min, inout float x) {
-  float x0,x1;
-  if (!quadratic0(A,B,C,x0,x1)) return false;
-  bool result = true;
-  if (min <= x0 && x0 < x) x = x0;
-  else if (min <= x1 && x1 < x) x = x0;
-  else result = false;
-  return result;
-}
-
-float evalcubic(float x, float a, float b, float c, float d) {
-  return d+x*(c+x*(b+x*a));
-}
-
-bool cubic(float a, float b, float c, float d, float min, inout float x) {
-  // Need a good way of handling ill-conditioned coeffs.
-  if (a == 0.0) return quadratic(b,c,d,min,x);
-  //if (abs(c) < 1e-3) highlight = true;
-  if (abs(a) < 5e-3) {
-    //highlight = true;
-#if 1
-    // This helps, but isn't perfect.
-    float A = 1.0;
-    float C = d/b;
-    float B = (c-a*C)/b;
-    return quadratic(A,B,C,min,x);
-#endif
-  }
-#if 0
-  int nroots = 1;
-  {
-    // Find out how many roots to expect by looking at
-    // the derivative. We have three roots just when
-    // there are 2 different points where the derivative
-    // goes to zero.
-    float x0, x1;
-    if (quadratic0(3.0*a,2.0*b,c,x0,x1) && x0 != x1) {
-      float y0 = evalcubic(x0,a,b,c,d);
-      float y1 = evalcubic(x1,a,b,c,d);
-      if (y0*y1 < 0.0) nroots = 3;
-    }
-  }
-#endif
+   x[0] = (-B-D)/(2.0*A);
+   x[1] = C/(A*x[0]);
+   return 2;
+}
+
+// NR algorithm for solving cubic equation
+int cubic0(float a, float b, float c, float d, out vec3 x) {
+  if (a == 0.0) return quadratic(b,c,d,x);
   float tmp = a; a = b/tmp; b = c/tmp; c = d/tmp;
   // solve x^3 + ax^2 + bx + c = 0
-  // return smallest +ve solution
   float Q = (a*a-3.0*b)/9.0;
   float R = (2.0*a*a*a - 9.0*a*b + 27.0*c)/54.0;
   float R2 = R*R, Q3 = Q*Q*Q;
-  bool found = false;
   if (R2 < Q3) {
-    //if (nroots != 3) highlight = true;
     float X = clamp(R/sqrt(Q3),-1.0,1.0);
     float theta = acos(X);
     float S = sqrt(Q); // Q must be positive since 0 <= R2 < Q3
-    float t = -2.0*S*cos(theta/3.0) - a/3.0;
-    if (t >= min && (t < x || !found)) { x = t; found = true; }
-    t = -2.0*S*cos((theta+2.0*PI)/3.0) - a/3.0;
-    if (t >= min && (t < x || !found)) { x = t; found = true; }
-    t = -2.0*S*cos((theta+4.0*PI)/3.0) - a/3.0;
-    if (t >= min && (t < x || !found)) { x = t; found = true; }
+    x[0] = -2.0*S*cos(theta/3.0) - a/3.0;
+    x[1] = -2.0*S*cos((theta+2.0*PI)/3.0) - a/3.0;
+    x[2] = -2.0*S*cos((theta+4.0*PI)/3.0) - a/3.0;
+    return 3;
   } else {
-    //if (nroots != 1) highlight = true;
     float A = -sign(R)*pow(abs(R)+sqrt(R2-Q3),0.3333);
     float B = A == 0.0 ? 0.0 : Q/A;
-    float t = (A+B) - a/3.0;
-    if (t >= min && t < x) { x = t; found = true; }
-  }
-  return found;
-}
-
-const int NLINES = 27;
-vec4 lines[NLINES*2];
-int colors[NLINES];
-
-void init() {
-  lines[0] = vec4(1,-1,0,0); lines[1] = vec4(0,0,1,0);
-  lines[2] = vec4(1,-1,0,0); lines[3] = vec4(0,0,0,1);
-  lines[4] = vec4(1,0,-1,0); lines[5] = vec4(0,1,0,0);
-  lines[6] = vec4(1,0,-1,0); lines[7] = vec4(0,0,0,1);
-  lines[8] = vec4(1,0,0,-1); lines[9] = vec4(0,1,0,0);
-  lines[10] = vec4(1,0,0,-1); lines[11] = vec4(0,0,1,0);
-
-  lines[12] = vec4(0,1,-1,0); lines[13] = vec4(1,0,0,0);
-  lines[14] = vec4(0,1,-1,0); lines[15] = vec4(0,0,0,1);
-  lines[16] = vec4(0,1,0,-1); lines[17] = vec4(1,0,0,0);
-  lines[18] = vec4(0,1,0,-1); lines[19] = vec4(0,0,1,0);
-  lines[20] = vec4(0,0,1,-1); lines[21] = vec4(1,0,0,0);
-  lines[22] = vec4(0,0,1,-1); lines[23] = vec4(0,1,0,0);
-
-  lines[24] = vec4(1,-1,0,0); lines[25] = vec4(0,0,1,-1);
-  lines[26] = vec4(1,0,-1,0); lines[27] = vec4(0,1,0,-1);
-  lines[28] = vec4(1,0,0,-1); lines[29] = vec4(0,1,-1,0);
-
-  lines[30] = vec4(1,PHI,-1,0); lines[31] = vec4(PHI,1,0,-1);
-  lines[32] = vec4(1,PHI,0,-1); lines[33] = vec4(PHI,1,-1,0);
-  lines[34] = vec4(1,-1,PHI,0); lines[35] = vec4(PHI,0,1,-1);
-  lines[36] = vec4(1,0,PHI,-1); lines[37] = vec4(PHI,-1,1,0);
-  lines[38] = vec4(1,-1,0,PHI); lines[39] = vec4(PHI,0,-1,1);
-  lines[40] = vec4(1,0,-1,PHI); lines[41] = vec4(PHI,-1,0,1);
-
-  lines[42] = vec4(-1,1,PHI,0); lines[43] = vec4(0,PHI,1,-1);
-  lines[44] = vec4(0,1,PHI,-1); lines[45] = vec4(-1,PHI,1,0);
-  lines[46] = vec4(-1,1,0,PHI); lines[47] = vec4(0,PHI,-1,1);
-  lines[48] = vec4(0,1,-1,PHI); lines[49] = vec4(-1,PHI,0,1);
-  lines[50] = vec4(-1,0,1,PHI); lines[51] = vec4(0,-1,PHI,1);
-  lines[52] = vec4(0,-1,1,PHI); lines[53] = vec4(-1,0,PHI,1);
-
-  colors[0] = colors[2] = colors[4] = colors[6] = colors[8] = colors[10] = 0;
-  colors[1] = colors[3] = colors[5] = colors[7] = colors[9] = colors[11] = 1;
-  colors[12] = colors[13] = colors[14] = 2;
-  colors[15] = colors[17] = colors[19] = colors[21] = colors[23] = colors[25] = 3;
-  colors[16] = colors[18] = colors[20] = colors[22] = colors[24] = colors[26] = 4;
+    x[0] = (A+B) - a/3.0;
+    return 1;
+  }
+}
+
+int cubic(float A, float B, float C, float D, out vec3 x) {
+  int nroots;
+  // Some ill-conditioned coeffs can cause problems
+  // The worst is fixed by solving for reciprocal
+  if (abs(A) > abs(D)) {
+    nroots = cubic0(A,B,C,D,x);
+  } else {
+    nroots = cubic0(D,C,B,A,x);
+    for (int i = 0; i < 3; i++) {
+      x[i] = 1.0/x[i];
+    }
+  }
+  return nroots;
+}
+
+// Initialize array of lines & colors
+void initlines() {
+  lines[0] =  m4inv*vec4(1,-1,0,0); lines[1] =  m4inv*vec4(0,0,1,0);
+  lines[2] =  m4inv*vec4(1,-1,0,0); lines[3] =  m4inv*vec4(0,0,0,1);
+  lines[4] =  m4inv*vec4(1,0,-1,0); lines[5] =  m4inv*vec4(0,1,0,0);
+  lines[6] =  m4inv*vec4(1,0,-1,0); lines[7] =  m4inv*vec4(0,0,0,1);
+  lines[8] =  m4inv*vec4(1,0,0,-1); lines[9] =  m4inv*vec4(0,1,0,0);
+  lines[10] = m4inv*vec4(1,0,0,-1); lines[11] = m4inv*vec4(0,0,1,0);
+
+  lines[12] = m4inv*vec4(0,1,-1,0); lines[13] = m4inv*vec4(1,0,0,0);
+  lines[14] = m4inv*vec4(0,1,-1,0); lines[15] = m4inv*vec4(0,0,0,1);
+  lines[16] = m4inv*vec4(0,1,0,-1); lines[17] = m4inv*vec4(1,0,0,0);
+  lines[18] = m4inv*vec4(0,1,0,-1); lines[19] = m4inv*vec4(0,0,1,0);
+  lines[20] = m4inv*vec4(0,0,1,-1); lines[21] = m4inv*vec4(1,0,0,0);
+  lines[22] = m4inv*vec4(0,0,1,-1); lines[23] = m4inv*vec4(0,1,0,0);
+
+  lines[24] = m4inv*vec4(1,-1,0,0); lines[25] = m4inv*vec4(0,0,1,-1);
+  lines[26] = m4inv*vec4(1,0,-1,0); lines[27] = m4inv*vec4(0,1,0,-1);
+  lines[28] = m4inv*vec4(1,0,0,-1); lines[29] = m4inv*vec4(0,1,-1,0);
+
+  lines[30] = m4inv*vec4(1,PHI,-1,0); lines[31] = m4inv*vec4(PHI,1,0,-1);
+  lines[32] = m4inv*vec4(1,PHI,0,-1); lines[33] = m4inv*vec4(PHI,1,-1,0);
+  lines[34] = m4inv*vec4(1,0,PHI,-1); lines[35] = m4inv*vec4(PHI,-1,1,0);
+  lines[36] = m4inv*vec4(1,-1,PHI,0); lines[37] = m4inv*vec4(PHI,0,1,-1);
+  lines[38] = m4inv*vec4(1,-1,0,PHI); lines[39] = m4inv*vec4(PHI,0,-1,1);
+  lines[40] = m4inv*vec4(1,0,-1,PHI); lines[41] = m4inv*vec4(PHI,-1,0,1);
+
+  lines[42] = m4inv*vec4(-1,1,PHI,0); lines[43] = m4inv*vec4(0,PHI,1,-1);
+  lines[44] = m4inv*vec4(0,1,PHI,-1); lines[45] = m4inv*vec4(-1,PHI,1,0);
+  lines[48] = m4inv*vec4(-1,1,0,PHI); lines[49] = m4inv*vec4(0,PHI,-1,1);
+  lines[46] = m4inv*vec4(0,1,-1,PHI); lines[47] = m4inv*vec4(-1,PHI,0,1);
+  lines[50] = m4inv*vec4(-1,0,1,PHI); lines[51] = m4inv*vec4(0,-1,PHI,1);
+  lines[52] = m4inv*vec4(0,-1,1,PHI); lines[53] = m4inv*vec4(-1,0,PHI,1);
+
+  colors[0] = colors[5] = colors[9] = 5;
+  colors[1] = colors[3] = colors[7] = 2;
+  colors[2] = colors[4] = colors[11] = 7;
+  colors[6] = colors[8] = colors[10] = 4;
+  colors[12] = colors[13] = colors[14] = 6;
+  // This colouring shows a "double six" for 12 lines.
+  colors[15] = colors[17] = colors[19] = colors[21] = colors[23] = colors[25] = 1;
+  colors[16] = colors[18] = colors[20] = colors[22] = colors[24] = colors[26] = 3;
 }
 
 vec4 getColor(int i) {
   if (i == 0) return vec4(0.5,0.5,1,1);
   if (i == 1) return vec4(1,0,0,1);
   if (i == 2) return vec4(1,1,0,1);
-  if (i == 3) return vec4(0,1,0,1);
-  if (i == 4) return vec4(0,0,1,1);
+  if (i == 3) return vec4(0,0,1,1);
+  if (i == 4) return vec4(0,1,0,1);
   if (i == 5) return vec4(0,1,1,1);
   if (i == 6) return vec4(1,0,1,1);
-  return vec4(0.5,0.5,0.5,1);
-}
-
-vec4 gradient(vec4 p) {
-  float x = p.x; float y = p.y; float z = p.z; float w = p.w;
-  float A = dot(p,I);
-  return -3.0*vec4((A+x)*(A-x),(A+y)*(A-y),(A+z)*(A-z),(A+w)*(A-w));
-}
-
-mat4 m;
-mat4 minv;
-
-// 2x^2 - (x+y+z+w)^3
-// x^3 + 3x^2(y+z+w) + 3x(yy+yz+yw+zz+zw+ww)
+  if (i == 7) return vec4(0.5,0.5,0.1,1);
+  return vec4(0,0,0,1);
+}
+
+// p is in world coordinates
+float eval(vec4 p) {
+  p = m4*p;
+  float pI = dot(p,I);
+  return dot(p*p,p) - pI*pI*pI;
+}
+
 // Equation is: pp.p - (p.I)^3 = 0
-bool surface(vec4 p, vec4 r, float min, inout float t, out vec3 normal) {
-  p = m*p;
-  r = m*r;
+// A = rr.r - (r.I)^3
+// B = 3(rr.p -(p.I)(r.I)^2)
+// C = 3(pp.r -(p.I)^2(r.I))
+// D = pp.p - (p.I)^3
+bool surface(vec4 p0, vec4 r0, float min, inout float t0, out vec3 normal) {
+  vec4 p = m4*p0;
+  vec4 r = m4*r0;
   float pI = dot(p,I);
   float rI = dot(r,I);
-#if 0
-  float A = dot(p*p,p)-pI*pI*pI;
-  float B = 3.0*(dot(p*p,r) - pI*pI*rI);
-  float C = 3.0*(dot(r*r,p) - pI*rI*rI);
-  float D = dot(r*r,r) - rI*rI*rI;
-#else
-  // Possibly this is a better way.
-  float A = dot(p*p - pI*pI*I, p);
-  float B = 3.0*dot(p*p - pI*pI*I, r);
-  float C = 3.0*dot(r*r - rI*rI*I, p);
-  float D = dot(r*r - rI*rI*I, r);
-#endif
-  if (!cubic(D,C,B,A,min,t)) return false;
-  vec4 tmp = p+t*r; // Actual surface point
-  vec4 n4 = minv*gradient(tmp);  // Rotate n4 here
-  normal = normalize(n4.xyz); // E3 component of P4 normal
-  return true;
+  float A = dot(r*r,r) - rI*rI*rI;
+  float B = 3.0*(dot(r*r,p) - pI*rI*rI);
+  float C = 3.0*(dot(p*p,r) - pI*pI*rI);
+  float D = dot(p*p,p)-pI*pI*pI;
+  vec3 x;
+  int nroots = cubic(A,B,C,D,x);
+  bool found = false;
+  for (int i = 0; i < 3; i++) {
+    if (i == nroots) break;
+    if (min < x[i] && x[i] < t0) {
+      t0 = x[i];
+      found = true;
+    }
+  }
+  // Compute normal by finding gradient
+  vec4 p1 = p0+t0*r0;
+  float x0 = eval(p1);
+  vec2 eps = vec2(1e-2,0);
+  normal = vec3 (eval(p1+eps.xyyy)-x0,
+                 eval(p1+eps.yxyy)-x0,
+                 eval(p1+eps.yyxy)-x0);
+  normal = normalize(normal);
+  return found;
 }
     
-// A = pp.p - (p.I)^3
-// B = 3(pp.r -(p.I)^2(r.I))
-// C = 3(rr.p -(p.I)(r.I)^2)
-// D = rr.r - (r.I)^3
-
-vec3 light;
-float ambient;
-float diffuse;
+// Find the distance to the line in R3
+void tryline(vec4 q, vec4 p, vec4 r, int color) {
+  const float u = 0.025;
+  const float u2 = u*u;
+  if (abs(p.w) < 1e-2) return;
+   r = p.w*r - r.w*p; // r.w = 0
+   q = p.w*q-p;
+   float qr = dot(q,r);
+   float d2 = (dot(q,q) - (qr*qr/dot(r,r)))/(p.w*p.w);
+   if (d2 < min(u2,d2best)) {
+      d2best = d2;
+      colorindex = color;
+   }
+}
 
 vec4 solve(vec3 p, vec3 r, float min) {
-  float t = 1e10;
-  float u = 0.015;
-  float u2 = u*u;
-  vec3 q0, s0;
-  int colorindex = -1;
   vec3 normal = vec3(1,0,0);
   float theta = iGlobalTime*0.1; 
-  mat4 minv0 = 0.25*mat4(-1,1,1,1, 1,-1,1,1, 1,1,-1,1, -1,-1,-1,1);
-  mat4 m0 = 0.25*mat4(-1,1,1,-1, 1,-1,1,-1, 1,1,-1,-1, 1,1,1,1);
-#if 0
-  m = mat4(1.0, 0.0, 0.0, 0.0,
-                0.0, cos(theta), 0.0, -sin(theta),
-                0.0, 0.0, 1.0, 0.0,
-                0.0, sin(theta), 0.0, cos(theta));
-  minv = mat4(1.0, 0.0, 0.0, 0.0,
-                0.0, cos(theta), 0.0, sin(theta),
-                0.0, 0.0, 1.0, 0.0,
-                0.0, -sin(theta), 0.0, cos(theta));
-#elif 0
-  m = mat4(cos(theta), 0.0, -sin(theta), 0.0,
-           0.0, 1.0, 0.0, 0.0,
-           sin(theta), 0.0, cos(theta), 0.0,
-                0.0, 0.0, 0.0, 1.0);
-  minv = mat4(cos(theta), 0.0, sin(theta), 0.0,
-                0.0, 1.0, 0.0, 0.0,
-                -sin(theta), 0.0, cos(theta), 0.0,
-                0.0, 0.0, 0.0, 1.0);
-#else
-  m = qmat(vec4(0.0,0.0,sin(theta),cos(theta)));
-  minv = qmat(vec4(0.0,0.0,sin(theta),-cos(theta)));
-#endif
-  //m = m0*m; minv = minv*minv0;
-
-  bool dosurface = true;
-  bool dolines = true;
-  if (dosurface && surface(vec4(p,1),vec4(r,0),min,t,normal)) {
+  vec3 a = normalize(AXIS);
+  m4 = qmat(vec4(sin(theta)*a,cos(theta)));
+  m4inv = qmat(vec4(sin(-theta)*a,cos(-theta)));
+  initlines();
+
+  float t = 1e10;
+  if (surface(vec4(p,1),vec4(r,0),min,t,normal)) {
     colorindex = 0;
-  }
-  if (dolines) {
+    vec4 q = vec4(p+t*r,1);
     for (int i = 0; i < NLINES; i++) {
-      vec3 q,s;
-      {
-        // Convert 2 projective points into an E3
-        // point+direction line. I'm sure there is
-        // a cleaner way to do this, but this seems
-        // to work.
-        vec4 p0 = minv*lines[2*i];
-        vec4 p1 = minv*lines[2*i+1];
-        float eps = 1e-1;
-        if (abs(p0.w) < eps) p0 += p1;
-        if (abs(p1.w) < eps) p1 += p0;
-        if (abs(p0.w) < eps*eps) continue;
-        if (abs(p1.w) < eps*eps) continue;
-        q = project(p0);
-        s = project(p1);
-        s = normalize(s-q); // Line direction
-      }
-      //q = vec3(0,0,0);
-      //s = vec3(0,1,0);
-      vec3 p1 = p-q; // Move p to line space
-      float rs = dot(r,s);
-      float ps = dot(p1,s);
-      float pr = dot(p1,r);
-      float pp = dot(p1,p1);
-      float A = 1.0 - rs*rs;
-      float B = 2.0*(pr - ps*rs);
-      float C = pp - ps*ps - u2;
-      if (quadratic(A,B,C,min,t)) {
-        q0 = q; s0 = s; colorindex = 1+colors[i];
-        p1 += t*r; // Final point in line space
-        p1 -= dot(p1,s)*s; // Subtract parallel part
-        normal = normalize(p1);
+      vec4 p = lines[2*i];
+      vec4 r = lines[2*i+1];
+      if (abs(p.w) < abs(r.w)) {
+        tryline(q,r,p,colors[i]);
+      } else {
+        tryline(q,p,r,colors[i]);
       }
     }
   }
-  if (highlight) colorindex = 1;
   if (colorindex == -1) return vec4(0,0,0,1);
+  if (dot(r,normal) > 0.0) {
+    normal = -normal; // Face forwards
+  }
   vec4 baseColor = getColor(colorindex);
   vec3 color;
   color = baseColor.xyz*(ambient+(1.0-ambient)*dot(light,normal));

  ambient = 0.6;
  diffuse = 1.0-ambient;

  init();
  vec2 uv = 2.0*fragCoord.xy/iResolution.xy - 1.0;
  vec3 p = vec3(0.0, 0.0, 10.0);
  vec3 r = normalize(vec3(iResolution.x/iResolution.y * uv.x, uv.y, -3.0));
  float k = -dot(p,r);
  //k = 0.0;
  p += k*r;
  fragColor = solve(p,r,-k);
}

const float PI =  3.141592654;
const float PHI = 1.618033989;
const vec4 I = vec4(1,1,1,1);
const vec3 AXIS = vec3(1,1,1); // axis for rotation
const int NLINES = 27;
vec4 lines[NLINES*2];
int colors[NLINES];

// Rotation & inverse in R4
mat4 m4;
mat4 m4inv;

// Lighting params
vec3 light;
float ambient;
float diffuse;

float d2best = 1e10; // Smallest distance to line
int colorindex = -1; // Selected color

// Quaternion multiplication as a matrix.
// w coordinate is real element of quaternion
