#version 130

#include "Progressive2D.frag"

#group TriangularMandelbrot

uniform int Iterations; slider[1,10,1000]
uniform float EscapeRadius; slider[0,2,100]
// where to put p.xyz into the 6D C in "Z \to Z^2+C"
uniform int SliceA; slider[0,0,5]
uniform int SliceB; slider[0,3,5]
// no float[6] sliders...
uniform float CenterU; slider[-64.0,0.0,64.0]
uniform float CenterV; slider[-64.0,0.0,64.0]
uniform float CenterW; slider[-64.0,0.0,64.0]
uniform float CenterX; slider[-64.0,0.0,64.0]
uniform float CenterY; slider[-64.0,0.0,64.0]
uniform float CenterZ; slider[-64.0,0.0,64.0]

// commutative addition
float[6] add(float[6] a, float[6] b) {
  return float[6]
    ( a[0] + b[0]
    , a[1] + b[1]
    , a[2] + b[2]
    , a[3] + b[3]
    , a[4] + b[4]
    , a[5] + b[5]
    );
}

// non-commutative multiplication
float[6] mul(float[6] a, float[6] b) {
  return float[6]
    ( a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[4] + a[4] * b[3] + a[5] * b[5]
    , a[0] * b[1] + a[1] * b[0] + a[2] * b[3] + a[3] * b[5] + a[4] * b[2] + a[5] * b[4]
    , a[0] * b[2] + a[1] * b[4] + a[2] * b[0] + a[3] * b[1] + a[4] * b[5] + a[5] * b[3]
    , a[0] * b[3] + a[1] * b[5] + a[2] * b[1] + a[3] * b[0] + a[4] * b[4] + a[5] * b[2]
    , a[0] * b[4] + a[1] * b[2] + a[2] * b[5] + a[3] * b[3] + a[4] * b[0] + a[5] * b[1]
    , a[0] * b[5] + a[1] * b[3] + a[2] * b[4] + a[3] * b[2] + a[4] * b[1] + a[5] * b[0]
    );
}

// attempt at exterior distance estimate
vec3 color(vec2 p) {
  // prepare C
  float c[6] = float[6]
    ( CenterU
    , CenterV
    , CenterW
    , CenterX
    , CenterY
    , CenterZ
    );
  c[SliceA] += p.x;
  c[SliceB] += p.y;
  // prepare Z and DZ
  float zero[6] = float[6](0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
  float one[6] = float[6](1.0, 0.0, 0.0, 0.0, 0.0, 0.0);
  float z[6] = zero; // iterate
  float dc[6] = zero; // derivative
  float normz = 0.0, normdc = 0.0;
  for (int i = 0; i < Iterations; ++i) {
    normz = 0.0;
    normdc = 0.0;
    for (int j = 0; j < 6; ++j) {
      normz += z[j] * z[j];
      normdc += dc[j] * dc[j];
    }
    if (normz > EscapeRadius * EscapeRadius) {
      // escaped
      break;
    }
    // apply product rule, mul is non-commutative
    dc = add(add(mul(z, dc), mul(dc, z)), one);
    z = add(mul(z, z), c);
  }
  float r = sqrt(normz);
  float dr = sqrt(normdc);
  float de = 0.5 * r * log(r) / dr;
  float px = length(dFdy(p));
  return r < EscapeRadius ? vec3(1.0/(r + 1.0),0.0,0.0) : vec3(tanh(de / px));
}