Processing3; for exploring conformal branched self-covers of the Poincaré disc. it's a bit rough, at the moment, but usable.

cxarith.pde

/************ basic complex functions ***
* I once tried to implement "double" versions of these, too; something hung-up...
*/

float[] cOne = {1.0,0.0};
float[] cZro = {0.0,0.0};
float[] cmOne = { -1.0, 0.0 };

float[] cmul ( float[] c1 , float[] c2 ) {
    float ans[] = new float[2];
    ans[0] = c1[0] * c2[0] - c1[1] * c2[1];
    ans[1] = c1[0] * c2[1] + c1[1] * c2[0];
    return ans; 
}

float[] rmul ( float r , float[] c1){
  float ans[] = new float[2];
  ans[0] = c1[0] * r ;
  ans[1] = c1[1] * r ;
  return ans;
}

float[] cadd ( float[] c1, float[] c2 ) {
  float ans[] = new float[2];
  ans[0] = c1[0] + c2[0];
  ans[1] = c1[1] + c2[1];
  return ans;
}

float[] cdif( float[] c1, float[] c2) {
  float ans[] = new float[2];
  ans[0] = c1[0] - c2[0];
  ans[1] = c1[1] - c2[1];
  return ans;
}

float[] conj ( float[] c1 ) {
  float[] ans = new float[2];
  ans[0] = c1[0];
  ans[1] = -c1[1];
  return ans;
}

float cnorm( float[] c1 ) {
  return (c1[0] * c1[0] + c1[1] * c1[1]);
}

float[] cinv ( float[] c1){
  return rmul ( 1/cnorm(c1), conj(c1) );
}

float[] cdiv ( float[] c1, float[] c2) {
  float ans[] = new float[2];
  float R = cnorm(c2);
  ans[0] = (c1[0]*c2[0] + c1[1]*c2[1])/R ;
  ans[1] = (-c1[0]*c2[1]+c1[1]*c2[0])/R;
  return ans;
}

float[] cexp ( float[] c1 ){
  float ans[] = new float[2];
  float R = exp(c1[0]);
  ans[0] = R * cos(c1[1]);
  ans[1] = R * sin(c1[1]);
  return ans;
}

float[] clog( float[] c1){
  float ans[] = new float[2];
  float R = sqrt(cnorm(c1));
  ans[0] = log(R);
  ans[1] = 2 * atan(c1[1]/(R + c1[0]));
  return ans;
}

float[] ipower(int n, float[] z){
  int j = n;
  if ( j < 0 ){
    j = -j;
    z = cinv(z);
  }
  float scale[] = { z[0] , z[1] };
  float ans[] = { 1 , 0 };
  while(j > 0){
    if ( (j & 1) != 0 ){
      ans = cmul (ans , scale);
    }
    j = j/2;
    scale = cmul(scale,scale);
  }
  return ans;
}

float[] roughPower(float[] a , float[] z){
  return cexp( cmul ( a , clog(z)));
}

float[] roughPower(float a , float[] z){
  return cexp( rmul ( a , clog(z)));
}

discFolding.pde

import interfascia.*; // we don't do anything with that, yet.

hdivisor ds[] = new hdivisor[124]; //what's an hdivisor? I'm So Glad You Asked!
int topOfTheStack = -1;

float[] thePolyEval(float[] zed) {
  float[] ans = { zed[0] , zed[1] };
  int j;
  for (j = topOfTheStack ; j >= 0 ; j--){
    ans = ds[j].eval(ans);
  }
  return ans;
}

color lemnColor(float[] zed){ // because we need to choose colors *Somehow*, and programmatically is better than Pixelated.
  float twoz[] = rmul(2,zed);
  float p = cnorm (cdif (twoz,cOne));
  float q = cnorm (cdif (twoz,cmOne));
  if ( p * q < 1 ) {
    if (p < q) return color(240,0,0);
    else return color(0,0,240);
  } 
  else 
  return color(255,255,255);
}

void setup(){
  size(400,400);
  loadPixels();
  noFill();
}

void draw(){
  float[] xy = new float[2];
  float[] pz;
  int rc;
  int cx, cy;
  int thePixel = 0 ;
  for ( int row = 0 ; row < 400 ; row++ ) {
    xy[1] = ( 200.0 - row ) / 200; 
    for ( int col = 0 ; col < 400 ; col++ ) {
      xy[0] = ( col - 200.0 ) / 200;
      pz = thePolyEval(xy);
      if ( cnorm(pz) > 3 ) pixels[thePixel] = color(255,255,255);
      else {
        cx = floor(400.0 + xy[1] * 200);
        cy = floor(400.0 - xy[0] * 200);
        pixels[thePixel] = lemnColor(pz);
      } 
      thePixel++;
    }
  }
  updatePixels();
  ellipse(200,200,400,400);
}

void mouseClicked(){
  int order = 1;
  float[] rr;
  float[] tgt = {(mouseX - 200.0)/200 , (200.0 - mouseY)/200 };
  if (topOfTheStack < 123)
  {
    switch ( mouseButton ) {
    case CENTER:
      order++;
    case RIGHT:
      order++;
    case LEFT:
      order++;
    }  
    topOfTheStack++;
    rr = rmul(-1, roughPower( 1.0 /order , rmul(-1, tgt)));
    ds[topOfTheStack] = new hdivisor( order, rr[0], rr[1] );
  }
}

hdivisor.pde

class hdivisor{
  // for branched covers of The Poincaré disc.
  // our disc will have radius 2, just 'cauz.
  private
  float pair[];
  float untwist[];
  public
  float re;
  float im;
  int degree;
  hdivisor( int dg , float rr , float ii){
    pair = new float[2];
    re = rr;
    im = ii;
    degree =dg;
    pair[0] = re;
    pair[1] = im;
    untwist = ipower(degree,cdiv( cdif(cZro,pair), cdif ( cOne , cmul(cZro,conj(pair)))));
  }
  public float[] eval( float[] z ) { 
    float[] mov1 = cdiv( cdif(z,pair), cdif ( cOne ,cmul(z,conj(pair))));
    float[] mov2 = ipower(degree,mov1);
    return       cdiv( cdif(mov2,untwist),cdif(cOne, cmul(mov2,conj(untwist))));
  }
}