#define PI 3.141592
#define TWOPI 6.283185
#define S2inv 0.707106
#define R vec2(1.,0.)
#define I vec2(0.,1.)
#define iters 6

vec2 cmul(vec2 z, vec2 c) {
    return vec2(z.x * c.x - z.y * c.y, z.x * c.y + z.y * c.x);
}

vec2 cdiv(vec2 z, vec2 c) {
    float r = dot(c, c);
    return vec2(z.x * c.x + z.y * c.y, z.y * c.x - z.x * c.y) / r;
}


vec2 mobius(vec2 z, vec2 a, vec2 b, vec2 c, vec2 d) {
    return cdiv(cmul(a,z)+b,cmul(c,z)+d);
}

vec2 cpow(vec2 z, vec2 p) {
    float a = atan(z.y, z.x);
    float lnr = 0.5 * log(dot(z,z));
    float m = exp(p.x * lnr - p.y * a);
    float angle = p.x * a + p.y * lnr + TWOPI;
    return vec2(cos(angle), sin(angle)) * m;
}

vec2 clog(vec2 z) {
    return vec2(0.5 * log(dot(z,z)), atan(z.y,z.x));
}

vec2 catan(vec2 z) {
    return clog(mobius(z,R,R,-R,R));
}

vec2 casin(vec2 z) {
    return clog(z+cpow(cmul(z,z)+R,R*0.5));
}

vec2 cagm(vec2 ari, vec2 geo) {
    vec2 a = ari;
    vec2 g = geo;
    for(int i = 0; i < 10; i++) {
        a = cdiv(a + g,2.*R);
        g = cpow(cmul(a,g),R*0.5);
    }
    return g;
}

vec2 cahm(vec2 ari, vec2 har) {
    vec2 a = ari;
    vec2 h = har;
    for(int i = 0; i < 10; i++) {
        a = cdiv(a + h,2.*R);
        h = cdiv(2.*R,cdiv(R,a)+cdiv(R,h));
    }
    return h;
}

vec2 cghm(vec2 geo, vec2 har) {
    vec2 g = geo;
    vec2 h = har;
    for(int i = 0; i < 10; i++) {
        g = cpow(cmul(g,h),R*0.5);
        h = cdiv(2.*R,cdiv(R,g)+cdiv(R,h));
    }
    return g;
}

float l(float r) {
    return 2. / PI * atan(r);
}

float hue2rgb(float p, float q, float t) {
    do{
      if(t < 0.0) t += 1.0;
      if(t > 1.0) t -= 1.0;
    } while (t < 0.0 || t > 1.0);

  if(t < 1.0 / 6.0) return p + (q - p) * 6.0 * t;
  if(t < 1.0 / 2.0) return q;
  if(t < 2.0 / 3.0) return p + (q - p) * (2.0 / 3.0 - t) * 6.0;
  return p;
}

vec3 hslToRgb(float h, float s, float l) {
  float r, g, b;

  if(s == 0.0) {
    r = g = b = l; // achromatic
  } else {
    float q = l < 0.5 ? l * (1.0 + s) : l + s - l * s;
    float p = 2.0 * l - q;

    r = hue2rgb(p, q, h + 1.0 / 3.0);
    g = hue2rgb(p, q, h);
    b = hue2rgb(p, q, h - 1.0 / 3.0);
  }

  return vec3(r,g,b);
}

vec3 domainColoring(vec2 z) {
    float H = mod(iTime*0.0125,1.)-atan(z.y,z.x)/ TWOPI - TWOPI / 3.0;
    float S = 1.0;
    float L = l(length(z));
    return hslToRgb(H,S,L);
}

vec4 fC(vec2 fragCoord )
{
    vec2 uv = 2.0 * (fragCoord.xy - 0.5*iResolution.xy) / -iResolution.y;

    float t = iTime * 0.125;

    uv.x += 0.5;
    uv = dot(uv,uv)<1.0?uv:vec2(-uv.x,uv.y);
    uv = mobius(uv,R,-R,R,R);
    uv.x = -abs(uv.x);
    uv = mobius(uv,R,R,-R,R);

    uv = mobius(uv,R,0.75*R,0.75*R,R);

    uv = cmul(uv,vec2(cos(t),sin(t)));

    uv *= 14.01;

    uv = cdiv(R,uv);

    for(int i = 0; i < 20; i++){
        uv = cahm(cmul(uv,uv),cdiv(R,cmul(uv,uv)));
    }

    vec2 nuv = uv/(0.1*PI);

    float grid = mod(floor(nuv.x)+floor(nuv.y),2.);

    vec3 col = -0.5 * length(vec2(grid)) + domainColoring(uv);

    // Output to screen
    return vec4(col,1.0);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    fragColor = vec4(0);
    float A = 2.,
          s = 1./A, x, y;

    for (x=-.5; x<.5; x+=s) for (y=-.5; y<.5; y+=s) fragColor += min ( fC(vec2(x,y)+fragCoord), 1.0);

	fragColor /= A*A;
}