How does Chapel do its complex arithmetic

Does it do its own or does it map to some complex arithmetic handler in LLVM?


Different operations can be implemented differently, but I will describe one just to give you a sense. Taking a quick look at +, it's implemented as calling a primitive which is code-generated to call complexAdd128 (for 128-bit complex numbers) which is present in the runtime headers and just contains return c1+c2;. The result should be that we're relying on clang to convert + on complex numbers into LLVM IR, but that then we are inlining this into the generated LLVM IR from the Chapel program.

Are you curious about a specific operation?

Sorry for the delay. I missed this.

So, it is all done in the back-end, not Chapel per se.

I want to figure out how to achieve

any real component * an imaginary component of 0.0 = 0.0 i


real component of 0.0 * any imaginary component = 0.0 i

where any includes an Infinity.

Note that in IEEE 754

Intinity (as a real number) * 0.0 (as a real number) = NaN

So it looks like it requires a change to the back-end to do this in Chapel. Thanks.

Hi Damian --

We can implement such things in Chapel module code or in C code (that we inline). So if you are interested in making some particular changes, you might want to look at trying to write some of your own functions in Chapel module code to do it. If we get to the point of wanting to include such things in the standard library, we'll need to investigate to what degree using them performs better or worse than what we have today.

I'm also curious if you think this would necessarily be a different routine from the normal * because of the difference from IEEE 754 behavior that you noted.

I do not know yet.

HP's (now defunct?) C compiler implemented true imaginary numbers properly as defined as optional in the C standard. Nobody else's compiler did/does this properly. We now have the situation where Chapel also does this correctly after the recent fix.

However, given complex number such as

const x = (a, b):complex(128);
const y = (c, d):complex(128);
const z = x * y;

then will be NaN if say a is a zero and d is an Infinity because of existing IEEE 754 rules of multiplying two floating point numbers. Remember that IEEE 754 does not deal with complex numbers.

Unfortunately languages need to. So for the case where either the real of imaginary component is empty, there is no multiplication to be done (by Infinity or anything else). Every language I know of produces the mathematically incorrect number because it has way of thinking that 0.0 is nothingness. There is discussion going on currently about correcting that in some languages, and even introducing complex numbers into IEEE 754 although such discussion note and do not underestimate the magnitude of that latter task (and note that is may preclude it being done).

It would be nice if we had an option in Chapel to produce the correct mathematical result for both practical purposes for some types of physical problems and to show that it can be done in a modern language. I have not thought exactly how best to do it but it would be as an adjunct done after the multiply. For now, I just wanted to know what went on under the hood so to speak so as to better inform my own thinking or discussions with others.

Ultimately, we will have an option to a hardware multiply or fma, or a separate instruction to mandate that the operation

+/- Infinity * +/- Zero

is always 0.0 rather than NaN, i.e. the Zero is treated as emptiness rather than a numeric value. But do not expect that soon, and most probably not before the next IEEE 754 revision.

I think I know enough for now to think about it more intelligently.

I think what I am talking about is the concept of missing. But I am so NOT a person who works with missing that I might stop here. When I work in statistics which deals in that stuff, I work with a statistics guru to keep me from falling off a cliff.