Mitch Richling: Apollonian Gasket

Author:	Mitch Richling
Updated:	2022-10-05

1. Apollonian Gaskets
2. Generating Circles and Spheres
3. Rendering
4. An IFS Approach
5. References

1. Apollonian Gaskets

An Apollonian gasket fractal is generated by starting with a triple of mutually tangent circles, and successively filling in more circles, each tangent to another three. The result is sometimes called the curvilinear Sierpinski gasket. You can read all about it on wikipedia

A single Apollonian gasket is fun, but four is better. ;) For many sets of three mutually tangent circles one may find an enclosing circle that is also mutually tangent. Now one may create four different Apollonian gaskets, one for each possible 3-tuple of circles. Here are a couple examples:

If we turn each of the circles in a gasket into spheres, then we get the construction at the top of the page. Here are some examples:

2. Generating Circles and Spheres

The first step in creating images like the ones above is to generate descriptions of the geometric objects (circle/sphere centers and radii). Practically speaking this means creating a files in some "standard" format with circle data for 2D renderings and spheres data for 3D renderings. The following little ruby script will generate SVG files for 2D data:

#!/bin/bash /home/richmit/bin/ruby
# -*- Mode:Ruby; Coding:us-ascii-unix; fill-column:132 -*-
####################################################################################################################################
##
# @file      gen.rb
# @author    Mitch Richling https://www.mitchr.me
# @date      2017-05-16
# @brief     Generate geometry files (SVG or Pov-Ray) representing an Apollonian Gasket.@EOL
# @keywords  fractal geometry SVG Pov-Ray Apollonian Gasket
# @std       Ruby 2.0
####################################################################################################################################

#-----------------------------------------------------------------------------------------------------------------------------------
require 'cmath'
require 'optparse'
require 'optparse/time'

#-----------------------------------------------------------------------------------------------------------------------------------
# Command line parse and generation parameters

outputFileName = 'apollonianFractal'
$doSVG         = FALSE
$doPOV         = TRUE
r1             = 1.0
r2             = 1.0
r3             = 1.0
$curvLimit     = 20000
$povClipRad    = 0.5
$povMinFakeRad = 0.003

debug = 0
opts = OptionParser.new do |opts|
  opts.banner = "Usage: agen.rb [options]"
  opts.separator ""
  opts.separator "Options:"
  opts.on("-h",        "--help",             "Show this message")                      { puts opts; exit                 }
  opts.on(             "--r1 RADIUS",        "Radius for first circle")                { |v| r1=v.to_f                   }
  opts.on(             "--r2 RADIUS",        "Radius for second circle")               { |v| r2=v.to_f                   }
  opts.on(             "--r3 RADIUS",        "Radius for third circle")                { |v| r3=v.to_f                   }
  opts.on(             "--doSVG T/F",        "Dump SVG file")                          { |v| $doSVG=v.match(/^[YyTt1]/)  }
  opts.on(             "--doPOV T/F",        "Dump Povray include file")               { |v| $doPOV=v.match(/^[YyTt1]/)  }
  opts.on(             "--output BASE",      "Name for output file(s)")                { |v| outputFileName=v            }
  opts.on(             "--crvLim VALUE",     "Drop circles with large curvature")      { |v| $curvLimit=v.to_f           }
  opts.on(             "--povClipRad VALUE", "Drop spheres bigger than this")          { |v| $povClipRad=v.to_f          }
  opts.on(             "--povMinRad VALUE",  "Puff up small spheres to this size")     { |v| $povMinFakeRad=v.to_f       }
  opts.separator "                                                                             "
  opts.separator "Generate geometry files (SVG or Pov-Ray) representing an Apollonian Gasket.  "
  opts.separator "Inputs are the starting radii of the initial interior circles.  The exterior "
  opts.separator "circle is automatically generated.  Note that not all radii will work, and   "
  opts.separator "we don't do any error checking.  Circle generation stops when radii become   "
  opts.separator "too small.  The Pov-Ray output is in the form of union object with a version "
  opts.separator "declaration of 3.7.  The SVG output is 1 pix wide, so use a converter that   "
  opts.separator "will do the right thing or adjust the width.                                 "
  opts.separator "                                                                             "
  opts.separator "Examples:                                                                    "
  opts.separator "   * Classic apollonian fractal starting with three circles of equal size.   "
  opts.separator "     * agen.rb --doSVG F --doPOV T --output apf2 --r1 1.0    --r1 1.0    --r1 1.0     --crvLim 20000.0 --povClipRad 0.5 --povMinRad 0.003"
  opts.separator "     * agen.rb --doSVG T --doPOV F --output apf2 --r1 1000.0 --r1 1000.0 --r1 1000.0  --crvLim 1.0"
  opts.separator "   * The tow eye apollonian fractal starting with two bit circles of equal size centered in the enclosing circle."
  opts.separator "     * agen.rb --doSVG T --doPOV F --output apf3 --r1 1000.0 --r1 1000.0 --r1 666.6666666666  --crvLim 1.0"
end
opts.parse!(ARGV)

if( !($doSVG || $doPOV)) then
  puts("ERROR: No output will be produced (see --doPOV & --doSVG)")
  exit
end

#-----------------------------------------------------------------------------------------------------------------------------------

$circles = Hash.new   # The primary data structure with all the geometry
$quads   = Array.new  # Used for holding lists of tangent circles

#-----------------------------------------------------------------------------------------------------------------------------------
# Given three circles, produce the two tangent ones -- if possible.  No error checking, so be careful.  We only use this for the
# enclosing circle.  Method is that of Descartes.
#
# https://en.wikipedia.org/wiki/Descartes'_theorem

def newCircles (circle1, circle2, circle3)
  k1, k2, k3 = circle1[0], circle2[0], circle3[0]
  c1, c2, c3 = circle1[1], circle2[1], circle3[1]
  # Compute the inner and outer circle radius
  tmp1 = k1+k2+k3
  tmp2 = 2*Math::sqrt(k1*k2+k2*k3+k3*k1)
  iRad = tmp1+tmp2
  oRad = tmp1-tmp2
  if(iRad < oRad) then
    iRad, oRad = oRad, iRad
  end
  # Compute the coordinates
  tmp1 = c1*k1+c2*k2+c3*k3
  tmp2 = 2*CMath::sqrt(k1*k2*c1*c2+k2*k3*c2*c3+k3*k1*c1*c3)

  circles = [ [iRad, (tmp1+tmp2)/iRad ],
              [oRad, (tmp1-tmp2)/oRad ]
            ]
  return circles  
end

#-----------------------------------------------------------------------------------------------------------------------------------
# Given four circles, find an alternate solution for the first one.  This is a matter of finding a second root to a polynomial once
# you already know one.

def otherSol (circle1, circle2, circle3, circle4)
  k1, k2, k3, k4 = circle1[0], circle2[0], circle3[0], circle4[0]
  c1, c2, c3, c4 = circle1[1], circle2[1], circle3[1], circle4[1]
  newK1 = 2*(k2+k3+k4) - k1
  newC1 = (2*(k2*c2+k3*c3+k4*c4) - k1*c1) / newK1
  return [ newK1, newC1 ]
end

#-----------------------------------------------------------------------------------------------------------------------------------
# Add a circle to our master list.  We return a hash for the circle if the add was successful.  It won't be successful if the circle
# was already on the list or if the circle radius was too small.  Note the custom has key -- the idea is to call circles the same if
# radius and coordinates match to 6 digits.

def addAcircle(circle)
  theKey = sprintf("%0.5f/%0.5f/%0.5f", circle[0], circle[1].real, circle[1].imag);
  if( ($circles.member?(theKey)) || (circle[0].abs>$curvLimit) ) then
    return nil
  else
    $circles[theKey] = circle
    return theKey
  end
end

#-----------------------------------------------------------------------------------------------------------------------------------
# Create first quad element

tmp = (r1*r1+r1*r3+r1*r2-r2*r3)/(r1+r2)

$quads.push(Array.new)

$quads[0].push(addAcircle([ 1/r1, Complex(    0,                             0) ]))
$quads[0].push(addAcircle([ 1/r2, Complex(r1+r2,                             0) ]))
$quads[0].push(addAcircle([ 1/r3, Complex(  tmp, Math::sqrt((r1+r3)**2-tmp**2)) ]))

outsideCircle = addAcircle(newCircles(*($circles.values_at(*($quads[0]))))[1])

$quads[0].unshift(outsideCircle)

#-----------------------------------------------------------------------------------------------------------------------------------
# Compute circles

while(curQuad = $quads.shift) do
  a, b, c, d = curQuad
  tmp=addAcircle(otherSol($circles[a], $circles[b], $circles[c], $circles[d])); if(!tmp.nil?) then $quads.push([tmp, b, c, d]); end
  tmp=addAcircle(otherSol($circles[d], $circles[a], $circles[b], $circles[c])); if(!tmp.nil?) then $quads.push([tmp, a, b, c]); end
  tmp=addAcircle(otherSol($circles[b], $circles[a], $circles[c], $circles[d])); if(!tmp.nil?) then $quads.push([tmp, a, c, d]); end
  tmp=addAcircle(otherSol($circles[c], $circles[a], $circles[b], $circles[d])); if(!tmp.nil?) then $quads.push([tmp, a, b, d]); end
end

#-----------------------------------------------------------------------------------------------------------------------------------
# Dump SVG

if($doSVG) then
  open(outputFileName + '.svg', 'w') do |file|
    file.puts("<svg width='#{1}px' height='#{1}px'>")
    $circles.each do |key, circle|
      k, c = circle
      file.puts("    <circle cx='#{500+c.real}' cy='#{500+c.imag}' r='#{(1/k).abs}' fill='none' stroke='red'/>")
    end
    file.puts('</svg>')
  end
end

#-----------------------------------------------------------------------------------------------------------------------------------
# Dump Povray include file

if($doPOV) then
  open(outputFileName + '.inc', 'w') do |file|
    file.puts('#version 3.7;')
    file.puts('#declare apollonianfractal = union {')
    ox,oy = $circles[outsideCircle][1].rect
    $circles.each do |key, circle|
      k, c = circle
      r = (1/k).abs
      if(r < $povMinFakeRad) then
        r = $povMinFakeRad
      end
      if(r<=$povClipRad) then
        file.puts("  sphere {<#{c.real-ox},#{c.imag-oy},#{0}> #{r}}")
      end
    end
    file.puts('}')
  end
end

3. Rendering

One need not explicitly render SVG files to view them because many software applications do that behind the scenes; however, one can obtain nice raster images from SVG files using ImageMagick which uses inkscape internally, to convert them into PNG files. For the 3D scenes, I used POV-Ray with the following thre input files:

// -*- Mode:pov; Coding:us-ascii-unix; fill-column:132 -*-
/***********************************************************************************************************************************/
/**
 @file      ap.pov
 @author    Mitch Richling https://www.mitchr.me
 @brief     Apollonian Gasket.@EOL
 @std       Povray_3.7
************************************************************************************************************************************/

//----------------------------------------------------------------------------------------------------------------------------------
#version 3.7;

//----------------------------------------------------------------------------------------------------------------------------------
#include "colors.inc"
#include "textures.inc"
#include "metals.inc"
#include "finish.inc"

//----------------------------------------------------------------------------------------------------------------------------------
global_settings { 
  assumed_gamma   1
  max_trace_level 2
}

//----------------------------------------------------------------------------------------------------------------------------------
camera {
  location <5, 7, 4.5>    
  up       <0,0,1>
  sky      <0,0,1>
  right    <0,4/2,0>
  look_at  <0,0,-.2>         
  angle    28
}

//----------------------------------------------------------------------------------------------------------------------------------
light_source { < 2.0, 5.0,  4.0> color White*0.55 }
light_source { < 4.0, 7.0, -2.0> color White*0.55 }
light_source { < 3.0, 2.0,  4.5> color White*0.55 }
light_source { < 5.0, 8.0,  5.4> color White*0.55 }
light_source { < 5.0, 0.0,  0.0> color White*0.55 }

//----------------------------------------------------------------------------------------------------------------------------------
background { color Black }

//----------------------------------------------------------------------------------------------------------------------------------
object{
  apollonianfractal
  material {
    texture {
      pigment { rgbf <0.1, 0.1, 1.0, 0.0> quick_color Red }
      finish  {
        ambient 0.1
        diffuse 0.4
        specular 0.2
        reflection { 0.0 0.5 }
        roughness 0.01
        phong 0.5
        phong_size 8
      }
    }
  }
}

// -*- Mode:pov; Coding:us-ascii-unix; fill-column:132 -*-
/***********************************************************************************************************************************/
/**
 @file      ap3.pov
 @author    Mitch Richling https://www.mitchr.me
 @brief     Apollonian Gasket.@EOL
 @std       Povray_3.7
************************************************************************************************************************************/

//----------------------------------------------------------------------------------------------------------------------------------
#version 3.7;

//----------------------------------------------------------------------------------------------------------------------------------
#include "colors.inc"
#include "textures.inc"
#include "metals.inc"
#include "finish.inc"

//----------------------------------------------------------------------------------------------------------------------------------
global_settings { 
  assumed_gamma   1
  max_trace_level 2
}

//----------------------------------------------------------------------------------------------------------------------------------
camera {
  location <5, 7, 4.5>    
  up       <0,0,1>
  sky      <0,0,1>
  right    <0,4/2,0>
  look_at  <0,0,-.1>         
  angle    22
}

//----------------------------------------------------------------------------------------------------------------------------------
light_source { < 2.0, 5.0,  4.0> color White*0.55 }
light_source { < 4.0, 7.0, -2.0> color White*0.55 }
light_source { < 3.0, 2.0,  4.5> color White*0.55 }
light_source { < 5.0, 8.0,  5.4> color White*0.55 }
light_source { < 5.0, 0.0,  0.0> color White*0.55 }

//----------------------------------------------------------------------------------------------------------------------------------
background { color Black }

//----------------------------------------------------------------------------------------------------------------------------------
object{
  apollonianfractal
  rotate 200*z
  material {
    texture {
      pigment { rgbf <0.1, 0.1, 1.0, 0.0> quick_color Red }
      finish  {
        ambient 0.1
        diffuse 0.4
        specular 0.2
        reflection { 0.0 0.5 }
        roughness 0.01
        phong 0.5
        phong_size 8
      }
    }
  }
}

// -*- Mode:pov; Coding:us-ascii-unix; fill-column:132 -*-
/***********************************************************************************************************************************/
/**
 @file      ap.pov
 @author    Mitch Richling https://www.mitchr.me
 @brief     Apollonian Gasket.@EOL
 @std       Povray_3.7
************************************************************************************************************************************/

//----------------------------------------------------------------------------------------------------------------------------------
#version 3.7;

//----------------------------------------------------------------------------------------------------------------------------------
#include "colors.inc"
#include "textures.inc"
#include "metals.inc"
#include "finish.inc"

//----------------------------------------------------------------------------------------------------------------------------------
global_settings { 
  assumed_gamma   1
  max_trace_level 2
}

//----------------------------------------------------------------------------------------------------------------------------------
camera {
  location <5, 7, 4.5>    
  up       <0,0,1>
  sky      <0,0,1>
  right    <0,4/2,0>
  look_at  <0,0,-.2>         
  angle    28
}

//----------------------------------------------------------------------------------------------------------------------------------
light_source { < 3.0, 5.0,  2.0> color White*0.55 }
light_source { < 6.0, 5.0,  1.5> color White*0.55 }
light_source { < 5.0, 8.0,  3.4> color White*0.55 }

//----------------------------------------------------------------------------------------------------------------------------------
background { color White }

//----------------------------------------------------------------------------------------------------------------------------------
object{
  apollonianfractal
  rotate 90*clock*z
  material {
    texture {
      pigment { rgbf <0.1, 0.1, 1.0, 0.0> quick_color Red }
      finish  {
        ambient 0.1
        diffuse 0.4
        specular 0.2
        reflection { 0.0 0.5 }
        roughness 0.01
        phong 0.5
        phong_size 8
      }
    }
  }
}

Note these are the primary scene description files, and the sphere data is housed in "include" files generated by the script. For examples of how to use the rendering tools mentioned above, see the makefile:

# -*- Mode:Makefile; Coding:utf-8; fill-column:132 -*-
####################################################################################################################################
##
# @file      makefile
# @author    Mitch Richling https://www.mitchr.me
# @brief     Apollonian Gasket Related Images. @EOL
# @std       GNUmake
####################################################################################################################################

#-----------------------------------------------------------------------------------------------------------------------------------
WRES := 4
HRES := 2
RESM := 000

#-----------------------------------------------------------------------------------------------------------------------------------
all : www

www : apollonian_gasket_01w_1000.jpg apollonian_gasket_01w_500.jpg apollonian_gasket_01_1000.jpg  apollonian_gasket_01_250.jpg apollonian_gasket_02_1000.jpg  apollonian_gasket_02_250.jpg apollonian_gasket_2D_01_1000.jpg  apollonian_gasket_2D_01_250.jpg apollonian_gasket_2D_02_1000.jpg apollonian_gasket_2D_02_250.jpg apollonian_gasket_03_1000.jpg apollonian_gasket_03_350.jpg

clean :
    rm -f *.jpg *.png *.inc *.svg *~ *.bak

#-----------------------------------------------------------------------------------------------------------------------------------
apf1.inc : agen.rb
    time ruby agen.rb --doSVG F --doPOV T --output apf1 --r1 1.0    --r2 1.0    --r3 1.0               --crvLim 20000.0 --povClipRad 0.5 --povMinRad 0.003

apf2.inc : agen.rb
    time ruby agen.rb --doSVG F --doPOV T --output apf2 --r1 1.0    --r2 1.0    --r3 0.66666666666666  --crvLim 20000.0 --povClipRad 0.5 --povMinRad 0.003

apf3.inc : agen.rb
    time ruby agen.rb --doSVG F --doPOV T --output apf3 --r1 1.0    --r2 0.6    --r3 0.6 --crvLim 20000.0 --povClipRad 0.5 --povMinRad 0.003

apf3.png : ap3.pov apf3.inc
    povray -W$(WRES)$(RESM) -H$(HRES)$(RESM) -Q11 -P -D +A0.4 -AM2 -R5 +J3 -Oapf3.png -HIapf3.inc -Iap3.pov

apf1.png : ap.pov apf1.inc
    povray -W$(WRES)$(RESM) -H$(HRES)$(RESM) -Q11 -P -D +A0.4 -AM2 -R5 +J3 -Oapf1.png -HIapf1.inc -Iap.pov

apw1.png : apw.pov apf1.inc
    povray -W$(WRES)$(RESM) -H$(HRES)$(RESM) -Q11 -P -D +A0.4 -AM2 -R5 +J3 -Oapw1.png -HIapf1.inc -Iapw.pov

apf2.png : ap3.pov apf2.inc
    povray -W$(WRES)$(RESM) -H$(HRES)$(RESM) -Q11 -P -D +A0.4 -AM2 -R5 +J3 -Oapf2.png -HIapf2.inc -Iap.pov

apollonian_gasket_01.png : apf1.png
    convert apf1.png -pointsize 100 -draw "gravity southeast fill white text 1,150 '©2017 Mitch Richling'" -pointsize 120 -draw "gravity northwest fill white text 1,10 'Apollonian Gasket'" -quality 100 apollonian_gasket_01.png

apollonian_gasket_01w.png : apw1.png
    convert apw1.png -pointsize 50 -draw "gravity southeast fill black text 200,300 '©2017 Mitch Richling'" -quality 100 apollonian_gasket_01w.png

apollonian_gasket_01w_1000.jpg : apollonian_gasket_01w.png
    convert -resize 1000 apollonian_gasket_01w.png apollonian_gasket_01w_1000.jpg

apollonian_gasket_01w_500.jpg : apollonian_gasket_01w.png
    convert -resize 500 apollonian_gasket_01w.png apollonian_gasket_01w_500.jpg

apollonian_gasket_01_1000.jpg : apollonian_gasket_01.png
    convert -resize 1000 apollonian_gasket_01.png apollonian_gasket_01_1000.jpg

apollonian_gasket_01_250.jpg : apollonian_gasket_01.png
    convert -resize 250 apollonian_gasket_01.png apollonian_gasket_01_250.jpg

apollonian_gasket_02_1000.jpg : apollonian_gasket_02.png
    convert -resize 1000 apollonian_gasket_02.png apollonian_gasket_02_1000.jpg

apollonian_gasket_02_250.jpg : apollonian_gasket_02.png
    convert -resize 250 apollonian_gasket_02.png apollonian_gasket_02_250.jpg

apollonian_gasket_03_1000.jpg : apollonian_gasket_03.png
    convert -resize 1000 apollonian_gasket_03.png apollonian_gasket_03_1000.jpg

apollonian_gasket_03_350.jpg : apollonian_gasket_03.png
    convert -resize 250 apollonian_gasket_03.png apollonian_gasket_03_250.jpg

apollonian_gasket_03.png : apf3.png
    convert apf3.png -pointsize 100 -draw "gravity southeast fill white text 1,150 '©2017 Mitch Richling'" -pointsize 120 -draw "gravity northwest fill white text 1,10 'Apollonian Gasket'" -quality 100 apollonian_gasket_03.png


apollonian_gasket_02.png : apf2.png
    convert apf2.png -pointsize 100 -draw "gravity southeast fill white text 1,150 '©2017 Mitch Richling'" -pointsize 120 -draw "gravity northwest fill white text 1,10 'Apollonian Gasket'" -quality 100 apollonian_gasket_02.png

apf1.svg : agen.rb
    time ruby agen.rb --doSVG T --doPOV F --output apf1 --r1 1000.0    --r2 1000.0    --r3 1000.0               --crvLim 1.0

apf1sw.png : apf1.svg
    convert -background white -flatten apf1.svg apf1sw.png

apf1sb.png : apf1.svg
    convert -background black -flatten apf1.svg apf1sb.png

apollonian_gasket_2D_01_1000.jpg : apf1sb.png
    convert -resize 1000 apf1sb.png apollonian_gasket_2D_01_1000.jpg

apollonian_gasket_2D_01_250.jpg : apf1sw.png
    convert -resize 250 apf1sw.png apollonian_gasket_2D_01_250.jpg

apf2.svg : agen.rb
    time ruby agen.rb --doSVG T --doPOV F --output apf2 --r1 1000.0    --r2 1000.0    --r3 666.6666666666       --crvLim 1.0

apf2sw.png : apf2.svg
    convert -background white -flatten apf2.svg apf2sw.png

apf2sb.png : apf2.svg
    convert -background black -flatten apf2.svg apf2sb.png

apollonian_gasket_2D_02_1000.jpg : apf2sb.png
    convert -resize 1000 apf2sb.png apollonian_gasket_2D_02_1000.jpg

apollonian_gasket_2D_02_250.jpg : apf2sw.png
    convert -resize 250 apf2sw.png apollonian_gasket_2D_02_250.jpg

4. An IFS Approach

Above we have taken a constructive approach using a direct description of the circles for 2D images of the Apollonian Gasket. This is not the only way! One interesting alternative is an IFS frequently attributed to Kravchenko Alexei and Mekhontsev Dmitriy (I don't have a reference).

\[\begin{array}{rcl} f_1(z) &=& f(z) \\ f_2(z) &=& \frac{-1 + i\sqrt{3}}{2f(z)} \\ f_3(z) &=& \frac{-1 - i\sqrt{3}}{2f(z)} \\ \end{array}\]

\[ f(z) = \frac{3}{1-z+\sqrt{3}} - \frac{1+\sqrt{3}}{2+\sqrt{3}} \]

We can generate an image with a bit of code – C++ this time:

#include "ramCanvas.hpp"

int main() {
  std::chrono::time_point<std::chrono::system_clock> startTime = std::chrono::system_clock::now();
  std::cout << "apollony start" << std::endl;
  const int   CSIZE = 1080*2;          // Quad HD
  const int ITRTOSS = static_cast<int>(std::pow(2, 10)); // Throw away first iterations
  const long NUMITR = static_cast<int>(std::pow(2, 29)); // Needs to be big

  mjr::ramCanvas3c8b theRamCanvas(CSIZE, CSIZE, -4.0, 4.0, -4.0, 4.0);
  theRamCanvas.setDrawMode(mjr::ramCanvas3c8b::drawModeType::ADDCLAMP);

#if SUPPORT_DRAWING_MODE
  mjr::ramCanvas3c8b::colorType aColor[] = { mjr::ramCanvas3c8b::colorType(1, 0, 0),
                                             mjr::ramCanvas3c8b::colorType(0, 1, 0), 
                                             mjr::ramCanvas3c8b::colorType(0, 0, 1) };
#else
  mjr::ramCanvas3c8b::colorType aColor[] = { mjr::ramCanvas3c8b::colorType(255, 0, 0), 
                                             mjr::ramCanvas3c8b::colorType(0, 255, 0), 
                                             mjr::ramCanvas3c8b::colorType(0, 0, 255) };
#endif

  std::random_device rd;
  std::minstd_rand0 rEng(rd()); // Fast is better than high quality for this application.

  const double s = 1.73205080757;
//  const std::complex<double> i(0, 1);
  const std::complex<double> si(0, s);
  const std::complex<double> c1 = (1+s)/(2+s);
  const std::complex<double> c2 =  0.5*(si-1.0);
  const std::complex<double> c3 = -0.5*(si+1.0);
  const std::complex<double> c4 = 1.0+s;
  const std::complex<double> c5(3.0, 0.0);
  std::complex<double> z(0.1, 0.2);
  for (long n=0;n<NUMITR;n++) {
    std::complex<double> zNxt;
    if ((n % (NUMITR/100)) == 0) {
      if ((n % (NUMITR/10)) == 0)
        std::cout << "|" << std::flush;
      else
        std::cout << "." << std::flush;
    }
    std::complex<double> f = c5/(c4-z)-c1;
    std::minstd_rand0::result_type rn = rEng()%3;
    switch (rn) {
      case 0:
        zNxt = f;    break;
      case 1:
        zNxt = c2/f; break;
      case 2:
        zNxt = c3/f; break;
    }
    z = zNxt;
    if (n > ITRTOSS)
      theRamCanvas.drawPoint(z, aColor[rn]);
  }
  std::cout << "|" << std::endl;
  std::cout << "apollony dump" << std::endl;
  theRamCanvas.applyHomoPixTfrm(&mjr::ramCanvas3c8b::colorType::tfrmStdPow, 1/5.0);
  theRamCanvas.writeTIFFfile("apollony.tiff");
  std::cout << "apollony finish" << std::endl;
  std::chrono::duration<double> runTime = std::chrono::system_clock::now() - startTime;
  std::cout << "Total Runtime " << runTime.count() << " sec" << std::endl;
  return 0;
}

The results

5. References

All the code (both C++ and POV-Ray) used to generate everything on this page may be found on github.

Check out the fractals section of my reading list.