Complex Numbers

As mentioned elsewhere, I realized I wanted to use complex numbers in PB.
Complex numbers have the form x+iy where x is a real number, i is the square root of -1, and y is a real number multiplier of i.
More than this in explanation needed?
Here’s an intro

[long and live… you can also see umpteen other Youtube videos on complex numbers]

I found a javascript library for it, GitHub - infusion/Complex.js: A complex numbers library
but PB isn’t really JS, so lots of it needed to be ported and adjusted.

I’m not quite done (I didn’t bother porting the more obscure functions [yet]), and I’m pretty sure division by zero and other edge cases might not be right… but it works.

/**
 * Ported to PB's not quite JS by ScruffyNerf in 2021
 * based on Complex.js v2.0.13 12/05/2020 https://github.com/infusion/Complex.js/
 * used under the MIT license
 * 
 * This code allows the manipulation of complex numbers.
 * unlike the original code, you _must_ use a 2 element array
 *
 * The order of functions below is changed from the original, and I haven't ported all of the functions (yet?) 
 * I've moved needed core functions ahead of functions that depend on them.
 */

export var complex = array(2);
export var complex2 = array(2);
export var complex3 = array(2);
var re = 0;
var im = 1;
var rarray = array(2);
var varray = array(2);

var complexone = array(2);
complexone[re] = 1
complexone[im] = 0
var complexzero = array(2);
complexzero[re] = 0
complexzero[im] = 0
var complexi = array(2);
complexzero[re] = 0
complexzero[im] = 1
var complexpi = array(2);
complexzero[re] = PI
complexzero[im] = 0
var complexe = array(2);
complexzero[re] = E
complexzero[im] = 0

//missing from PB normal math functions

function cosh(v) {
    return (exp(v) + exp(-v)) / 2;
}

function sinh(v) {
    return (exp(v) - exp(-v)) / 2;
}

// hypot is in already PB v3  - not ported, v2 will require adding a shim for this

function logHypot(a, b) {
    return log(hypot(a, b));
}

/**
 * Calculate the magnitude (aka radius) of the complex number
 * @returns {number}
 */
function cabs(varray) {
    return hypot(varray[re], varray[im]);
}

/**
 * Calculate the angle of the complex number
 * @returns {number}
 */
function carg(varray) {
    return atan2(varray[im], varray[re]);
}

/**
 * Calculates the sign of a complex number, which is a normalized complex
 * @returns {Complex}
 */
function csign(varray) {
    var vabs = cabs(varray);
    rarray[re] = varray[re] / vabs;
    rarray[im] = varray[im] / vabs;
    return rarray;
}

/**
 * Adds two complex numbers
 * @returns {Complex}
 */
function cadd(a, b) {
	// infinity handling removed
    rarray[re] = a[re] + b[re];
    rarray[im] = a[im] + b[im];
    return rarray;
}

/**
 * Subtracts two complex numbers
 * @returns {Complex}
 */
function csub(a, b) {
    rarray[re] = a[re] - b[re];
    rarray[im] = a[im] - b[im];
    return rarray;
}

/**
 * Multiplies two complex numbers
 * @returns {Complex}
 */

function cmul(a, b) {
    // Short circuit for real values
    if (a[im] == 0 && b[im] == 0) {
        rarray[re] = a[re] * b[re];
        rarray[im] = 0;
    } else {
        rarray[re] = a[re] * b[re] - a[im] * b[im];
        rarray[im] = a[re] * b[im] + a[im] * b[re];
    }
    return rarray;
}

/**
 * Divides two complex numbers
 * @returns {Complex}
 */
function cdiv(a, b) {
    // needs some zero handling still... TBD
    var t, x;

    if (b[im] == 0) {
        // Divisor is real
        rarray[re] = a[re] / b[re];
        rarray[im] = a[im] / b[re];
        return rarray;
    }

    if (abs(b[re]) < abs(b[im])) {
        x = b[re] / b[im];
        t = b[re] * x + b[im];
        rarray[re] = (a[re] * x + a[im]) / t;
        rarray[im] = (a[im] * x - a[re]) / t;
        return rarray;
    } else {
        x = b[im] / b[re];
        t = b[im] * x + b[re];
        rarray[re] = (a[re] + a[im] * x) / t;
        rarray[im] = (a[im] - a[re] * x) / t;
        return rarray;
    }
}

/**
 * Calculate the power of two complex numbers
 * @returns {Complex}
 */
function cpow(a, b) {
    if (b[re] == 0 && b[im] == 0) {
        return complexone;
    }
	
    // If the exponent is real
    if (b[im] == 0) {
        if (a[im] == 0 && a[re] > 0) {
            rarray[re] = pow(a[re], b[re]);
            rarray[im] = 0;
            return rarray;
        } else if (a[re] == 0) { // If base is fully imaginary
            imgcase = ((b[re] % 4 + 4) % 4); // no switch case in PB, so we check each
            if (imgcase == 0) {
                rarray[re] = pow(a[im], b[re]);
                rarray[im] = 0;
                return rarray;
            }
            if (imgcase == 1) {
                rarray[re] = 0;
                rarray[im] = pow(a[im], b[re]);
                return rarray;
            }
            if (imgcase == 2) {
                rarray[re] = -pow(a[im], b[re]);
                rarray[im] = 0;
                return rarray;
            }
            if (imgcase == 3) {
                rarray[re] = 0;
                rarray[im] = -pow(a[im], b[re]);
                return rarray;
            }
        }
    }

    if (a[re] == 0 && a[im] == 0 && b[re] > 0 && b[im] >= 0) {
        return complexzero;
    }

    var arg = carg(a);
    var loh = logHypot(a[re], a[im]);
    v1 = exp(b[re] * loh - b[im] * arg);
    v2 = b[im] * loh + b[re] * arg;
    rarray[re] = v1 * cos(v2);
    rarray[im] = v1 * sin(v2);
    return rarray;
}

/**
 * Calculate the complex square root
 * @returns {Complex}
 */
function csqrt(varray) {

    var a = varray[re];
    var b = varray[im];
    var r = cabs(varray);
    var sre, sim;

    if (a >= 0) {
        if (b == 0) {
            rarray[re] = sqrt(a);
            rarray[im] = 0;
            return rarray;
        }
        sre = 0.5 * sqrt(2 * (r + a));
    } else {
        sre = abs(b) / sqrt(2 * (r - a));
    }

    if (a <= 0) {
        sim = 0.5 * sqrt(2 * (r - a));
    } else {
        sim = abs(b) / sqrt(2 * (r + a));
    }

    rarray[re] = sre;
    rarray[im] = b < 0 ? -sim : sim;
    return rarray;
}

/**
 * Calculate the complex exponent
 * @returns {Complex}
 */
function cexp(varray) {
    var tmp = exp(varray[re]);
    rarray[re] = tmp * cos(varray[im]);
    rarray[im] = tmp * sin(varray[im]);
    return rarray;
}

/**
 * Calculate the natural log
 * @returns {Complex}
 */
function clog(varray) {
    rarray[re] = logHypot(varray[re],varray[im]);
    rarray[im] = carg(varray);
    return rarray;
}

/**
 * Calculate the sine of the complex number
 * @returns {Complex}
 */
function csin(varray) {
    rarray[re] = sin(varray[re]) * cosh(varray[im]);
    rarray[im] = cos(varray[re]) * sinh(varray[im]);
    return rarray;
}

/**
 * Calculate the cosine
 * @returns {Complex}
 */
function ccos(varray) {
    rarray[re] = cos(varray[re]) * cosh(varray[im]);
    rarray[im] = -sin(varray[re]) * sinh(varray[im]);
    return rarray;
}

/**
 * Calculate the tangent
 * @returns {Complex}
 */
function ctan(varray) {
    var tre = 2 * varray[re];
    var tim = 2 * varray[im];
    var td = cos(tre) + cosh(tim);
    rarray[re] = sin(tre) / td;
    rarray[im] = sinh(tim) / td;
    return rarray;
}

// more to come...

// some quick demo code, with sliders, to adjust.  Very flowy
export var time1 = 1
export function sliderTime1(v) {
    time1 = v * 5
}
export var time2 = 1
export function sliderTime2(v) {
    time2 = v * 5
}
export var time3 = 1
export function sliderTime3(v) {
    time3 = v * 5
}

export var t1,t2,t3
export function beforeRender(delta) {
    t1 = wave(time(time1)) * 2 - 1
    t2 = wave(time(time2)) * 2 - 1
    t3 = wave(time(time3)) * 2 - 1
}
export var h
export function render2D(index, x, y) {
    complex[re] = (x - .5 + t1)
    complex[im] = (y - .5 + t2)
    complex2[re] = t1
    complex2[im] = t2
    complex3 = cmul(complex, complex2)
    h = cabs(complex3) * t3
    s = 1
    v = 1
    hsv(h, s, v)
}
3 Likes

some examples of the sort of graphs possible:

interactive:
https://talbrenev.com/complexgrapher/

(^ I’m going to port this functionally to PB, ala tixy, it’s really fast and looks great )

interactive:
https://jutanium.github.io/ComplexNumberGrapher

http://fs2.american.edu/lcrone/www/ComplexPlot.html

https://vqm.uni-graz.at/pages/complex/index1.html

http://www.nucalc.com/ComplexFunctions.html

https://www.jedsoft.org/fun/complex

I keep adding to this list, as I find neat cool examples and explainations

https://www.dynamicmath.xyz/domain-coloring/

and finally (the only one that seems to preview)

1 Like

@scruffynerf, I can put this librarary right to work! Thank you!

1 Like

Hey, @wizard , technical questions related to improving this library…

How does PB handle division by zero?
Do we have NaN or Infinity, like JavaScript?(guessing no)

Since all numbers are 16.16, trying to figure out the best to handle those cases. Since these are artificial constructs as a 2 element array, and we don’t have objects or strings, trying to figure out the best way (if possible) to define NaN and/or Infinity, or I might just have to leave them out.

How do you handle it for tan() and other functions where it comes up?

Awesome! I can use this, too!

@wizard native support would offer more reasons for me to need a heat sink on my PB!

1 Like

Thought I’d run into this before, so I just retested – dividing anything, including zero by zero gives you zero.

Though not mathematically correct, it is probably the least painful behavior. I’d document it and keep it this way for the time being.

1 Like

Division by zero in Pixelblaze results in a zero, as do other undefined results. In JavaScript you get Infinity which is a special value in IEEE 754 floating point. For fixed point math, there’s no where to hide extra bits of data to indicate Infinity or NaN.

A floating point Pixelblaze implementation might happen, its on my wishlist, in which case I’d emulate JavaScript.

2 Likes