The mandelbulb fractal shader from https://www.shadertoy.com/view/MdfGRr ported to C++ with GLUT. Uses W/A/S/D to rotate and up/down/left/right to move. Compiles on OS X with: g++ main.cpp gl.cpp -framework GLUT -framework OpenGL -Wno-deprecated-declarations

gl.cpp

#include "gl.h"

mat4 &mat4::transpose() {
    std::swap(m01, m10); std::swap(m02, m20); std::swap(m03, m30);
    std::swap(m12, m21); std::swap(m13, m31); std::swap(m23, m32);
    return *this;
}

mat4 &mat4::rotateX(float degrees) {
    float radians = degrees * (M_PI / 180);
    float s = sinf(radians), c = cosf(radians);
    float t01 = m01, t02 = m02;
    float t11 = m11, t12 = m12;
    float t21 = m21, t22 = m22;
    float t31 = m31, t32 = m32;
    m01 = c * t01 - s * t02;
    m02 = c * t02 + s * t01;
    m11 = c * t11 - s * t12;
    m12 = c * t12 + s * t11;
    m21 = c * t21 - s * t22;
    m22 = c * t22 + s * t21;
    m31 = c * t31 - s * t32;
    m32 = c * t32 + s * t31;
    return *this;
}

mat4 &mat4::rotateY(float degrees) {
    float radians = degrees * (M_PI / 180);
    float s = sinf(radians), c = cosf(radians);
    float t02 = m02, t00 = m00;
    float t12 = m12, t10 = m10;
    float t22 = m22, t20 = m20;
    float t32 = m32, t30 = m30;
    m02 = c * t02 - s * t00;
    m00 = c * t00 + s * t02;
    m12 = c * t12 - s * t10;
    m10 = c * t10 + s * t12;
    m22 = c * t22 - s * t20;
    m20 = c * t20 + s * t22;
    m32 = c * t32 - s * t30;
    m30 = c * t30 + s * t32;
    return *this;
}

mat4 &mat4::rotateZ(float degrees) {
    float radians = degrees * (M_PI / 180);
    float s = sinf(radians), c = cosf(radians);
    float t00 = m00, t01 = m01;
    float t10 = m10, t11 = m11;
    float t20 = m20, t21 = m21;
    float t30 = m30, t31 = m31;
    m00 = c * t00 - s * t01;
    m01 = c * t01 + s * t00;
    m10 = c * t10 - s * t11;
    m11 = c * t11 + s * t10;
    m20 = c * t20 - s * t21;
    m21 = c * t21 + s * t20;
    m30 = c * t30 - s * t31;
    m31 = c * t31 + s * t30;
    return *this;
}

mat4 &mat4::rotate(float degrees, float x, float y, float z) {
    float radians = degrees * (M_PI / 180);
    float d = sqrtf(x*x + y*y + z*z);
    float s = sinf(radians), c = cosf(radians), t = 1 - c;
    x /= d; y /= d; z /= d;
    return *this *= mat4(
        x*x*t + c, x*y*t - z*s, x*z*t + y*s, 0,
        y*x*t + z*s, y*y*t + c, y*z*t - x*s, 0,
        z*x*t - y*s, z*y*t + x*s, z*z*t + c, 0,
        0, 0, 0, 1);
}

mat4 &mat4::scale(float x, float y, float z) {
    m00 *= x; m01 *= y; m02 *= z;
    m10 *= x; m11 *= y; m12 *= z;
    m20 *= x; m21 *= y; m22 *= z;
    m30 *= x; m31 *= y; m32 *= z;
    return *this;
}

mat4 &mat4::translate(float x, float y, float z) {
    m03 += m00 * x + m01 * y + m02 * z;
    m13 += m10 * x + m11 * y + m12 * z;
    m23 += m20 * x + m21 * y + m22 * z;
    m33 += m30 * x + m31 * y + m32 * z;
    return *this;
}

mat4 &mat4::ortho(float l, float r, float b, float t, float n, float f) {
    return *this *= mat4(
        2 / (r - l), 0, 0, (r + l) / (l - r),
        0, 2 / (t - b), 0, (t + b) / (b - t),
        0, 0, 2 / (n - f), (f + n) / (n - f),
        0, 0, 0, 1);
}

mat4 &mat4::frustum(float l, float r, float b, float t, float n, float f) {
    return *this *= mat4(
        2 * n / (r - l), 0, (r + l) / (r - l), 0,
        0, 2 * n / (t - b), (t + b) / (t - b), 0,
        0, 0, (f + n) / (n - f), 2 * f * n / (n - f),
        0, 0, -1, 0);
}

mat4 &mat4::perspective(float fov, float aspect, float near, float far) {
    float y = tanf(fov * M_PI / 360) * near, x = y * aspect;
    return frustum(-x, x, -y, y, near, far);
}

mat4 &mat4::invert() {
    float t00 = m00, t01 = m01, t02 = m02, t03 = m03;
    *this = mat4(
        m11*m22*m33 - m11*m32*m23 - m12*m21*m33 + m12*m31*m23 + m13*m21*m32 - m13*m31*m22,
        -m01*m22*m33 + m01*m32*m23 + m02*m21*m33 - m02*m31*m23 - m03*m21*m32 + m03*m31*m22,
        m01*m12*m33 - m01*m32*m13 - m02*m11*m33 + m02*m31*m13 + m03*m11*m32 - m03*m31*m12,
        -m01*m12*m23 + m01*m22*m13 + m02*m11*m23 - m02*m21*m13 - m03*m11*m22 + m03*m21*m12,

        -m10*m22*m33 + m10*m32*m23 + m12*m20*m33 - m12*m30*m23 - m13*m20*m32 + m13*m30*m22,
        m00*m22*m33 - m00*m32*m23 - m02*m20*m33 + m02*m30*m23 + m03*m20*m32 - m03*m30*m22,
        -m00*m12*m33 + m00*m32*m13 + m02*m10*m33 - m02*m30*m13 - m03*m10*m32 + m03*m30*m12,
        m00*m12*m23 - m00*m22*m13 - m02*m10*m23 + m02*m20*m13 + m03*m10*m22 - m03*m20*m12,

        m10*m21*m33 - m10*m31*m23 - m11*m20*m33 + m11*m30*m23 + m13*m20*m31 - m13*m30*m21,
        -m00*m21*m33 + m00*m31*m23 + m01*m20*m33 - m01*m30*m23 - m03*m20*m31 + m03*m30*m21,
        m00*m11*m33 - m00*m31*m13 - m01*m10*m33 + m01*m30*m13 + m03*m10*m31 - m03*m30*m11,
        -m00*m11*m23 + m00*m21*m13 + m01*m10*m23 - m01*m20*m13 - m03*m10*m21 + m03*m20*m11,

        -m10*m21*m32 + m10*m31*m22 + m11*m20*m32 - m11*m30*m22 - m12*m20*m31 + m12*m30*m21,
        m00*m21*m32 - m00*m31*m22 - m01*m20*m32 + m01*m30*m22 + m02*m20*m31 - m02*m30*m21,
        -m00*m11*m32 + m00*m31*m12 + m01*m10*m32 - m01*m30*m12 - m02*m10*m31 + m02*m30*m11,
        m00*m11*m22 - m00*m21*m12 - m01*m10*m22 + m01*m20*m12 + m02*m10*m21 - m02*m20*m11
    );
    float det = m00 * t00 + m10 * t01 + m20 * t02 + m30 * t03;
    for (int i = 0; i < 16; i++) m[i] /= det;
    return *this;
}

mat4 &mat4::operator *= (const mat4 &t) {
    *this = mat4(
        m00*t.m00 + m01*t.m10 + m02*t.m20 + m03*t.m30,
        m00*t.m01 + m01*t.m11 + m02*t.m21 + m03*t.m31,
        m00*t.m02 + m01*t.m12 + m02*t.m22 + m03*t.m32,
        m00*t.m03 + m01*t.m13 + m02*t.m23 + m03*t.m33,

        m10*t.m00 + m11*t.m10 + m12*t.m20 + m13*t.m30,
        m10*t.m01 + m11*t.m11 + m12*t.m21 + m13*t.m31,
        m10*t.m02 + m11*t.m12 + m12*t.m22 + m13*t.m32,
        m10*t.m03 + m11*t.m13 + m12*t.m23 + m13*t.m33,

        m20*t.m00 + m21*t.m10 + m22*t.m20 + m23*t.m30,
        m20*t.m01 + m21*t.m11 + m22*t.m21 + m23*t.m31,
        m20*t.m02 + m21*t.m12 + m22*t.m22 + m23*t.m32,
        m20*t.m03 + m21*t.m13 + m22*t.m23 + m23*t.m33,

        m30*t.m00 + m31*t.m10 + m32*t.m20 + m33*t.m30,
        m30*t.m01 + m31*t.m11 + m32*t.m21 + m33*t.m31,
        m30*t.m02 + m31*t.m12 + m32*t.m22 + m33*t.m32,
        m30*t.m03 + m31*t.m13 + m32*t.m23 + m33*t.m33
    );
    return *this;
}

vec4 mat4::operator * (const vec4 &v) {
    return vec4(
        m00*v.x + m01*v.y + m02*v.z + m03*v.w,
        m10*v.x + m11*v.y + m12*v.z + m13*v.w,
        m20*v.x + m21*v.y + m22*v.z + m23*v.w,
        m30*v.x + m31*v.y + m32*v.z + m33*v.w
    );
}

vec4 operator * (const vec4 &v, const mat4 &t) {
    return vec4(
        t.m00*v.x + t.m10*v.y + t.m20*v.z + t.m30*v.w,
        t.m01*v.x + t.m11*v.y + t.m21*v.z + t.m31*v.w,
        t.m02*v.x + t.m12*v.y + t.m22*v.z + t.m32*v.w,
        t.m03*v.x + t.m13*v.y + t.m23*v.z + t.m33*v.w
    );
}

std::ostream &operator << (std::ostream &out, const mat4 &t) {
    return out << "mat4("
        << t.m00 << ", " << t.m01 << ", " << t.m02 << ", " << t.m03 << ",\n     "
        << t.m10 << ", " << t.m11 << ", " << t.m12 << ", " << t.m13 << ",\n     "
        << t.m20 << ", " << t.m21 << ", " << t.m22 << ", " << t.m23 << ",\n     "
        << t.m30 << ", " << t.m31 << ", " << t.m32 << ", " << t.m33 << ")";
}

Texture &Texture::create(int w, int h, int internalFormat, int format, int type, int filter, int wrap, void *data) {
    width = w;
    height = h;
    if (!id) glGenTextures(1, &id);
    bind();
    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, filter);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, filter);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, wrap);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, wrap);
    glTexImage2D(GL_TEXTURE_2D, 0, internalFormat, w, h, 0, format, type, data);
    unbind();
    return *this;
}

void Texture::swapWith(Texture &other) {
    std::swap(id, other.id);
    std::swap(width, other.width);
    std::swap(height, other.height);
}

void FBO::bind() {
    glBindFramebuffer(GL_FRAMEBUFFER, id);
    if (resizeViewport) {
        glGetIntegerv(GL_VIEWPORT, oldViewport);
        glViewport(newViewport[0], newViewport[1], newViewport[2], newViewport[3]);
    }
}

void FBO::unbind() {
    glBindFramebuffer(GL_FRAMEBUFFER, 0);
    if (resizeViewport) {
        glViewport(oldViewport[0], oldViewport[1], oldViewport[2], oldViewport[3]);
    }
}

FBO &FBO::attachColor(const Texture &texture, unsigned int attachment) {
    newViewport[2] = texture.width;
    newViewport[3] = texture.height;
    if (!id) glGenFramebuffers(1, &id);
    bind();

    // Bind a 2D texture
    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0 + attachment, GL_TEXTURE_2D, texture.id, 0);

    // Need to call glDrawBuffers() for OpenGL to draw to multiple attachments
    if (attachment >= drawBuffers.size()) drawBuffers.resize(attachment + 1, GL_NONE);
    drawBuffers[attachment] = GL_COLOR_ATTACHMENT0 + attachment;
    glDrawBuffers(drawBuffers.size(), drawBuffers.data());

    unbind();
    return *this;
}

FBO &FBO::detachColor(unsigned int attachment) {
    bind();

    // Update the draw buffers
    if (attachment < drawBuffers.size()) {
        drawBuffers[attachment] = GL_NONE;
        glDrawBuffers(drawBuffers.size(), drawBuffers.data());
    }

    unbind();
    return *this;
}

FBO &FBO::check() {
    bind();
    if (autoDepth) {
        if (renderbufferWidth != newViewport[2] || renderbufferHeight != newViewport[3]) {
            renderbufferWidth = newViewport[2];
            renderbufferHeight = newViewport[3];
            if (!renderbuffer) glGenRenderbuffers(1, &renderbuffer);
            glBindRenderbuffer(GL_RENDERBUFFER, renderbuffer);
            glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT32, renderbufferWidth, renderbufferHeight);
            glBindRenderbuffer(GL_RENDERBUFFER, 0);
        }
        glFramebufferRenderbuffer(GL_FRAMEBUFFER, GL_DEPTH_ATTACHMENT, GL_RENDERBUFFER, renderbuffer);
    }
    switch (glCheckFramebufferStatus(GL_FRAMEBUFFER)) {
        case GL_FRAMEBUFFER_COMPLETE: break;
        case GL_FRAMEBUFFER_UNDEFINED: printf("GL_FRAMEBUFFER_UNDEFINED\n"); exit(0);
        case GL_FRAMEBUFFER_UNSUPPORTED: printf("GL_FRAMEBUFFER_UNSUPPORTED\n"); exit(0);
        case GL_FRAMEBUFFER_INCOMPLETE_ATTACHMENT: printf("GL_FRAMEBUFFER_INCOMPLETE_ATTACHMENT\n"); exit(0);
        case GL_FRAMEBUFFER_INCOMPLETE_DRAW_BUFFER: printf("GL_FRAMEBUFFER_INCOMPLETE_DRAW_BUFFER\n"); exit(0);
        case GL_FRAMEBUFFER_INCOMPLETE_READ_BUFFER: printf("GL_FRAMEBUFFER_INCOMPLETE_READ_BUFFER\n"); exit(0);
        case GL_FRAMEBUFFER_INCOMPLETE_MULTISAMPLE: printf("GL_FRAMEBUFFER_INCOMPLETE_MULTISAMPLE\n"); exit(0);
        case GL_FRAMEBUFFER_INCOMPLETE_FORMATS_EXT: printf("GL_FRAMEBUFFER_INCOMPLETE_FORMATS_EXT\n"); exit(0);
        case GL_FRAMEBUFFER_INCOMPLETE_DIMENSIONS_EXT: printf("GL_FRAMEBUFFER_INCOMPLETE_DIMENSIONS_EXT\n"); exit(0);
        case GL_FRAMEBUFFER_INCOMPLETE_MISSING_ATTACHMENT: printf("GL_FRAMEBUFFER_INCOMPLETE_MISSING_ATTACHMENT\n"); exit(0);
        default: printf("Unknown glCheckFramebufferStatus error"); exit(0);
    }
    unbind();
    return *this;
}

Shader::~Shader() {
    glDeleteProgram(id);
    for (size_t i = 0; i < stages.size(); i++) {
        glDeleteShader(stages[i]);
    }
}

static void error(const char *type, const char *error, const char *source = NULL) {
    printf("----- %s -----\n", type);
    printf("%s\n", error);
    if (source) {
        printf("----- source code -----\n");
        printf("%s\n", source);
    }
    exit(0);
}

Shader &Shader::shader(int type, const char *source) {
    // Compile shader
    unsigned int shader = glCreateShader(type);
    glShaderSource(shader, 1, &source, NULL);
    glCompileShader(shader);
    stages.push_back(shader);

    // Check for errors
    char buffer[512];
    int length;
    glGetShaderInfoLog(shader, sizeof(buffer), &length, buffer);
    if (length) error("compile error", buffer, source);

    // Allow chaining
    return *this;
}

void Shader::link() {
    // Create and link program
    if (!id) id = glCreateProgram();
    for (size_t i = 0; i < stages.size(); i++) {
        glAttachShader(id, stages[i]);
    }
    glLinkProgram(id);

    // Check for errors
    char buffer[512];
    int length;
    glGetProgramInfoLog(id, sizeof(buffer), &length, buffer);
    if (length) error("link error", buffer);
}

void VAO::check() {
    if (vertices && vertices->currentTarget() != GL_ARRAY_BUFFER) {
        printf("expected vertices to have GL_ARRAY_BUFFER, got 0x%04X\n", vertices->currentTarget());
        exit(0);
    }
    if (indices && indices->currentTarget() != GL_ELEMENT_ARRAY_BUFFER) {
        printf("expected indices to have GL_ELEMENT_ARRAY_BUFFER, got 0x%04X\n", indices->currentTarget());
        exit(0);
    }
    if (offset != stride) {
        printf("expected size of attributes (%d bytes) to add up to size of vertex (%d bytes)\n", offset, stride);
        exit(0);
    }
}

gl.h

#ifndef GL_H
#define GL_H

#ifdef __APPLE_CC__
#define GL_DO_NOT_WARN_IF_MULTI_GL_VERSION_HEADERS_INCLUDED
#include <OpenGL/gl3.h>
#include <OpenGL/glu.h>
#else
#include <GL/gl.h>
#include <GL/glu.h>
#endif

#include <stdlib.h>
#include <stdio.h>
#include <ostream>
#include <vector>
#include <math.h>

struct vec2 {
    union {
        struct { float x, y; };
        struct { float s, t; };
        struct { float r, g; };
        float xy[2];
        float st[2];
        float rg[2];
    };

    vec2() : x(), y() {}
    vec2(float x, float y) : x(x), y(y) {}
    vec2(const vec2 &xy) : x(xy.x), y(xy.y) {}
    explicit vec2(float f) : x(f), y(f) {}

    vec2 operator + () const { return vec2(+x, +y); }
    vec2 operator - () const { return vec2(-x, -y); }

    vec2 operator + (const vec2 &vec) const { return vec2(x + vec.x, y + vec.y); }
    vec2 operator - (const vec2 &vec) const { return vec2(x - vec.x, y - vec.y); }
    vec2 operator * (const vec2 &vec) const { return vec2(x * vec.x, y * vec.y); }
    vec2 operator / (const vec2 &vec) const { return vec2(x / vec.x, y / vec.y); }
    vec2 operator + (float s) const { return vec2(x + s, y + s); }
    vec2 operator - (float s) const { return vec2(x - s, y - s); }
    vec2 operator * (float s) const { return vec2(x * s, y * s); }
    vec2 operator / (float s) const { return vec2(x / s, y / s); }

    friend vec2 operator + (float s, const vec2 &vec) { return vec2(s + vec.x, s + vec.y); }
    friend vec2 operator - (float s, const vec2 &vec) { return vec2(s - vec.x, s - vec.y); }
    friend vec2 operator * (float s, const vec2 &vec) { return vec2(s * vec.x, s * vec.y); }
    friend vec2 operator / (float s, const vec2 &vec) { return vec2(s / vec.x, s / vec.y); }

    vec2 &operator += (const vec2 &vec) { return *this = *this + vec; }
    vec2 &operator -= (const vec2 &vec) { return *this = *this - vec; }
    vec2 &operator *= (const vec2 &vec) { return *this = *this * vec; }
    vec2 &operator /= (const vec2 &vec) { return *this = *this / vec; }
    vec2 &operator += (float s) { return *this = *this + s; }
    vec2 &operator -= (float s) { return *this = *this - s; }
    vec2 &operator *= (float s) { return *this = *this * s; }
    vec2 &operator /= (float s) { return *this = *this / s; }

    bool operator == (const vec2 &vec) const { return x == vec.x && y == vec.y; }
    bool operator != (const vec2 &vec) const { return x != vec.x || y != vec.y; }

    friend float length(const vec2 &v) { return sqrtf(v.x * v.x + v.y * v.y); }
    friend float dot(const vec2 &a, const vec2 &b) { return a.x * b.x + a.y * b.y; }
    friend float max(const vec2 &v) { return fmaxf(v.x, v.y); }
    friend float min(const vec2 &v) { return fminf(v.x, v.y); }
    friend vec2 max(const vec2 &a, const vec2 &b) { return vec2(fmaxf(a.x, b.x), fmaxf(a.y, b.y)); }
    friend vec2 min(const vec2 &a, const vec2 &b) { return vec2(fminf(a.x, b.x), fminf(a.y, b.y)); }
    friend vec2 floor(const vec2 &v) { return vec2(floorf(v.x), floorf(v.y)); }
    friend vec2 ceil(const vec2 &v) { return vec2(ceilf(v.x), ceilf(v.y)); }
    friend vec2 abs(const vec2 &v) { return vec2(fabsf(v.x), fabsf(v.y)); }
    friend vec2 fract(const vec2 &v) { return v - floor(v); }
    friend vec2 normalized(const vec2 &v) { return v / length(v); }

    friend std::ostream &operator << (std::ostream &out, const vec2 &v) {
        return out << "vec2(" << v.x << ", " << v.y << ")";
    }
};

struct vec3 {
    union {
        struct { float x, y, z; };
        struct { float s, t, p; };
        struct { float r, g, b; };
        float xyz[3];
        float stp[3];
        float rgb[3];
    };

    vec3() : x(), y(), z() {}
    vec3(float x, float y, float z) : x(x), y(y), z(z) {}
    vec3(const vec2 &xy, float z) : x(xy.x), y(xy.y), z(z) {}
    vec3(float x, const vec2 &yz) : x(x), y(yz.x), z(yz.y) {}
    vec3(const vec3 &xyz) : x(xyz.x), y(xyz.y), z(xyz.z) {}
    explicit vec3(float f) : x(f), y(f), z(f) {}

    vec3 operator + () const { return vec3(+x, +y, +z); }
    vec3 operator - () const { return vec3(-x, -y, -z); }

    vec3 operator + (const vec3 &vec) const { return vec3(x + vec.x, y + vec.y, z + vec.z); }
    vec3 operator - (const vec3 &vec) const { return vec3(x - vec.x, y - vec.y, z - vec.z); }
    vec3 operator * (const vec3 &vec) const { return vec3(x * vec.x, y * vec.y, z * vec.z); }
    vec3 operator / (const vec3 &vec) const { return vec3(x / vec.x, y / vec.y, z / vec.z); }
    vec3 operator + (float s) const { return vec3(x + s, y + s, z + s); }
    vec3 operator - (float s) const { return vec3(x - s, y - s, z - s); }
    vec3 operator * (float s) const { return vec3(x * s, y * s, z * s); }
    vec3 operator / (float s) const { return vec3(x / s, y / s, z / s); }

    friend vec3 operator + (float s, const vec3 &vec) { return vec3(s + vec.x, s + vec.y, s + vec.z); }
    friend vec3 operator - (float s, const vec3 &vec) { return vec3(s - vec.x, s - vec.y, s - vec.z); }
    friend vec3 operator * (float s, const vec3 &vec) { return vec3(s * vec.x, s * vec.y, s * vec.z); }
    friend vec3 operator / (float s, const vec3 &vec) { return vec3(s / vec.x, s / vec.y, s / vec.z); }

    vec3 &operator += (const vec3 &vec) { return *this = *this + vec; }
    vec3 &operator -= (const vec3 &vec) { return *this = *this - vec; }
    vec3 &operator *= (const vec3 &vec) { return *this = *this * vec; }
    vec3 &operator /= (const vec3 &vec) { return *this = *this / vec; }
    vec3 &operator += (float s) { return *this = *this + s; }
    vec3 &operator -= (float s) { return *this = *this - s; }
    vec3 &operator *= (float s) { return *this = *this * s; }
    vec3 &operator /= (float s) { return *this = *this / s; }

    bool operator == (const vec3 &vec) const { return x == vec.x && y == vec.y && z == vec.z; }
    bool operator != (const vec3 &vec) const { return x != vec.x || y != vec.y || z != vec.z; }

    friend float length(const vec3 &v) { return sqrtf(v.x * v.x + v.y * v.y + v.z * v.z); }
    friend float dot(const vec3 &a, const vec3 &b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
    friend float max(const vec3 &v) { return fmaxf(fmaxf(v.x, v.y), v.z); }
    friend float min(const vec3 &v) { return fminf(fminf(v.x, v.y), v.z); }
    friend vec3 max(const vec3 &a, const vec3 &b) { return vec3(fmaxf(a.x, b.x), fmaxf(a.y, b.y), fmaxf(a.z, b.z)); }
    friend vec3 min(const vec3 &a, const vec3 &b) { return vec3(fminf(a.x, b.x), fminf(a.y, b.y), fminf(a.z, b.z)); }
    friend vec3 floor(const vec3 &v) { return vec3(floorf(v.x), floorf(v.y), floorf(v.z)); }
    friend vec3 ceil(const vec3 &v) { return vec3(ceilf(v.x), ceilf(v.y), ceilf(v.z)); }
    friend vec3 abs(const vec3 &v) { return vec3(fabsf(v.x), fabsf(v.y), fabsf(v.z)); }
    friend vec3 fract(const vec3 &v) { return v - floor(v); }
    friend vec3 normalized(const vec3 &v) { return v / length(v); }
    friend vec3 cross(const vec3 &a, const vec3 &b) { return vec3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x); }

    friend std::ostream &operator << (std::ostream &out, const vec3 &v) {
        return out << "vec3(" << v.x << ", " << v.y << ", " << v.z << ")";
    }
};

struct vec4 {
    union {
        struct { float x, y, z, w; };
        struct { float s, t, p, q; };
        struct { float r, g, b, a; };
        float xyzw[4];
        float stpq[4];
        float rgba[4];
    };

    vec4() : x(), y(), z(), w() {}
    vec4(float x, float y, float z, float w) : x(x), y(y), z(z), w(w) {}
    vec4(const vec2 &xy, float z, float w) : x(xy.x), y(xy.y), z(z), w(w) {}
    vec4(float x, const vec2 &yz, float w) : x(x), y(yz.x), z(yz.y), w(w) {}
    vec4(float x, float y, const vec2 &zw) : x(x), y(y), z(zw.x), w(zw.y) {}
    vec4(const vec2 &xy, const vec2 &zw) : x(xy.x), y(xy.y), z(zw.x), w(zw.y) {}
    vec4(const vec3 &xyz, float w) : x(xyz.x), y(xyz.y), z(xyz.z), w(w) {}
    vec4(float x, const vec3 &yzw) : x(x), y(yzw.x), z(yzw.y), w(yzw.z) {}
    vec4(const vec4 &xyzw) : x(xyzw.x), y(xyzw.y), z(xyzw.z), w(xyzw.w) {}
    explicit vec4(float f) : x(f), y(f), z(f), w(f) {}

    vec4 operator + () const { return vec4(+x, +y, +z, +w); }
    vec4 operator - () const { return vec4(-x, -y, -z, -w); }

    vec4 operator + (const vec4 &vec) const { return vec4(x + vec.x, y + vec.y, z + vec.z, w + vec.w); }
    vec4 operator - (const vec4 &vec) const { return vec4(x - vec.x, y - vec.y, z - vec.z, w - vec.w); }
    vec4 operator * (const vec4 &vec) const { return vec4(x * vec.x, y * vec.y, z * vec.z, w * vec.w); }
    vec4 operator / (const vec4 &vec) const { return vec4(x / vec.x, y / vec.y, z / vec.z, w / vec.w); }
    vec4 operator + (float s) const { return vec4(x + s, y + s, z + s, w + s); }
    vec4 operator - (float s) const { return vec4(x - s, y - s, z - s, w - s); }
    vec4 operator * (float s) const { return vec4(x * s, y * s, z * s, w * s); }
    vec4 operator / (float s) const { return vec4(x / s, y / s, z / s, w / s); }

    friend vec4 operator + (float s, const vec4 &vec) { return vec4(s + vec.x, s + vec.y, s + vec.z, s + vec.w); }
    friend vec4 operator - (float s, const vec4 &vec) { return vec4(s - vec.x, s - vec.y, s - vec.z, s - vec.w); }
    friend vec4 operator * (float s, const vec4 &vec) { return vec4(s * vec.x, s * vec.y, s * vec.z, s * vec.w); }
    friend vec4 operator / (float s, const vec4 &vec) { return vec4(s / vec.x, s / vec.y, s / vec.z, s / vec.w); }

    vec4 &operator += (const vec4 &vec) { return *this = *this + vec; }
    vec4 &operator -= (const vec4 &vec) { return *this = *this - vec; }
    vec4 &operator *= (const vec4 &vec) { return *this = *this * vec; }
    vec4 &operator /= (const vec4 &vec) { return *this = *this / vec; }
    vec4 &operator += (float s) { return *this = *this + s; }
    vec4 &operator -= (float s) { return *this = *this - s; }
    vec4 &operator *= (float s) { return *this = *this * s; }
    vec4 &operator /= (float s) { return *this = *this / s; }

    bool operator == (const vec4 &vec) const { return x == vec.x && y == vec.y && z == vec.z && w == vec.w; }
    bool operator != (const vec4 &vec) const { return x != vec.x || y != vec.y || z != vec.z || w != vec.w; }

    friend float length(const vec4 &v) { return sqrtf(v.x * v.x + v.y * v.y + v.z * v.z + v.w * v.w); }
    friend float dot(const vec4 &a, const vec4 &b) { return a.x * b.x + a.y * b.y + a.z * b.z + a.w * a.w; }
    friend float max(const vec4 &v) { return fmaxf(fmaxf(v.x, v.y), fmaxf(v.z, v.w)); }
    friend float min(const vec4 &v) { return fminf(fminf(v.x, v.y), fminf(v.z, v.w)); }
    friend vec4 max(const vec4 &a, const vec4 &b) { return vec4(fmaxf(a.x, b.x), fmaxf(a.y, b.y), fmaxf(a.z, b.z), fmaxf(a.w, b.w)); }
    friend vec4 min(const vec4 &a, const vec4 &b) { return vec4(fminf(a.x, b.x), fminf(a.y, b.y), fminf(a.z, b.z), fminf(a.w, b.w)); }
    friend vec4 floor(const vec4 &v) { return vec4(floorf(v.x), floorf(v.y), floorf(v.z), floorf(v.w)); }
    friend vec4 ceil(const vec4 &v) { return vec4(ceilf(v.x), ceilf(v.y), ceilf(v.z), ceilf(v.w)); }
    friend vec4 abs(const vec4 &v) { return vec4(fabsf(v.x), fabsf(v.y), fabsf(v.z), fabsf(v.w)); }
    friend vec4 fract(const vec4 &v) { return v - floor(v); }
    friend vec4 normalized(const vec4 &v) { return v / length(v); }

    friend std::ostream &operator << (std::ostream &out, const vec4 &v) {
        return out << "vec4(" << v.x << ", " << v.y << ", " << v.z << ", " << v.w << ")";
    }
};

struct mat4 {
    union {
        struct {
            float m00, m01, m02, m03;
            float m10, m11, m12, m13;
            float m20, m21, m22, m23;
            float m30, m31, m32, m33;
        };
        float m[16];
    };

    mat4() :
        m00(1), m01(), m02(), m03(),
        m10(), m11(1), m12(), m13(),
        m20(), m21(), m22(1), m23(),
        m30(), m31(), m32(), m33(1) {}
    mat4(const vec4 &r0, const vec4 &r1, const vec4 &r2, const vec4 &r3) :
        m00(r0.x), m01(r0.y), m02(r0.z), m03(r0.w),
        m10(r1.x), m11(r1.y), m12(r1.z), m13(r1.w),
        m20(r2.x), m21(r2.y), m22(r2.z), m23(r2.w),
        m30(r3.x), m31(r3.y), m32(r3.z), m33(r3.w) {}
    mat4(
        float m00, float m01, float m02, float m03,
        float m10, float m11, float m12, float m13,
        float m20, float m21, float m22, float m23,
        float m30, float m31, float m32, float m33) :
        m00(m00), m01(m01), m02(m02), m03(m03),
        m10(m10), m11(m11), m12(m12), m13(m13),
        m20(m20), m21(m21), m22(m22), m23(m23),
        m30(m30), m31(m31), m32(m32), m33(m33) {}

    mat4 &transpose();
    mat4 &rotateX(float degrees);
    mat4 &rotateY(float degrees);
    mat4 &rotateZ(float degrees);
    mat4 &rotate(float degrees, float x, float y, float z);
    mat4 &rotate(float degrees, const vec3 &v) { return rotate(degrees, v.x, v.y, v.z); }
    mat4 &scale(float x, float y, float z);
    mat4 &scale(const vec3 &v) { return scale(v.x, v.y, v.z); }
    mat4 &translate(float x, float y, float z);
    mat4 &translate(const vec3 &v) { return translate(v.x, v.y, v.z); }
    mat4 &ortho(float l, float r, float b, float t, float n, float f);
    mat4 &frustum(float l, float r, float b, float t, float n, float f);
    mat4 &perspective(float fov, float aspect, float near, float far);
    mat4 &invert();

    mat4 &operator *= (const mat4 &t);
    vec4 operator * (const vec4 &v);
    friend vec4 operator * (const vec4 &v, const mat4 &t);
    friend std::ostream &operator << (std::ostream &out, const mat4 &t);
};

// When rendering back and forth between two textures (ping-ponging), it is
// easiest to just call swapWith() after rendering.
struct Texture {
    unsigned int id;
    int width, height;

    Texture() : id(), width(), height() {}
    ~Texture() { glDeleteTextures(1, &id); }

    void bind(int unit = 0) const { glActiveTexture(GL_TEXTURE0 + unit); glBindTexture(GL_TEXTURE_2D, id); }
    void unbind(int unit = 0) const { glActiveTexture(GL_TEXTURE0 + unit); glBindTexture(GL_TEXTURE_2D, 0); }

    Texture &create(int width, int height, int internalFormat, int format, int type, int filter, int wrap, void *data = NULL);

    // Swap the members of this texture with the members of other.
    void swapWith(Texture &other);
};

// A framebuffer object that can take color attachments. Draw calls between
// bind() and unbind() are drawn to the attached textures.
//
// Usage:
//
//     FBO fbo;
//
//     fbo.attachColor(texture).check().bind();
//     // draw stuff
//     fbo.unbind();
//
struct FBO {
    unsigned int id;
    unsigned int renderbuffer;
    bool autoDepth;
    bool resizeViewport;
    int newViewport[4], oldViewport[4];
    int renderbufferWidth, renderbufferHeight;
    std::vector<unsigned int> drawBuffers;

    FBO(bool autoDepth = true, bool resizeViewport = true) : id(), renderbuffer(), autoDepth(autoDepth),
        resizeViewport(resizeViewport), newViewport(), oldViewport(), renderbufferWidth(), renderbufferHeight() {}
    ~FBO() { glDeleteFramebuffers(1, &id); glDeleteRenderbuffers(1, &renderbuffer); }

    // Draw calls between these will be drawn to attachments. If resizeViewport
    // is true this will automatically resize the viewport to the size of the
    // last attached texture.
    void bind();
    void unbind();

    // Draw to texture 2D in the indicated attachment location (or a 2D layer of
    // a 3D texture).
    FBO &attachColor(const Texture &texture, unsigned int attachment = 0);

    // Stop drawing to the indicated color attachment
    FBO &detachColor(unsigned int attachment = 0);

    // Call after all attachColor() calls, validates attachments.
    FBO &check();
};

// Use this macro to pass raw GLSL to Shader::shader()
#define glsl(x) #x

// Wraps a GLSL shader program and all attached shader stages. Meant to be used
// with the glsl() macro.
//
// Usage:
//
//     // Initialization
//     Shader shader;
//     shader.vertexShader(glsl(
//         attribute vec4 vertex;
//         void main() {
//             gl_Position = vertex;
//         }
//     )).fragmentShader(glsl(
//         void main() {
//             gl_FragColor = vec4(1.0, 0.0, 0.0, 1.0);
//         }
//     )).link();
//
//     // Rendering
//     shader.use();
//     // Draw stuff
//     shader.unuse();
//
struct Shader {
    unsigned int id;
    std::vector<unsigned int> stages;

    Shader() : id() {}
    ~Shader();

    Shader &shader(int type, const char *source);
    Shader &vertexShader(const char *source) { return shader(GL_VERTEX_SHADER, source); }
    Shader &fragmentShader(const char *source) { return shader(GL_FRAGMENT_SHADER, source); }

    void link();
    void use() const { glUseProgram(id); }
    void unuse() const { glUseProgram(0); }

    unsigned int attribute(const char *name) const { return glGetAttribLocation(id, name); }
    unsigned int uniform(const char *name) const { return glGetUniformLocation(id, name); }

    void uniformInt(const char *name, int i) const { glUniform1i(uniform(name), i); }
    void uniformFloat(const char *name, float f) const { glUniform1f(uniform(name), f); }
    void uniform(const char *name, const vec2 &v) const { glUniform2fv(uniform(name), 1, v.xy); }
    void uniform(const char *name, const vec3 &v) const { glUniform3fv(uniform(name), 1, v.xyz); }
    void uniform(const char *name, const vec4 &v) const { glUniform4fv(uniform(name), 1, v.xyzw); }
    void uniform(const char *name, const mat4 &m) const { glUniformMatrix4fv(uniform(name), 1, true, m.m); }
};

// A vertex buffer containing a certain type. For example, the simplest vertex
// buffer would be Buffer<vec3> but more complex formats could use Buffer<Vertex>
// with struct Vertex { vec3 position; vec2 coord; }. Each Buffer is meant to be
// used with a VAO.
//
// Index buffers should be specified using an integer type directly, i.e.
// Buffer<unsigned short>, since that's what VAO expects.
//
// Usage:
//
//     Buffer<vec3> vertices;
//     vertices << vec3(1, 0, 0) << vec3(0, 1, 0) << vec3(0, 0, 1);
//     vertices.upload();
//
//     Buffer<char> indices;
//     indices << 0 << 1 << 2;
//     indices.upload(GL_ELEMENT_ARRAY_BUFFER);
//
template <typename T>
struct Buffer {
    std::vector<T> data;
    unsigned int id;
    int currentTarget;

    Buffer() : id(), currentTarget() {}
    ~Buffer() {
        glDeleteBuffers(1, &id);
    }

    void bind() const { glBindBuffer(currentTarget, id); }
    void unbind() const { glBindBuffer(currentTarget, 0); }

    void upload(int target = GL_ARRAY_BUFFER, int usage = GL_STATIC_DRAW) {
        if (!id) glGenBuffers(1, &id);
        currentTarget = target;
        bind();
        glBufferData(currentTarget, data.size() * sizeof(T), data.data(), usage);
        unbind();
    }

    unsigned int size() const { return data.size(); }
    Buffer<T> &operator << (const T &t) { data.push_back(t); return *this; }
};

// Convert a C++ type to an OpenGL type enum using TypeToOpenGL<T>::value
template <typename T> struct TypeToOpenGL {};
template <> struct TypeToOpenGL<bool> { enum { value = GL_BOOL }; };
template <> struct TypeToOpenGL<float> { enum { value = GL_FLOAT }; };
template <> struct TypeToOpenGL<double> { enum { value = GL_DOUBLE }; };
template <> struct TypeToOpenGL<int> { enum { value = GL_INT }; };
template <> struct TypeToOpenGL<char> { enum { value = GL_BYTE }; };
template <> struct TypeToOpenGL<short> { enum { value = GL_SHORT }; };
template <> struct TypeToOpenGL<unsigned int> { enum { value = GL_UNSIGNED_INT }; };
template <> struct TypeToOpenGL<unsigned char> { enum { value = GL_UNSIGNED_BYTE }; };
template <> struct TypeToOpenGL<unsigned short> { enum { value = GL_UNSIGNED_SHORT }; };

// Groups a Buffer for vertices with a set of attributes for drawing. Can
// optionally include a Buffer for indices too.
//
// Usage:
//
//     // Initialization
//     VAO vao;
//     vao.create(shader, vertices, indices).attribute<float>("vertices", 3).check();
//
//     // Rendering
//     vao.draw(GL_TRIANGLE_STRIP);
//
struct VAO {
    // A holder for a Buffer so we can query information (i.e. number of vertices
    // for drawing). Type erasure is used so the VAO class doesn't need to be
    // templated.
    struct BufferHolder {
        virtual ~BufferHolder() {}
        virtual int currentTarget() const = 0;
        virtual unsigned int size() const = 0;
    };
    template <typename T>
    struct BufferHolderImpl : BufferHolder {
        const Buffer<T> &buffer;
        BufferHolderImpl(const Buffer<T> &buffer) : buffer(buffer) {}
        int currentTarget() const { return buffer.currentTarget; }
        unsigned int size() const { return buffer.size(); }
    };

    // You should not need to access these
    unsigned int id;
    int stride, offset, indexType;
    const Shader *shader;
    const BufferHolder *vertices;
    const BufferHolder *indices;

    VAO() : id(), stride(), offset(), indexType(), shader(), vertices(), indices() {}
    ~VAO() { glDeleteVertexArrays(1, &id); delete vertices; delete indices; }

    // You should not need to bind a VAO directly
    void bind() const { glBindVertexArray(id); }
    void unbind() const { glBindVertexArray(0); }

    // Create a vertex array object referencing a shader and a vertex buffer.
    // The shader is used to query the location of attributes in attribute()
    // and the vertex buffer is used to determine the number of elements to
    // draw in draw() and drawInstanced().
    template <typename Vertex>
    VAO &create(const Shader &shader, const Buffer<Vertex> &vbo) {
        delete vertices;
        delete indices;

        this->shader = &shader;
        vertices = new BufferHolderImpl<Vertex>(vbo);
        indices = NULL;
        stride = sizeof(Vertex);
        indexType = GL_INVALID_ENUM;

        if (!id) glGenVertexArrays(1, &id);
        bind();
        vbo.bind();
        unbind();

        return *this;
    }

    // Create a vertex array object referencing a shader, a vertex buffer, and
    // an index buffer. The shader is used to query the location of attributes
    // in attribute() and the index buffer is used to determine the number of
    // elements to draw in draw() and drawInstanced().
    template <typename Vertex, typename Index>
    VAO &create(const Shader &shader, const Buffer<Vertex> &vbo, const Buffer<Index> &ibo) {
        delete vertices;
        delete indices;

        this->shader = &shader;
        vertices = new BufferHolderImpl<Vertex>(vbo);
        indices = new BufferHolderImpl<Index>(ibo);
        stride = sizeof(Vertex);
        indexType = TypeToOpenGL<Index>::value;

        if (!id) glGenVertexArrays(1, &id);
        bind();
        vbo.bind();
        ibo.bind();
        unbind();

        return *this;
    }

    // Define an attribute called name in the provided shader. This attribute
    // has count elements of type T. If normalized is true, integer types are
    // mapped from 0 to 2^n-1 (where n is the bit count) to values from 0 to 1.
    //
    // Attributes should be declared in the order they are declared in the
    // vertex struct (assuming interleaved data). Call check() after declaring
    // all attributes to make sure the vertex struct is packed.
    template <typename T>
    VAO &attribute(const char *name, int count, bool normalized = false) {
        int location = shader->attribute(name);
        bind();
        glEnableVertexAttribArray(location);
        glVertexAttribPointer(location, count, TypeToOpenGL<T>::value, normalized, stride, (char *)NULL + offset);
        unbind();
        offset += count * sizeof(T);
        return *this;
    }

    // Validate VBO modes and attribute byte sizes
    void check();

    // Draw the attached VBOs
    void draw(int mode = GL_TRIANGLES) const {
        bind();
        if (indices) glDrawElements(mode, indices->size(), indexType, NULL);
        else glDrawArrays(mode, 0, vertices->size());
        unbind();
    }
};

#endif // GL_H

main.cpp

#ifdef __APPLE_CC__
#include <GLUT/glut.h>
#else
#include <GL/glut.h>
#endif
#include "gl.h"

Shader shader;
Buffer<vec2> quad;
VAO layout;

void setup() {
    shader.vertexShader(glsl(
      attribute vec2 vertex;
      varying vec2 oTexCoord;

      void main() {
        oTexCoord = 0.5 + 0.5 * vertex;
        gl_Position = vec4(vertex, 0.0, 1.0);
      }
   )).fragmentShader(glsl(
      uniform vec3 eye;
      uniform vec3 axisX;
      uniform vec3 axisY;
      uniform vec3 axisZ;
      varying vec2 oTexCoord;

      // Values that were scattered throughout the Oculus world demo
      const vec4 HmdWarpParam = vec4(1.0, 0.22, 0.24, 0.0); // For the 7-inch device
      const vec4 ChromAbParam = vec4(0.996, -0.004, 1.014, 0.0);
      const float HMD_HResolution = 1280.0;
      const float HMD_VResolution = 800.0;
      const float HMD_HScreenSize = 0.14976;
      const float HMD_VScreenSize = HMD_HScreenSize / (HMD_HResolution / HMD_VResolution);
      const float HMD_InterpupillaryDistance = 0.064;
      const float HMD_LensSeparationDistance = 0.0635;
      const float HMD_EyeToScreenDistance = 0.041;
      const float lensOffset = HMD_LensSeparationDistance * 0.5;
      const float lensShift = HMD_HScreenSize * 0.25 - lensOffset;
      const float Distortion_XCenterOffset = 4.0 * lensShift / HMD_HScreenSize;
      const float DistortionFitX = -1.0;
      const float DistortionFitY = 0.0;
      const float stereoAspect = 0.5 * HMD_HResolution / HMD_VResolution;
      const float dx = DistortionFitX - Distortion_XCenterOffset;
      const float dy = DistortionFitY / stereoAspect;
      const float fitRadiusSquared = dx * dx + dy * dy;
      const float Distortion_Scale =
        HmdWarpParam.x +
        HmdWarpParam.y * fitRadiusSquared +
        HmdWarpParam.z * fitRadiusSquared * fitRadiusSquared +
        HmdWarpParam.w * fitRadiusSquared * fitRadiusSquared * fitRadiusSquared;

      // This will be linked in later
      vec4 child_main(vec3 eye, vec3 ray);

      // This samples from a single unwarped [0,0]x[1,1] box containing two views
      // side-by-side that have been rendered using normal perspective projection.
      // The left eye takes up [0,0]x[0.5,1] and the right eye takes up [0.5,0]x[1,1].
      vec4 sample(vec2 point, vec3 eye, vec2 LensCenter) {
        float scale = 2.0;
        float dx = (point.x - LensCenter.x) * 2.0 * scale;
        float dy = (point.y - LensCenter.y) * scale;
        vec3 ray = axisX * dx + axisY * dy + axisZ;
        return child_main(eye, ray);
      }

      void main() {
        // Compute the viewport size
        bool isRight = oTexCoord.x > 0.5;
        float x = isRight ? 0.5 : 0.0;
        float y = 0.0;
        float w = 0.5;
        float h = 1.0;

        // Set up values for the shader
        float XCenterOffset = isRight ? -Distortion_XCenterOffset : Distortion_XCenterOffset;
        vec2 LensCenter = vec2(x + (w + XCenterOffset * 0.5) * 0.5, y + h * 0.5);
        vec2 ScreenCenter = vec2(x + w * 0.5, y + h * 0.5);
        float scaleFactor = 1.0 / Distortion_Scale;
        vec2 Scale = vec2(w * 0.5 * scaleFactor, h * 0.5 * scaleFactor * stereoAspect);
        vec2 ScaleIn = vec2(2.0 / w, 2.0 / h / stereoAspect);

        // Compute the warp
        vec2 theta = (oTexCoord - LensCenter) * ScaleIn;
        float rSq = dot(theta, theta);
        vec2 theta1 = theta * dot(HmdWarpParam, vec4(1.0, rSq, rSq * rSq, rSq * rSq * rSq));
        vec2 tc = LensCenter + theta1 * Scale;

        // Trace the warped ray
        vec3 localEye = eye + axisX * HMD_InterpupillaryDistance * 0.5 * (isRight ? 1.0 : -1.0);
        gl_FragColor = sample(tc, localEye, LensCenter);
      }

      // From: https://www.shadertoy.com/view/MdfGRr
      // Created by inigo quilez - iq/2013
      // License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.

      bool isphere(in vec4 sph, in vec3 ro, in vec3 rd, out vec2 t) {
        vec3 oc = ro - sph.xyz;
        float b = dot(oc, rd);
        float c = dot(oc, oc) - sph.w * sph.w;

        float h = b * b - c;
        if (h < 0.0)
          return false;

        float g = sqrt(h);
        t.x = -b - g;
        t.y = -b + g;

        return true;
      }

      const int NumIte = 7;
      const float Bailout = 1000.0;

      bool iterate(in vec3 p, in vec3 CC, out float resPot, out vec4 resColor) {
        vec3 zz = p;
        vec4 trap = vec4(abs(zz.xyz), dot(zz, zz));
        float dz = 1.0;

        for (int i = 0; i < NumIte; i++) {
          float m = dot(zz, zz);
          if (m > Bailout) {
            resColor = trap;
            resPot = 0.25 * log(m) * sqrt(m) / dz;
            return false;
          }

          dz = 8.0 * pow(m, 3.5) * dz;

          float x = zz.x;
          float x2 = x * x;
          float x4 = x2 * x2;
          float y = zz.y;
          float y2 = y * y;
          float y4 = y2 * y2;
          float z = zz.z;
          float z2 = z * z;
          float z4 = z2 * z2;

          float k3 = x2 + z2;
          float k2 = inversesqrt(k3 * k3 * k3 * k3 * k3 * k3 * k3);
          float k1 = x4 + y4 + z4 - 6.0 * y2 * z2 - 6.0 * x2 * y2 + 2.0 * z2 * x2;
          float k4 = x2 - y2 + z2;

          zz.x = CC.x + 64.0 * x * y * z * (x2 - z2) * k4 * (x4 - 6.0 * x2 * z2 + z4) * k1 * k2;
          zz.y = CC.y + -16.0 * y2 * k3 * k4 * k4 + k1 * k1;
          zz.z = CC.z + -8.0 * y * k4 * (x4 * x4 - 28.0 * x4 * x2 * z2 + 70.0 * x4 * z4 - 28.0 * x2 * z2 * z4 + z4 * z4) * k1 * k2;

          trap = min(trap, vec4(abs(zz.xyz), dot(zz, zz)));
        }

        resColor = trap;
        resPot = 0.0;
        return true;
      }

      bool ifractal(in vec3 ro, in vec3 rd, out float rest, in float maxt, out vec3 resnor, out vec4 rescol, float fov, vec3 ccc) {
        vec4 sph = vec4(0.0, 0.0, 0.0, 1.25);
        vec2 dis;

        if (!isphere(sph, ro, rd, dis))
          return false;

        // early skip
        if (dis.y < 0.001) return false;

        // clip to near!
        if (dis.x < 0.001) dis.x = 0.001;

        if (dis.y > maxt) dis.y = maxt;

        float dt;
        vec3 gra;
        vec4 color;

        float fovfactor = 1.0 / sqrt(1.0 + fov * fov);

        float t = dis.x;
        for (int i = 0; i < 80; i++) {
          vec3 p = ro + rd * t;

          float Surface = clamp(0.002 * t * fovfactor, 0.000001, 0.005);

          float eps = Surface * 0.1;
          vec4 col2;
          if (iterate(p, ccc, dt, color)) {
            rest = t;
            resnor = vec3(0.0, 0.0, 0.0);
            rescol = color;
            return true;
          }

          //gra = vec3(pot2-pot1, pot3-pot1, pot4-pot1);
          //dt = 0.01;
          if (dt < Surface) {
            rescol = color;

            vec4 tmp;
            float eps = Surface * 0.75;
            float p2;
            iterate(p + vec3(eps, 0.0, 0.0), ccc, p2, tmp);
            float p3;
            iterate(p + vec3(0.0, eps, 0.0), ccc, p3, tmp);
            float p4;
            iterate(p + vec3(0.0, 0.0, eps), ccc, p4, tmp);
            resnor = normalize(vec3(p2 - dt, p3 - dt, p4 - dt));
            rest = t;
            return true;
          }

          t += dt;
        }

        return false;
      }

      vec4 child_main(vec3 campos, vec3 rd) {
        float time = 3.5;

        vec3 light1 = vec3(0.577, 0.577, -0.577);
        vec3 light2 = vec3(-0.707, 0.000, 0.707);

        float fov = 1.5;

        vec3 cc = vec3(0.9 * cos(3.9 + 1.2 * time) - .3, 0.8 * cos(2.5 + 1.1 * time), 0.8 * cos(3.4 + 1.3 * time));
        if (length(cc) < 0.50) cc = 0.50 * normalize(cc);
        if (length(cc) > 0.95) cc = 0.95 * normalize(cc);

        vec3 nor;
        vec3 rgb;
        vec4 col;
        float t;

        if (!ifractal(campos, rd, t, 1e20, nor, col, fov, cc)) {
          rgb = 1.3 * vec3(1.0, .98, 0.9) * (0.7 + 0.3 * rd.y);

          rgb += vec3(0.8, 0.7, 0.5) * pow(clamp(dot(rd, light1), 0.0, 1.0), 32.0);
        } else {
          vec3 xyz = campos + t * rd;

          float dif1 = clamp(dot(light1, nor), 0.0, 1.0);
          float dif2 = clamp(0.5 + 0.5 * dot(light2, nor), 0.0, 1.0);
          float ao = clamp(1.5 * col.w - 0.9, 0.0, 1.0);
          float lt1;
          vec3 ln;
          vec4 lc;
          if (dif1 > 0.001)
            if (ifractal(xyz, light1, lt1, 1e20, ln, lc, fov, cc)) dif1 = 0.0;

          rgb = vec3(1.0, 1.0, 1.0) * 0.3;

          rgb = mix(rgb, vec3(1.0, 0.1, 0.0), sqrt(col.x));
          rgb = mix(rgb, vec3(1.0, 0.5, 0.2), sqrt(col.y));
          rgb = mix(rgb, vec3(1.0, 1.0, 1.0), col.z);

          vec3 brdf = 1.5 * vec3(0.17, 0.19, 0.20) * (0.6 + 0.4 * nor.y) * (0.1 + 0.9 * ao);
          brdf += 1.9 * vec3(1.00, 0.90, 0.60) * dif1 * (0.5 + 0.5 * ao);
          brdf += 1.1 * vec3(0.14, 0.14, 0.14) * dif2 * ao;

          rgb *= brdf;
        }

        rgb = sqrt(rgb);

        return vec4(rgb, 1.0);
      }
    )).link();

    quad << vec2(-1, -1) << vec2(1, -1) << vec2(-1, 1) << vec2(1, 1);
    quad.upload();
    layout.create(shader, quad).attribute<float>("vertex", 2).check();
}

static vec3 eye = vec3(0, 1, -5);
static vec3 axisX = vec3(1, 0, 0);
static vec3 axisY = vec3(0, 1, 0);
static vec3 axisZ = vec3(0, 0, 1);

static bool left = false;
static bool right = false;
static bool up = false;
static bool down = false;

static bool W = false;
static bool A = false;
static bool S = false;
static bool D = false;

static float angleX = 0;
static float angleY = 0;

static void draw() {
  glClear(GL_COLOR_BUFFER_BIT);
  shader.use();

  shader.uniform("eye", eye);
  shader.uniform("axisX", axisX);
  shader.uniform("axisY", axisY);
  shader.uniform("axisZ", axisZ);

  layout.draw(GL_TRIANGLE_STRIP);
  shader.unuse();
  glutSwapBuffers();
}

static void keyboard(unsigned char key, int x, int y) {
  switch (key) {
    case 27: exit(0); break; // Escape
    case 'w': case 'W': W = true; break;
    case 'a': case 'A': A = true; break;
    case 's': case 'S': S = true; break;
    case 'd': case 'D': D = true; break;
  }
}

static void keyboardUp(unsigned char key, int x, int y) {
  switch (key) {
    case 'w': case 'W': W = false; break;
    case 'a': case 'A': A = false; break;
    case 's': case 'S': S = false; break;
    case 'd': case 'D': D = false; break;
  }
}

static void specialKey(int key, int x, int y) {
  switch (key) {
    case 27: exit(0); break; // Escape
    case GLUT_KEY_LEFT: left = true; break;
    case GLUT_KEY_RIGHT: right = true; break;
    case GLUT_KEY_UP: up = true; break;
    case GLUT_KEY_DOWN: down = true; break;
  }
}

static void specialKeyUp(int key, int x, int y) {
  switch (key) {
    case GLUT_KEY_LEFT: left = false; break;
    case GLUT_KEY_RIGHT: right = false; break;
    case GLUT_KEY_UP: up = false; break;
    case GLUT_KEY_DOWN: down = false; break;
  }
}

static void centerCursor() {
  int screenWidth = glutGet(GLUT_SCREEN_WIDTH);
  int screenHeight = glutGet(GLUT_SCREEN_HEIGHT);
  glutWarpPointer(screenWidth / 2, screenHeight / 2);
}

static void tick() {
  float seconds = 0.01; // TODO

  if (left != right || up != down || W != S || A != D) {
    angleX += (S - W) * seconds;
    angleY += (D - A) * seconds;

    float sinX = sinf(angleX);
    float cosX = cosf(angleX);
    float sinY = sinf(angleY);
    float cosY = cosf(angleY);
    axisX = vec3(cosY, 0, -sinY);
    axisY = vec3(sinY * sinX, cosX, cosY * sinX);
    axisZ = vec3(sinY * cosX, -sinX, cosY * cosX);

    float deltaX = (right - left) * seconds;
    float deltaZ = (up - down) * seconds;
    float speed = 5;
    eye.x += (axisX.x * deltaX + axisZ.x * deltaZ) * speed;
    eye.y += (axisX.y * deltaX + axisZ.y * deltaZ) * speed;
    eye.z += (axisX.z * deltaX + axisZ.z * deltaZ) * speed;
  }

  glutPostRedisplay();
}

int main(int argc, char *argv[]) {
  printf("rotate: W/A/S/D\n");
  printf("move: up/down/left/right\n");
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
  glutCreateWindow("Example");
  glutDisplayFunc(draw);
  glutSpecialFunc(specialKey);
  glutSpecialUpFunc(specialKeyUp);
  glutKeyboardFunc(keyboard);
  glutKeyboardUpFunc(keyboardUp);
  glutIdleFunc(tick);
  glutFullScreen();
  setup();
  glutMainLoop();
  return 0;
}