Complex numbers

[Snektron]

This is a concrete proposal for adding complex numbers to Zig. The basic idea has already been greenlit informally, so this issue is here to describe the concrete roadmap that is required for implementing these.

The main idea is to add complex numbers as a first-class citizen instead of as library types. There are 2 main reasons for this:

• C interop: C99 introduces complex numbers (See C99 spec annex G) as a language built-in type. Unfortunately, their ABI is unspecified, and in practice different things are done depending on the target. This means that currently, it is almost impossible to interface with libraries that expose complex numbers in their ABIs. The easiest way to fix this up is within the compiler itself, and switch on the target when interfacing with those.
• Ease of use. Math operators are relatively annoying to use as free-standing functions, and being able to use the regular operators with complex numbers helps in these case. While this feature can be excused to keep the language grammar and semantics simpler, it seems appropriate to have it if there is a built-in complex number type anyway. Slightly related: Proposal: Short Vector Primitives #7295

This proposal seeks to deal with both of the above, though the primary argument for having built-in complex number is the former.

Scope

This proposal seeks to implement the following things:

• First-class complex numbers in the Zig language. We are specifically interested in the floating point versions, which are equivalent to _Complex float, _Complex double, and _Complex long double. We can extend the C family with to the usual Zig types; complex versions of f16, f32, f64, f80, f128, comptime_float.
• Supplementary built-in functions where required.
• Layout is fixed: It needs to be C-compatible, so the same memory layout as [2]T, where T is the backing float type. Layout is real first, then imaginary.
• Standard library support. This means implementing support for complex numbers in std.math where required. In time, std.math.Complex would be deprecated, and uses are replaced by the new type.
• compiler-rt support. This is something that is already required, and partly implemented using std.math.Complex. Note that the implementations here in particular suffer from ABI issues: the complex type is here passed as a struct but as stated before, this is not always correct.

Out of scope

• Integer complex numbers. These are not in standard C99, but are supported by some major compilers (GCC, Clang). By making the complex number type be generic over the backing type, we can leave the door open for an implementation of these in the future.
• An equivalent for the _Imaginary (see C99 spec) type specifier. This is an optional component of C99, and is in practice not really used [citation needed]. Besides, these types have the same ABI requirements as their backing type, and so these could be dealt with in the C translator without requiring additional changes in Zig.
• Quaternions, small vectors, etc.
• Complex numbers with representations that do not match C.

Language Changes

In order to support complex types, the following changes would need to be made:

• A new complex number type which is generic over the backing type that represents both the real and imaginary parts of the complex value. This variant would be added to builtin.Type.
• A type constructor complex types, which is generic over the backing type. This could be something like @Complex(T), that creates a complex number type where both fields are backed by T.
• A way to access the fields of a complex value. I have 2 main methods in mind:
  • Via a built-in function such as @real and @imag (similar to creal and cimag. This would also support the following syntax:

    var z: @Complex(f32) = ...;
    // Fetch real part
    const real = @real(z);
    // Fetch imaginary part
    const imag = @imag(z);
    // Assign real
    @real(z) = 123;
    // Assign imaginary
    @imag(z) = 456;

  • Via fields. This would make the complex type be more like a struct:

    var z: @Complex(f32) = ...;
    const real = z.real;
    const imag = z.imag;
    z.real = 123;
    z.imag = 456;

• A way to construct the complex numbers. I have considered a few methods, but im not 100% convinced of either of these:
  • A built-in function @I(), that returns 0+1i. This could be used to create more complex complex types, for example 1 + 1 * @I(). It looks a bit weird imo.
  • A type suffix, for example im: 1+1im would return the equivalent complex value. Since Zig does not have any other numeric suffixes, I think this option is not the best.
  • A built-in function to construct a complex value. For example: @complex(1, 1). The result would have the same type as the peer type of the arguments. For comptime values, this would be @Complex(comptime_float).
  • Support tuple initialization syntax: var z: @Complex(f32) = .{ 1, 1 };.
• Coercion rules between complex types. These come in two main forms:
  • Complex(T) -> Complex(U). This has the same rules as T -> U.
  • T -> Complex(T). This would initialize a complex value with the value x + 0i. The main motivation for this feature is to simplify syntax when dealing with real and complex values, and allow the user to write x + z. These operations are well-defined in mathematics.
• Peer-type resolution: @TypeOf(@Complex(T), U) would yield @Complex(@TypeOf(T, U)).
• Overloads for operators that apply on floats:
  • Binary arithmetic +, -, *, /.
    • Multiplication and division have the usual complex implementations.
    • Is % necessary? C does not define this one.
  • Assignment variants of the above.
  • Binary comparative ==, !=.
    • Only if both components are pairwise equal, a complex number compares equal to another.
  • Unary -.
    • This inverts both the real and imaginary parts.
  • Relational operators, <=, <, >, >=?
    • Relational operators are not supported in C. This usually compares the modulus, and would require the implementation of abs to be supported in the compiler. It can be implemented by comparing the hypot without computing the square root.
• Support for complex number in existing built-in functions:

  • @mulAdd.

  • @floatCast.

  • @max, @min.

  • @sqrt.

  • @sin, @cos, @tan.

  • @exp, @exp2.

  • @log, @log2, @log10.

  • @fabs. This would implement the hypot. GingerBill shared some notes on this:

    expand

    I've just had to do all of this today in fact: https://github.com/odin-lang/Odin/blob/master/core/math/cmplx/cmplx_trig.odin
    Odin/core/math/cmplx/cmplx_trig.odin at master · odin-lang/Odin
    Odin Programming Language. Contribute to odin-lang/Odin development by creating an account on GitHub.
    Odin/core/math/cmplx/cmplx_trig.odin at master · odin-lang/Odin
    A really simple example of complex numbers not working as you expect would be abs(c)
    You'd trivially thing it would be as simple as this:

    p, q := real(z), imag(z)
    return sqrt(pp + qq)

    But if any has ever implemented the hypot algorithm, you'll know this isn't accurate.
    A more accurate approach would be to do this:

    p, q := abs(real(z)), abs(imag(z))
    if p < q {
    p, q = q, p
    }
    if p == 0 {
    return 0
    }
    q = q / p
    return p * sqrt(1 + q*q)

    This might seem a little weird and also computationally more expensive too, but you are doing four different tricks here.
    First is not relying on pp to be the absolute square component in places.
    Second is ordering the values so that p is always the largest out of the two
    Third if p == 0, then so is q, and thus the entire result is 0
    Fourth is to change sqrt(pp + qq) into p * sqrt(1 + (qq)/(p*p))
    This makes the sqrt a lot more accurate because it will make the sqrt in the range of [1, 2]
    All of these little things that you may not have realized when coming to complex numbers and floating point.
    (and this is the simplest example too, and simplified too)

Open Questions

• Should vectors of complex numbers be supported? How would those be represented?
• Should we support relational operators?
  • Resolved: No, the semantics of this operation is unclear and not well defined.
• Should we support %?
  • Resolved: No, see previous.
• ...

Related Issues

• Proposal: make std.math.Complex ABI-compatible with C _Complex #8465
• std.math.complex: c_Complex for C maths libraries #11961

Please bikeshed on syntax below

[rofrol]

why bikeshed? I am not native speaker but it seems strange to ask for bikeshed:

Futile expenditure of time and energy in discussion of marginal technical issues

https://en.m.wiktionary.org/wiki/bikeshedding

Edit: Is my comment real bikeshedding?

[Snektron]

why bikeshed? I am not native speaker but it seems strange to ask for bikeshed:

Its a joke about that syntax discussions are often long with many drive-by comments

[rohlem]

I currently don't expect to use this feature, so not much weight behind my words, but anyway:
It seems worthwhile enough to have this just for the C interoperability, the rest seems like lower priority (although, as soon as someone confidently implements it, why not).

My only strong opinions reading this:

• Fields are strictly unsurprising in Zig. (Only maybe union fields, to people who don't know what unions are).
  In contrast, builtin @ functions are allowed to do many things.
  Therefore I'd strongly prefer to read code using "synthetic field names" .real and .imag/.imaginary for access, rather than @real and @imag.
  (I think I'd also prefer them for construction/initialization (EDIT: to make that work exactly like with structs, e.g.@Complex(f32){.real = 0, .imag = 3}), but that might be more subjective aesthetics.)
• I'm pretty surprised to hear about a definition of relational operators <, <=, >= and > for complex numbers.
  Okay, in mathematics it's not unheard of that you can redefine everything however you want.
  For a programming language I'd expect only a universally agreed-upon meaning to be exposed via builtins.
  I would also expect for == and != to mean full by-member equality in this case (and not be linked to abs or some other reductive operation).
  But having a situation where <=/>= might be true even though == isn't doesn't seem "obvious" / readable for newcomers.
  Therefore I'd rather advocate for equality to mean by-member and other relations to be provided only via named (std) functions.

[recursivetree]

My take on a few things:

• complex number creation via @complex(a,b) is in my opinion easier to understand and read than a+b*@I()
• I think accessing imaginary and real parts via @real and @imag makes it clear that it is a complex number. z.imag looks like a struct and not a complex number.
• I'd leave away any operator that's kind of vaguely defined, like % or greater than/less than comparisons. We can always add them later on. if you really need them, you can always implement them yourself using @real() and @imag() or whatever it will turn out to be.
• Random thought: What about complex numbers that use different types for real and imaginary parts? E.g. a f32 for the real part and a f64 for the imaginary part. I'm not sure where it will be useful and it is probably a bit out of scope, but I thought I'll throw it in.
• I agree with rohlem that == and != should be by member equality.

[ehaas]

For complex number literals / initializing complex numbers: My suggestions would be to have pub const I: @Complex(comptime_float) = .{0, 1}; in std.math (or something like that). Then people can put const I = std.math.I;at the top of their file and initialize complex numbers with var z: @Complex(f64) = 1 + 2*I;. This basically looks the same as C99, where I is a macro defined in <complex.h>

The tuple syntax and/or @complex builtin is also important because the method listed above can't be used to construct a complex number with a negative zero (-0.0) imaginary part without special casing (e.g. 5.0 + (-0.0)*I -> (5.0 + 0.0i) + (0.0 + -0.0i) -> 5.0 + 0.0i. Instead you would do @as(@Complex(f64), .{5, -0.0}) or @complex(5, -0.0). This corresponds to the CMPLX macro in C99.

The existence of a @complex builtin to construct complex numbers is nice from a "type fewer characters" perspective but I don't think it's necessary since it would do the same thing as @as(@Complex(f64), .{1, 2}) (and RLS would allow omitting the @as most of the time)

[recursivetree]

Another reason why I'm clearly not in favor of const z = a + b*I or const z = a + b*@I() is readability to people who don't know complex numbers. We can assume that anyone commenting in here has a basic understanding of what complex numbers are. Still, they are a fairly advanced topic that not everyone at my school gets taught. Therefore, they will most likely assume with const z = a + b*I that z is a float or int and not a complex. With @complex(a,b), they don't understand what's going on any better, but at least they know something is going on that they don't understand.

[BlueAlmost]

Thanks for the work on this. This looks quite nice. Some reflections:

*) Regarding construction: @complex(1, 1) has the cleanest and clearest look.

*) Regarding @real(x) vs. x.re access. Both seem quite reasonable.

The x.re syntax emphasises the "struct" nature of complex x.

• Does this struct emphasis introduce the possibility of any user confusion,
  with the more generic struct?? (i.e. any possible dangers there???).

• Or, does it help with clarity of the actual layout?

If we eventually choose to go with construction using @complex(1, 1) syntax, then using
the @real(x) harmonises quite nicely with that. Also, both accentuate the "built-in
nature of the complex type".

(The tuple-based construction var z: @complex(f32) = .{1, 1}; seems a bit ugly to me.)

*) Coercion rule: T -> Complex(T) looks very nice.

*) Regarding relational operators: Mathematically, complex numbers are NOT an ordered field.
In essence, it is impossible to define relational operators in a mathematically reasonable
way. Doing a "hidden" use of the modulus (a.k.a. abs) operator is dangerous IMHO. Forcing the programmer
to use the abs function, if that is what they intend for comparison, is the correct way.

Related to this, the hypot function's nasty computational burden is the sqrt part
of the computation. A good built-it for hypotSqr would relieve this heavy part. Or
perhaps call it abs2, or absSqr?

*) Regarding support for complex numbers:

• @max, @min have no meaning in a field without order. (see above).
  these must be flagged as a compile error.

[rohlem]

Quick note that we have precedence in slices []T being a builtin type that in usage works like a struct with members len and ptr, which are accessible as such. I don't see any reason for treating complex numbers differently in this regard.
IMHO a source of errors would be common syntax doing different things (than expected) depending on context. Sharing syntax for the same concept across types (builtin vs user-provided struct) seems more like a feature to me.
(But please provide counter-examples / footguns if you do happen to think of some.)

[Snektron]

Thanks all for your comments

I'm pretty surprised to hear about a definition of relational operators <, <=, >= and > for complex numbers.

But having a situation where <=/>= might be true even though == isn't doesn't seem "obvious" / readable for newcomers.

@max, @min have no meaning in a field without order. (see above).
these must be flagged as a compile error.

Ok, we should leave these away citing unclear mathematical definition.

I would also expect for == and != to mean full by-member equality in this case (and not be linked to abs or some other reductive operation).

I agree with rohlem that == and != should be by member equality.

Yep, that was the intention. In hindsight it seems obvious that these things interfere. I've adapted the proposal to highlight the intended semantics.

Random thought: What about complex numbers that use different types for real and imaginary parts? E.g. a f32 for the real part and a f64 for the imaginary part. I'm not sure where it will be useful and it is probably a bit out of scope, but I thought I'll throw it in.

Out of scope of this proposal. We are specifically interested in C99-compatible complex numbers. I think for anything that's more complicated, like mixed precision or SoA, the user should create a custom struct and convert that to a built-in complex number before processing (or write out the operation manually all together).

For complex number literals / initializing complex numbers: My suggestions would be to have pub const I: @complex(comptime_float) = .{0, 1};

That is a nice addition that is also a relatively minor change and fits well with the other constants that are already defined in std.math. That could even be added with the current complex number implementation :).

Related to this, the hypot function's nasty computational burden is the sqrt part
of the computation. A good built-it for hypotSqr would relieve this heavy part. Or
perhaps call it abs2, or absSqr?

Yes, I'm aware. I think it would be unexpected for the user if the main way to write sqrt in Zig does not work with complex numbers though, so implementation here feels appropriate. I think anything like hypotSqr should be implemented in library code, though.

Quick note that we have precedence in slices []T being a builtin type that in usage works like a struct with members len and ptr, which are accessible as such. I don't see any reason for treating complex numbers differently in this regard.

That's a good argument.

[presentfactory]

Is there any reason why quaternions are out of scope, since I do feel like if you're going to have logic for one type of extended number system like that then generalizing it so it would work for others like quaternions might be useful. Though I mostly just say that because quaternions are used a good bit in games and would be nice to have native support for like other kinds of vector math, I don't think anyone uses stuff like octonions for anything.

[Snektron]

Yes, the reason is that we are currently mainly focussing on C99 compatibility and sneaking in native complex operations as an added bonus. It does not rule out quaternions in Zig in general, just that they are not part of this proposal.

[BlueAlmost]

Small comment: I noticed that issue #16026, is now accepted for release 0.12, so perhaps our above use of @fabs, should instead be thought of as using @abs.

[silver-signal]

A builtin type function should return a type, not a typed value. I feel that @Complex(T: type) is a good fit, but @complex(re, im) is inconsistent. If a builtin is needed, it should have a different name like @complexVal(re, im) (or something else).

[expikr]

A builtin type function should return a type, not a typed value. I feel that @Complex(T: type) is a good fit, but @complex(re, im) is inconsistent. If a builtin is needed, it should have a different name like @complexVal(re, im) (or something else).

@Complex(T: type){re, im}

[expikr]

• Standard library support. This means implementing support for complex numbers in std.math where required. In time, std.math.Complex would be deprecated, and uses are replaced by the new type.

I have made a proof of concept here #19218

[expikr]

I’d like to propose a new language construct "builtin methods" for primitive types, which kills many birds with this one stone, including:

• vector swizzling
• @splat peer type resolution
• Builtin complex type field access and construction syntax

Rationale explained in the issue referenced above, quoted below:

@splat in its current form is useful for pre-computing a reusable SIMD vector that will be applied many times. Using it as a one-off scaling operation results in the verbose problem above.

But unequivocally there will certainly be many use cases of one-off scaling of a vector, and so we have a problem of using a tool in a situation for which it was not optimally designed for.

So really the crux of the issue comes down to needing a syntax to express a one-off operations on a vector as an expession.

My suggestion would be to add builtin methods to vector types, e.g. vec.@scale(s)

Notationally this introduces a new category of builtins easily distinguishable from @"name" style decls in the absence of quotes, indicating that it is part of the language itself rather than some custom definition.

This syntax might also be useful for built-in complex numbers (#16278 (comment)) where field access via z.@re and z.@im clearly distinguishes it from struct field access, which would otherwise indicate that they should be initialized with @Complex(F){.re=x, .im=y} rather than @Complex(F){x,y}

So the linked ray-tracing example might end up looking like:

    return (uv - n.@scale(dt)).@scale(ni_over_nt) - n.@scale(@sqrt(discriminant)) ;

Sample of what it might look like:

var c = @Complex(f32){0,0};

c.@re += 7;
c.@im += 8;

var cs = @Vector(@Complex(f32),3){
    .{1,2},
    .{3,4},
    .{5,6},
};

var ds = cs.@scale(c);
cs *= ds;

const es = cs.@swiz(1,1,2);

[rohlem]

@expikr Is there some rationale for why complex number fields .@re and .@im have the @ in your example when slice fields .len and .ptr don't have them in status-quo?
In my eyes as language-provided generic types they're in the same boat.

[expikr]

Hmm, true. I forgot about slice. In that case then I agree it's nicer to drop the @ for builtin attributes and reserve that notation for functions and methods that are language builtins.

[expikr]

A type suffix, for example im: 1+1im would return the equivalent complex value. Since Zig does not have any other numeric suffixes, I think this option is not the best.

Semantically 1.23456e-789 already sort-of has an annotating number suffix, so I think strictly speaking there is precedence.

[ilayn]

As a long time lurker, I just wanted to add that a typical operation on complex numbers is conjugation that is a + bi -> a-bi. Typically a method that goes by the name @conj. Having complex numbers as built-ins would be indeed fantastic for porting numerical code to Zig.