Processing Complex Numbers

Processing Complex Numbers

A complex number has a real part and an imaginary part — 3 + 4i — and every REAL in Sheerpower can hold one. There is no complex type to declare and no library to load: write the imaginary unit as a suffix on a number, and the arithmetic, the text forms and the functions follow. This page covers the values themselves (Part 1), the precision and the magnitude range they carry (Part 2), three worked stories (Part 3), complex arrays (Part 4), the transforms with their stories (Part 5), what the operations cost (Part 6), and the functions of a complex number (Part 7). The operators page (Mathematical and Logical Operators) introduced the basics; the array pages (Array Math Functions, Index Lists & Slices, Solve(), Sort(), and More) matter for the transforms that come later on this page.

Part 1: Writing and Reading Complex Numbers

4i is a literal, and 3 + 4i is a real plus an imaginary, which the compiler folds into one complex value. From computed parts use complex(re, im); from text, val(). The four arithmetic operators work as you expect, and a real next to a complex is promoted. The functions that take a complex apart: real(), imag(), abs() (the magnitude), conj() (the sign of the imaginary part flipped) and arg() (the angle in radians, from -pi to pi).

Two rules keep this simple. A value that is complex stays complex: z - 4i is 3+0i, not 3, so a program can always tell what it is holding; real(z) is the way back to an ordinary number. And a complex number has no order: z < w raises the catchable exception NOORDER, and so do max, min and clamp; compare abs(), real() or imag(), which is what you meant anyway.

Part 2: Precision and Magnitude

Each part of a complex number is a decimal floating-point value with 30 significant digits — twice a double's 15 — and its own exponent from 10^-8191 to 10⁸¹⁹¹. Every + - * / is computed exactly and rounded once, to the nearest 30-digit value (ties to even), so the result of one operation is always the correctly rounded answer, never the sum of several roundings. abs() is a correctly rounded square root of the exact sum of squares. Printing shows 16 digits per part by default, and sprintf$ shows as many as you ask for, up to 30.

The limits are loud, not silent. A part whose exponent would pass 8191 raises NUMOVER (the same exception an INTEGER overflow raises), and there is no complex NaN or infinity to absorb a mistake: a NaN operand gives NaN, an infinite one raises. Signed zeros are kept and shown, because the branch cuts of the transcendental functions depend on them: -complex(0, 0) prints -0-0i and conj(-4 + 0i) prints -4-0i, whose arg() is -pi where arg(-4 + 0i) is pi.

Part 3: Three Stories

An AC circuit. Impedance is a complex number: the resistance is its real part, the reactance (inductive minus capacitive) its imaginary part. Ohm's law then reads exactly as it does for direct current, current = volts / impedance, and the phase angle between voltage and current is arg(current). The apparent power is volts * conj(current): its real part is the real power in watts, its imaginary part the reactive power, its magnitude the volt-amperes.

A rotation in the plane. Multiplying by a complex number of magnitude 1 rotates a point about the origin, which is how a 2-D graphics or geometry program turns things without a matrix. 1i is a quarter turn, exactly; for any other angle build the unit number from cos and sin.

The Mandelbrot set. The classic loop: square a complex number and add the starting point, over and over, and count how many steps it takes to escape past magnitude 2. A point that never escapes belongs to the set. One routine, one abs() per step, and a line of the picture falls out of a for loop.

Part 4: Complex Arrays

An array of REALs holds complex elements the way a variable does, and the array pages' machinery carries over: element-wise arithmetic, the maps, masks, filter(), where(), the sorts and the linear algebra. abs(), real() and imag() turn a complex array into a real one, which is how the numeric statistics and comparisons get at it.

Two things to know. sort(z), unique() and isin() order complex values by real part and then imaginary part — a bookkeeping order that makes them work, not a size: to sort by magnitude, sort by abs(z) and gather (z[sortindex(abs(z))]). And the statistics that need an order — stats$max, min, median, percentile, argmax and their kin — raise NOORDER on a complex element, while the ones that need a square (stats$var, stddev, pcorr, ...) raise NOCOMPLEX; stats$sum and stats$mean are defined. Matrices with complex coefficients go through matmul(), solve() and determinant() with the same exact-then-round arithmetic, and cluster input reads a cell written as 3+4i into a REAL field as a complex value.

Part 5: Transforms — fft(), ifft() and fftfreq()

The discrete Fourier transform turns a signal sampled in time into the tones it is made of. fft(x) takes an array of numbers (or complex numbers) and gives back an array of the same length whose elements are complex: element k is the amplitude and phase of the tone that repeats k times over the whole sample. ifft(z) goes back. fftfreq(n, spacing) gives the frequency that goes with each element, so you can read a spectrum in hertz instead of bin numbers. Any length works: a power of two is fastest, anything else is handled too.

A vibration sensor. A thousand samples a second for one second from a machine: a 50 Hz hum from the mains, a weaker 120 Hz tone from a bearing, and noise. The power spectrum is abs(fft(x)) ^ 2, the frequency axis is fftfreq(n, 1 / rate), and stats$argmax over the first half (the second half mirrors it) names the strongest tone.

Touch tones. A telephone key is two tones at once: the 1 key is 697 Hz and 1209 Hz. Sampled at 8 kHz for a tenth of a second (800 samples), the two strongest bins are 700 and 1210 Hz — not 697 and 1209, because 800 samples at 8 kHz resolve frequencies to the nearest 10 Hz (rate / n). A longer sample resolves finer; that trade between length and resolution is the whole craft of reading a spectrum.

The weekly cycle in daily sales. Twelve weeks of daily sales have a level, a weekend bump that repeats every seven days, and a slow trend. The transform finds the seven-day cycle at once (the strongest bin, DC aside, sits at 1/7 cycles per day), and zeroing that bin and its mirror before transforming back removes the cycle and leaves the level and the trend — a seasonal adjustment in four statements.

Smoothing by convolution. A moving average is a convolution, and a convolution is a multiplication of spectra: transform both, multiply, transform back. (This is the circular convolution — the ends wrap — which is what the kernel of zeros padding is for.)

What the transform costs, this machine (the 2026-09-09 build):

The refusals are loud: a NaN or infinity in the input raises NUM_OUTOFRANGE naming the element (fill holes first with where(isnan(x), 0, x)), a value beyond a double's range raises NUMOVER and names digits: 30, which takes it, and a string array is a compile error.

Part 6: What the Operations Cost

A complex operation is exact-then-round-once at 30 digits, and it is still fast: the kernel does its products in 64-bit pieces and rounds with a reciprocal multiply, no division instruction anywhere. Measured in an ordinary for loop (200,000 iterations, one statement per iteration, including the statement's own overhead:

So a complex multiply or add costs about two real multiplies. Division is the exception at under a million a second — it runs the general exact path — so in a hot loop divide once and multiply by the result where the algebra allows. A million complex multiply-adds is about a third of a second. For scale, a radix-2 FFT written in Sheerpower itself runs about 700,000 butterflies a second (65,536 points in 0.75 s); the built-in transforms of the parts that follow do the same arithmetic at about 3.3 million butterflies a second, a million points in about three seconds.

Part 7: The Functions of a Complex Number

sqr(), exp(), log() (and log10, log2), the trig family (sin, cos, tan, csc, sec, cot, asin, acos, atn), the hyperbolic one (sinh, cosh, tanh) and a non-whole power all take a complex number. They are computed in doubles, like their REAL versions, with the standard branch cuts of the C language (Annex G): the principal square root, the logarithm's cut along the negative real axis, and the signed zero deciding which side of a cut you are on. A whole-number power stays exact.

The last line shows the precision: these are double-precision results (15 to 16 digits), so an identity comes back to that many digits and no more — exactly what exp(log(x)) does for a REAL x, whose exp and log also run in doubles. Two rules follow. A result that is not finite raises NUM_OUTOFRANGE naming the function (log(0 + 0i), csc(0 + 0i)), and a zero base to a negative or complex power raises EXPZERONEG. And a real argument keeps its real rules: sqr(-4) still raises, because promoting it silently would turn a typo into an imaginary number; write sqr(-4 + 0i) when you mean the complex root. The functions map over arrays like every other built-in.