Source file: complex-numbers.fut

Complex numbers

Futhark does not have complex numbers, but they are fairly easy to implement as a module. First we define a module type describing the interface of complex numbers:

module type complex = {

The type of the components of the complex number.

  type real

The type of complex numbers.

  type complex

Construct a complex number from real and imaginary components.

  val mk : real -> real -> complex

Construct a complex number from just the real component. The

imaginary part will be zero.

  val mk_re : real -> complex

Construct a complex number from just the imaginary component. The

real part will be zero.

  val mk_im : real -> complex

Conjugate a complex number.

  val conj : complex -> complex

The real part of a complex number.

  val re : complex -> real

The imaginary part of a complex number.

  val im : complex -> real

The magnitude (or modulus, or absolute value) of a complex number.

  val mag : complex -> real

The argument (or phase) of a complex number.

  val arg : complex -> real

  val (+) : complex -> complex -> complex
  val (-) : complex -> complex -> complex
  val (*) : complex -> complex -> complex
  val (/) : complex -> complex -> complex

  val sqrt : complex -> complex
  val exp : complex -> complex
  val log : complex -> complex
}

Then we define a parametric module that given a real module (in practice f32 or f64) will give us back a complex number module. We use local opens to perform operations using the operators from the T module. For examle, T.(x + y) means the same as x T.+ y, but is more readable. The definitions here are straight out of a textbook, so we won’t elaborate them.

module mk_complex (T: real)
  : (complex
     with real = T.t
     with complex = (T.t, T.t)) = {
  type real = T.t
  type complex = (T.t, T.t)

  def mk (a: real) (b: real) = (a, b)
  def mk_re (a: real) = (a, T.i32 0)
  def mk_im (b: real) = (T.i32 0, b)

  def conj ((a, b): complex) = T.((a, i32 0 - b))
  def re ((a, _b): complex) = a
  def im ((_a, b): complex) = b

  def (+) ((a, b): complex) ((c, d): complex) = T.((a + c, b + d))
  def (-) ((a, b): complex) ((c, d): complex) = T.((a - c, b - d))

  def (*) ((a, b): complex) ((c, d): complex) =
    T.(( a * c - b * d
       , b * c + a * d
       ))

  def (/) ((a, b): complex) ((c, d): complex) =
    T.(( (a * c + b * d) / (c * c + d * d)
       , (b * c - a * d) / (c * c + d * d)
       ))

  def mag ((a, b): complex) =
    T.(sqrt (a * a + b * b))

  def arg ((a, b): complex) =
    T.atan2 b a

  def sqrt ((a, b): complex) =
    let gamma = T.(sqrt ((a + sqrt (a * a + b * b)) / i32 2))
    let delta =
      T.(sgn b
         * sqrt (((i32 0 - a) + sqrt (a * a + b * b)) / i32 2))
    in (gamma, delta)

  def exp ((a, b): complex) =
    let expx = T.exp a
    in mk T.(expx * cos b) T.(expx * sin b)

  def log (z: complex) =
    mk (T.log (mag z)) (arg z)
}

We instantiate the complex number module using double-precision floats for the components:

module c64 = mk_complex f64

And try it out in the REPL:

> c64.exp (c64.mk 1 2)
(-1.1312043837568135f64, 2.4717266720048188f64)

The github.com/diku-dk/complex package (which you should use in real programs rather than rolling your own) is implemented in exactly this way.

Complex numbers

See also