The last time Hackerfall tried to access this page, it returned a not found error. A cached version of the page is below, or clickhereto continue anyway

Complex numbers in JavaScript / Observable

Complex numbers in JavaScript

Our goal for today is to perform the dot product of two complex vectors in JavaScript. To make it a bit more challenging we are going to make use of typed arrays for efficient data storage.

JavaScript typed arrays are array-like objects and provide a mechanism for accessing raw binary data. As you may already know, Array objects grow and shrink dynamically and can have any JavaScript value. JavaScript engines perform optimizations so that these arrays are fast. However, as web applications become more and more powerful, adding features such as audio and video manipulation, access to raw data using WebSockets, and so forth, it has become clear that there are times when it would be helpful for JavaScript code to be able to quickly and easily manipulate raw binary data in typed arrays.

Mozilla Developer Network

The basis of our implementation will be the zdotc routine from the Reference BLAS package by Netlib. We'll go through that a bit later, but let's start with a crash course in complex arithmetic.

A complex number is a number that can be expressed in the form , where and are real numbers, and is a solution of the equation .

Wikipedia

So, a complex number consists of two parts, one real part and one imaginary part .

We can also imagine it as a coordinate in a 2 dimensional space where is the x axis and is the y axis.

Enough chit chat, let's write it up!

class Complex16 extends Float64Array {
  constructor(...args) {
    super(...args);
    
    if (this.length !== 2) {
      throw new Error('invalid complex number');
    }
  }
  
  // setters
  set re(value) { this[0] = value; }
  set im(value) { this[1] = value; }
  
  // getters
  get re() { return this[0]; }
  get im() { return this[1]; }
  
  // for visualization purposes only
  toString() {
    const { im, re } = this;

    let real = '';
    if (re !== 0 || im === 0) {
      real = re;
    }
    
    let op = '';
    if (re !== 0 && im > 0) {
      op = '+';
    }
    
    let imaginary = '';
    if (im !== 0) {
      if (im === 1) {
        imaginary = `i`;
      } else if (im === -1) {
        imaginary = `-i`
      } else {
        imaginary = `${im}i`;
      }
    }

    return `${real}${op}${imaginary}`;
  }
}

For readability purposes I'm subclassing Float64Array, which let's us borrow some of its internal methods. We are ready to create our first complex number!

Hint: Try changing a number and press on your keyboard.

{
  const z = new Complex16(2);
  
  z.re = 4;
  z.im = 3;

  return tex.block`
    \begin{aligned}
    z &= ${z} \\
    \operatorname{Re}(z) &= ${z.re} \\
    \operatorname{Im}(z) &= ${z.im} \\
    \end{aligned}
  `;
}

Addition

To reach our final goal, we need to build a toolbox of functions that operate on complex numbers.

Our first member is complex addition and it is performed by adding the real and imaginary parts of two numbers separately.

function add(z1, z2) {
  return new Complex16([
    z1.re + z2.re,
    z1.im + z2.im
  ]);
}
{
  const z1 = new Complex16([2, 1]);
  const z2 = new Complex16([3, 1]);

  return tex.block`
    \begin{aligned}
    z_1 &= ${z1} \\
    z_2 &= ${z2} \\
    z_1+z_2 &= ${add(z1, z2)}
    \end{aligned}
  `;
}

Multiplication

The second tool we are going to need is complex multiplication. What we need to keep in mind here is that the square of equals . Otherwise it's just normal multiplication.

function mul(z1, z2) {
  return new Complex16([
    z1.re * z2.re - z1.im * z2.im,
    z1.re * z2.im + z1.im * z2.re
  ]);
}
{
  const z1 = new Complex16([2, -1]);
  const z2 = new Complex16([3, 1]);

  return tex.block`
    \begin{aligned}
    z_1 &= ${z1} \\
    z_2 &= ${z2} \\
    z_1z_2 &= (${new Complex16([z1.re * z2.re, 0])}) + (${new Complex16([0, z1.re * z2.im])}) + (${new Complex16([0, z1.im * z2.re])}) + (${new Complex16([0, z1.im * z2.im])}^2) \\
    &= ${mul(z1, z2)} \\
    \end{aligned}
  `;
}

Conjugation

Now we just need one more piece of arithmetic, the complex conjugate.

In mathematics, the complex conjugate of a complex number is the number with an equal real part and an imaginary part equal in magnitude but opposite in sign.

Wikipedia

Well, that's not even a challenge. We just have to negate the imaginary part.

function conjg(z) {
  return new Complex16([z.re, -z.im]);
}
{
  const z = new Complex16([1, 2]);

  return tex.block`
    \begin{aligned}
    z &= ${z} \\
    \overline{z} &= ${conjg(z)}
    \end{aligned}
  `;
}

Up until now we have been working with a single complex value, Complex16, for our calculations. To be able to compute the dot product of two complex vectors, we need to introduce a complex vector data structure.

If we were lazy we could just save a bunch of Complex16 instances in a JavaScript Array and loop through it. But we like it fast, right? Regular JavaScript arrays are not.

It turns out that an efficient way of packing complex vectors is to store one complex number after the other, like this:

That is good news for us, it means we can still use Float64Array as our base class, just in a slightly different way.

class Complex16Array extends Float64Array {
  constructor(...args) {
    super(...args);

    if (this.length % 2 !== 0) {
      throw new Error('you supplied an uneven amount of elements'); 
    }
  }

  // get the ith complex number from the vector
  get(i) {
    const pos = i * 2;
    this.checkBounds(pos);
    // see https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/TypedArray/subarray
    return new Complex16(this.subarray(pos, pos + 2));
  }
  
  // set the ith complex number in the vector
  set(i, z) {
    const pos = i * 2;
    this.checkBounds(pos);
    this[pos] = z.re;
    this[pos + 1] = z.im;
  }
  
  // helper to check if accessed index is out of bounds
  checkBounds(pos) {
    if (pos < 0 || pos > this.length - 2) {
      throw new Error('index out of bounds');
    }
  }
  
  // helper to get count of complex numbers instead of actual length of the array
  get size() {
    return this.length / 2;
  }
  
  // again, just for visualization purposes...
  toString() {
    const res = [];
    for (let i = 0; i < this.size; i++) {
      res.push(this.get(i));
    }
    
    return `(${res.join(', ')})`;
  }
}

Here is our complex vector in action:

{
  // allocate new complex vector with 4 memory slots
  const u = new Complex16Array(4);
 
  u.set(0, new Complex16([1, 2]));
  u.set(1, new Complex16([3, 4]));
  
  const out = [];
  for (let i = 0; i < u.size; i++) {
    out.push(`u_${i} &= ${u.get(i)} \\\\`);
  }
  
  return tex.block`
    \begin{aligned}
    \vec u &= ${u} \\
    ${out.join('')}
    \end{aligned}
  `;
}

Dot product

The complex dot product is (usually) defined as the sum of the th conjugated element of multiplied by the th element of .

I mentioned the zdotc routine in the beginning of the article, so we're going to take a look at it now. For those of you that have never heard of BLAS, it stands for Basic Linear Algebra Subprograms and they are the de facto standard low level routines for linear algebra libraries. They are written in Fortran.

To make it a bit shorter I took the liberty of ignoring increments (incx and incy) and removing comments. You can see the full source code here.

       COMPLEX*16 FUNCTION zdotc(N,ZX,ZY)
 *
 *     .. Scalar Arguments ..
       INTEGER N
 *     ..
 *     .. Array Arguments ..
       COMPLEX*16 ZX(*),ZY(*)
 *
 *     .. Local Scalars ..
       COMPLEX*16 ZTEMP
       INTEGER I
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC dconjg
 *     ..
       ztemp = (0.0d0,0.0d0)
       zdotc = (0.0d0,0.0d0)
       IF (n.LE.0) RETURN

       DO i = 1,n
          ztemp = ztemp + dconjg(zx(i))*zy(i)
       END DO

       zdotc = ztemp
       RETURN
       END

Not too hard to translate!

function zdotc(n, zx, zy) {
  let ztemp = new Complex16([0, 0]);
  if (n < 0) {
    return ztemp;
  }

  for (let i = 0; i < n; i++) {
    ztemp = add(ztemp, mul(conjg(zx.get(i)), zy.get(i)));
  }

  return ztemp;
}
{
  const n = 3;
  const u = new Complex16Array([1, 2, 3, 4, 5, 6]);
  const v = new Complex16Array([4, 3, 2, 1, 0, -1]);
  
  return tex.block`
    \begin{aligned}
    \vec{u} &= ${u} \\
    \vec{v} &= ${v} \\
    \vec{u}\cdot\vec{v} &= \sum_{i = 0}^n{\overline{u_i} v_i} \\
    &= ${zdotc(n, u, v)}
    \end{aligned}
  `;
}

If you have feedback, questions, or maybe even suggestions for a future post you'll find my contact info in my GitHub profile.

Continue reading on beta.observablehq.com