ƒ : ℂ → ℂ


The source code for this site is maintained in git, and can be found on GitHub. Feel free to view the full code, including the website core, there.

Complex Plotter

This file implements the core of the plotter core. You'll note that the parser we use just inserts the parsed input straight into C++ code. Yes, that's massively insecure. We know. Your input is passed through a parser that pretties it up and also tries its best to generate safe code; see below for the parser.

#define OUTNAME "%(o)s"
#define WIDTH   %(w)s
#define HEIGHT  %(h)s
#define RANGE   {{%(l)s, %(b)s}, {%(r)s, %(t)s}}
#define F1      %(f)s

#include "EasyBMP/EasyBMP.h"
#include <complex>
#include "ccfunc.cc"

#define OUTPACK(r, g, b) (int(256*r) << 16) + (int(256*g) << 8) + (int(256*b));

inline int hsi2rgb(double H, double I) {
     double r, g, b;

     if(H==1) H = 0;
     double z = floor(H*6); int i = int(z);
     double f = H*6 - z;
     double p = I*(1-I);
     double q = I*(1-I*f);

     if (I > .999) I = .999;
     if (! (i & 1)) {
          q = I + p - q;
     }
 
     switch(i){
     case 0: r=I, g=q, b=p; break;
     case 1: r=q, g=I, b=p; break;
     case 2: r=p, g=I, b=q; break;
     case 3: r=p, g=q, b=I; break;
     case 4: r=q, g=p, b=I; break;
     case 5: r=I, g=p, b=q; break;
     }

     return OUTPACK(r, g, b);
}

inline double decc(cc a) {
     if (imag(a) == 0) {
          return real(a);
     } else {
          return abs(a);
     }
}

inline cc f1(cc z) {
     return F1;
}

BMP image;
double win[2][2] = RANGE;
double xstep = (win[1][0] - win[0][0]) / WIDTH;
double ystep = (win[1][1] - win[0][1]) / HEIGHT;

int run(int c, int r) {
     cc z = cc(win[0][0] + xstep * c, win[1][1] - ystep * r);
     cc w = log(-f1(z));

     double in = atan(real(w)) / pi + .5;
     double hue = imag(w) / (2 * pi) + .5;
     return hsi2rgb(hue, in);
}
	
int main() {
     image.SetSize(WIDTH, HEIGHT);
     cc z = cc(win[0][0], win[1][1]);
     for(int c=0; c < WIDTH; c++) {
          for(int r=0; r < HEIGHT; r++) {
            
               int color = run(c, r);
               RGBApixel* p = image(c, r);
               p->Red = (color >> 16) & 0xff;
               p->Green = (color >> 8) & 0xff;
               p->Blue = color & 0xff;
          }
     }
	
     image.WriteToFile(OUTNAME);
     return 0;
}

Complex Function Library

While the core of the complex functionality in the plotter is provided by the C++ STL, this library does not provide many functions mathematicians are accustomed to, either because these functions are simply shortcuts (tanh, for example), or because they are so rarely useful to programmers (e.g. zeta). Thus, these functions are implemented in a seperate library.

#include <complex>
using std::abs; using std::arg; using std::conj; using std::imag;
using std::real; using std::exp; using std::log; using std::pow;
using std::sqrt; using std::sin; using std::cos; using std::tan;
using std::sinh; using std::cosh;
#include <cstdlib>
using std::rand;

typedef std::complex<double> cc;

/* Constants */
const cc i(0, 1);
const double pi   = 3.141592653590;
const double e    = exp(1);
const double masc = .5772156649015;

/* Functions */

// Random
extern inline cc ccrand() {return cc(rand(), rand()) / (double) RAND_MAX;}

// Gamma
// See also: Numerical Recipies
double gamma_coeff[] = {75122.6331530, 80916.6278952, 36308.2951477, \
    8687.24529705, 1168.92649479, 83.8676043424, 2.50662827511};
extern inline cc gamma(cc z) {
    cc ret1(0, 0);
    cc ret2(1, 0);
    for (int i = 0; i < 6; i++) {
        ret1 += gamma_coeff[i] * pow(z, i);
        ret2 *= z + double(i);
    }
    return ret1 / ret2 * pow(z + 5.5, z + 0.5) * exp(-(z + 5.5));
}

extern inline cc beta(cc z, cc w) {
    return gamma(z)*gamma(w)/gamma(z + w);
}

// Dirichlet series

// See also: Borwein's method of calculating the eta function
int eta_coeff[] = {0, 200, 6800, 92180, 640400, 2690440, 7349640, \
    13903240, 19473800, 22095240, 22619520};
extern inline cc eta(cc s) {
    int n = 8;
    cc t(0, 0);
    for (int k = 0; k < n; k++) {
        t += -(pow(-1.0, k) * (eta_coeff[k] - eta_coeff[n]) / \
                pow(k + 1.0, s)) / double(eta_coeff[n]);
    }
    return t;
}
extern inline cc zeta(cc s) {return eta(s) / (1.0 - pow(2.0, 1.0-s));}
extern inline cc xi(cc s) {return (s-1.0)*gamma(s/2.0 + 1.0)*zeta(s)/sqrt(pow(pi, s));}

// Trig
extern inline cc cot(cc z)  {return 1.0 / tan(z);}
extern inline cc sec(cc z)  {return 1.0 / cos(z);}
extern inline cc csc(cc z)  {return 1.0 / sin(z);}

// HyperTrig
extern inline cc tanh(cc z) {return sinh(z) / cosh(z);}
extern inline cc coth(cc z) {return cosh(z) / sinh(z);}
extern inline cc sech(cc z) {return 1.0 / cosh(z);}
extern inline cc csch(cc z) {return 1.0 / sinh(z);}

// Arctrig
extern inline cc asin(cc z) {return .5*i * log(i*z + sqrt(1.0 - z*z));}
extern inline cc acos(cc z) {return .5*pi  + i * log(i*z + sqrt(1.0 - z*z));}
extern inline cc atan(cc z) {return .5*i * (log(1.0 - i*z) - log(1.0 + i*z));}
extern inline cc asec(cc z) {return .5*pi + i*(log(1.0 - 1.0/(z*z)) + i/z);}
extern inline cc acsc(cc z) {return -i * (log(1.0 - 1.0/(z*z)) + i/z);}
extern inline cc acot(cc z) {return .5*i * (log((z-i)/z) - log((z+i)/z));}

// ArcHyperTrig
extern inline cc asinh(cc z) {return log(z + sqrt(z*z + 1.0));}
extern inline cc acosh(cc z) {return log(z + sqrt(z*z - 1.0));}
extern inline cc atanh(cc z) {return .5*(log(1.0 + z) - log(1.0 - z));}
extern inline cc acsch(cc z) {return log(sqrt(1.0/(z*z) + 1.0) + 1.0/z);}
extern inline cc asech(cc z) {return log(sqrt(1.0/(z*z) - 1.0) + 1.0/z);}
extern inline cc acoth(cc z) {return .5*(log(1.0 + 1.0/z) - log(1.0 - 1.0/z));}

extern inline cc hyper2f1(cc a, cc b, cc c, cc z) {
    cc w(1, 0);
    cc q(1, 0);
    for (int i = 1; i < 25; i++) {
        double fi = (double)i;
        q *= (a + fi) * (b + fi) / (c + fi) / (1.0 + fi);
        w += q * pow(z, fi);
    }
    return w;
}

mathparse.parser

Your input on the web is not spliced directly into C++; instead, it is first passed through a parser that expands calls like z^5, and also works to ensure that only safe code is inserted into the image generator.

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import ply.yacc as yacc
from .lexer import tokens, literals

class YaccError(Exception): pass

start = "expr"
precedence = (
    ("left", "OR"),
    ("left", "AND"),
    ("right", "NOT"),
    ("left", "LE", "GE", "NE", "EQ", "<", ">"),
    ("left", "-", "+"),
    ("left",  "/", "*", "MUL"),
    ("right", "UMINUS"),
    ("right", "^"),
    ("left", "!"),
    ("left", "FNCALL"),
    ("left", "PARENTHESES"),
    ("left", "|"),
)

def p_kernel_numeric(p):
    """kernel : DEC
              | INT
              | BOOL"""
    p[0] = "cc(%.1f, 0.0)" % float(p[1])

def p_kernel_ident(p):
    """kernel : IDENT"""
    p[0] = p[1]

def p_expr_or(p):
    """expr : expr OR expr"""
    p[0] = "(" + p[1] + ") || (" + p[3] + ")"

def p_expr_and(p):
    """expr : expr AND expr"""
    p[0] = "(" + p[1] + ") && (" + p[3] + ")"

def p_expr_not(p):
    """expr : NOT expr"""
    p[0] = "!(" + p[2] + ")"

def p_expr_eq(p):
    """expr : expr NE expr
            | expr EQ expr"""
    p[0] = "double((" + p[1] + ") " + p[2] + " (" + p[3] + "))"

def p_expr_math(p):
    """expr : expr '+' expr
            | expr '-' expr
            | expr '*' expr
            | expr '/' expr"""
    p[0] = "(" + p[1] + ") " + p[2] + " (" + p[3] + ")"

def p_kernel_impmult(p):
    """kernel : expr expr %prec MUL"""
    p[0] = "(" + p[1] + ") * (" + p[2] + ")"

def p_expr_comp(p):
    """expr : expr LE expr
            | expr '<' expr
            | expr '>' expr
            | expr GE expr"""
    p[0] = "double(decc(" + p[1] + ") " + p[2] + " decc(" + p[3] + "))"

def p_expr_exp(p):
    """expr : expr '^' expr"""
    p[0] = "pow(" + p[1] + ", " + p[3] + ")"

def p_expr_postcall(p):
    """expr : expr '!' IDENT"""
    p[0] = p[3] + "(" + p[1] + ")"

def p_fncall_start(p):
    """fncall : IDENT '(' arg %prec FNCALL"""
    p[0] = p[1] + "(" + p[3]

def p_expr_emptyfncall(p):
    """expr : IDENT '(' ')' %prec FNCALL"""
    p[0] = p[1] + "()"

def p_arg_expr(p):
    """arg : expr"""
    p[0] = p[1]

def p_fncall(p):
    """fncall : fncall ',' arg"""
    p[0] = p[1] + ", " + p[3]

def p_expr_fncall(p):
    """expr : fncall ')'
            | fncall ',' ')'"""
    p[0] = p[1] + ")"

def p_expr_parenthesised(p):
    """expr : '(' expr ')' %prec PARENTHESES"""
    p[0] = p[2]

def p_expr_unary(p):
    """expr : '-' expr %prec UMINUS
            | '+' expr %prec UMINUS"""
    p[0] = p[1] + p[2]

def p_expr_conj(p):
    """expr : '~' expr %prec UMINUS"""
    p[0] = "conj(" + p[2] + ")"

def p_expr_abs(p):
    """expr : '|' expr '|'"""
    p[0] = "abs(" + p[2] + ")"

def p_expr_kernel(p):
    """expr : kernel"""
    p[0] = p[1]

def p_error(p):
    raise YaccError("Error parsing input. Try checking for syntax" \
            " errors (unbalanced parentheses, unfinished " \
            "expressions).")

if __name__ == "__main__":
    parser = yacc.yacc(debug=1)
else:
    parser = yacc.yacc(debug=0)

Licensing

All code on this page, and all other code in this web site not entered by users or with its own license header, is licensed under the terms of the GPL, version 3 or higher. If you are unsure of whether or not you may use this code in a project of your own, feel free to contact the authors.

All code is © 2009 Pavel Panchekha and Sam Fingeret