implement complex div/mul libcalls by folkertdev · Pull Request #5319 · rust-lang/miri

folkertdev · 2026-09-17T21:05:14Z

This needs rust-lang/rust#162832 to be merged and synced for the Div/Mul instances.

Maybe we can iterate on the implementation already, or if not this will just sit for a couple of days.

RalfJung · 2026-09-17T21:17:35Z

+fn complex_binop<'tcx, S: Semantics>(
+    this: &mut MiriInterpCx<'tcx>,
+    op: ComplexFloatBinop<S>,
+    [a, b, c, d]: &[OpTy<'tcx>; 4],


These names are not very helpful... at the very least please document them. But I assume they correspond to the real and imaginary parts of the left and right input, respectively -- the names should reflect that to some extend, like [a_re, a_im, b_re, b_im] or so.

View changes since the review

RalfJung · 2026-09-17T21:20:48Z

+    let ac = (a * c).value;
+    let bd = (b * d).value;
+    let ad = (a * d).value;
+    let bc = (b * c).value;
+
+    let z = Complex::new((ac - bd).value, (ad + bc).value);


Isn't this duplicating naive_complex_mul?

View changes since the review

yes, but the ac etc. variables are used below.

RalfJung · 2026-09-17T21:21:13Z

+        mag.copy_sign(x)
+    };
+
+    let mut recalc = false;


Why is this so imperative? What's the point of this mutable variable?

View changes since the review

It is translated directly from compiler-rt, e.g.

https://github.com/llvm/llvm-project/blob/5bff640ac99dc29ab61cb55bdfc9b992d09b5a85/compiler-rt/lib/builtins/muldc3.c

RalfJung · 2026-09-18T09:33:27Z

+    z.re = (z.re * IeeeFloat::<S>::INFINITY).value;
+    z.im = (z.im * IeeeFloat::<S>::INFINITY).value;
+
+    z


There's a lot going on in this function and none of it is explained. What exactly is the spec being implemented here? Why is "at least one non-NaN" sufficient to know that the "fast-path" is fine? What is the point of all the logic after that NaN check?

View changes since the review

This is specified in annex G of the C23 standard

https://www.open-std.org/jtc1/sc22/wg14/www/docs/n3220.pdf#page=556

I'll port over the comments, those are nice.

So I guess this is about this part of the spec?

Please explicitly explain how the code achieves that spec. Tests should also reference that spec to ensure they cover all these corner cases.

That document also says

It is less than ideal in that undue overflow and underflow could occur.

What is that about?

It is less than ideal in that undue overflow and underflow could occur.

I found https://arxiv.org/pdf/1210.4539 which presents some cases where the C specification algorithm still produces incorrect results. There is also https://www.open-std.org/jtc1/sc22/wg14/www/docs/n1399.htm which presents more incorrect results. It proposes an overhaul to complex division that has clearly not been adopted.

So far I don't even know what "correct" and "incorrect" mean. What is the specification for these operations? That's what the doc comment should explain. Without that, discussing correctness is pointless.

However it seems like correctness is not actually precisely defined and the functions are approximations. Doesn't that mean we should randomize them like we randomize sin, log, etc?

RalfJung · 2026-09-19T08:19:12Z

+///
+/// a+bi * c+di = ((ac + bd) + (bc - ad)i) / (c*c + d*d)
+///
+/// But with recovery of infinities if the above expresion results in NaN + NaNi.


Please explain what "recovery of infinities" means, or reference a document that explains this.

View changes since the review

I've added some examples

RalfJung · 2026-09-19T08:19:25Z

+    mut a: IeeeFloat<S>,
+    mut b: IeeeFloat<S>,
+    mut c: IeeeFloat<S>,
+    mut d: IeeeFloat<S>,
+) -> Complex<IeeeFloat<S>> {
+    let ac = (a * c).value;
+    let bd = (b * d).value;
+    let ad = (a * d).value;
+    let bc = (b * c).value;


These names are still completely meaningless.^^

View changes since the review

They represent a + bi and c + di, and ac is the product a * c? They are abstract but relate to the doc comment and this is, I think, the standard way to write this down?

in other words: I don't think using lhs_re, lhs_re_mul_rhs_re etc is substantially clearer.

Hm, I guess we're going to link to https://www.open-std.org/jtc1/sc22/wg14/www/docs/n3220.pdf#page=556 which also uses a, b, c, d...

Yeah, I've linked that now. Having done more reading, it's really not "the best" way of implementing these operations (especially division), but for these libcalls the C reference algorithm is what we should implement.

RalfJung · 2026-09-19T08:20:47Z

+    let mut z = naive_complex_mul(a, b, c, d);
+    z.re = (z.re * IeeeFloat::<S>::INFINITY).value;
+    z.im = (z.im * IeeeFloat::<S>::INFINITY).value;


What is this doing?

View changes since the review

RalfJung · 2026-09-19T08:23:12Z

+    // Scale the inputs
+    let max = IeeeFloat::<S>::max(c.abs(), d.abs());
+    let mut neg_ilogbw = 0;
+    if max.is_finite() && max != IeeeFloat::<S>::ZERO {
+        neg_ilogbw = max.ilogb().checked_neg().unwrap();
+        c = c.scalbn(neg_ilogbw);
+        d = d.scalbn(neg_ilogbw);
+    }


What is happening here? Scale the inputs how and why?

View changes since the review

My theory is that it's about preventing overflow. I'll have to look into that more though.

Hmm, the most I could find about it is from

In the ISO/IEC 9899:TC3 C99 standard [9], section G.5.1, “Multiplicative operators”, the authors present a _Cdivd function in the C language which implements complex division. According to Kahan [10] (in the Appendix, “Over/Underflow Undermines Complex Number Division in Java”), this code is due to Jim Thomas and Fred Tydeman.

The algorithm is based on a prescaling of c and d, such that their base-2 exponent is made zero. This scales the denominator c² + d² and avoids most overflows or underflows. The post-scaling produces then the correct result. This scaling is based on a power of 2, which avoids rounding.

Only in the case of an IEEE exception, the algorithm recomputes the division, taking into account for NaNs and Infinities. The code does not defend against overflow and underflow in the calculation of the numerator.

We will show later in this paper that the C99 algorithm performs nearly as well as Smith’s algorithm, but fails in many cases.

RalfJung · 2026-09-19T18:23:46Z

+        // Replace infinities with a signed unit value, and NaN with zero.
+        // This avoids invalid infinity arithmetic.


We're doing the inf replacement on a, b and the NaN replacement on c, d. Why that? The comment makes it sounds like we're doing both replacements on all values.

View changes since the review

rustbot added the S-waiting-on-author Status: Waiting for the PR author to address review comments label Sep 17, 2026

RalfJung reviewed Sep 18, 2026

View reviewed changes

folkertdev force-pushed the complex-mul-div branch from d48a831 to 2a59386 Compare September 18, 2026 22:56

RalfJung reviewed Sep 19, 2026

View reviewed changes

folkertdev force-pushed the complex-mul-div branch from 2a59386 to 779d463 Compare September 19, 2026 18:07

RalfJung reviewed Sep 19, 2026

View reviewed changes

folkertdev force-pushed the complex-mul-div branch from 779d463 to 09b6791 Compare September 19, 2026 20:24

implement complex div/mul libcalls

2dd9e42

folkertdev force-pushed the complex-mul-div branch from 09b6791 to 2dd9e42 Compare September 19, 2026 20:28

		// Replace infinities with a signed unit value, and NaN with zero.
		// This avoids invalid infinity arithmetic.

Conversation

folkertdev commented Sep 17, 2026

Uh oh!

Uh oh!

RalfJung Sep 17, 2026 • edited by rustbot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

RalfJung Sep 17, 2026 • edited by rustbot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

RalfJung Sep 17, 2026 • edited by rustbot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

RalfJung Sep 18, 2026 • edited by rustbot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

RalfJung Sep 19, 2026 • edited by rustbot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

RalfJung Sep 19, 2026 • edited by rustbot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

RalfJung Sep 19, 2026 • edited by rustbot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

RalfJung Sep 19, 2026 • edited by rustbot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

RalfJung Sep 19, 2026 • edited by rustbot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

RalfJung Sep 17, 2026 •

edited by rustbot

Loading

RalfJung Sep 17, 2026 •

edited by rustbot

Loading

RalfJung Sep 17, 2026 •

edited by rustbot

Loading

RalfJung Sep 18, 2026 •

edited by rustbot

Loading

RalfJung Sep 19, 2026 •

edited by rustbot

Loading

RalfJung Sep 19, 2026 •

edited by rustbot

Loading

RalfJung Sep 19, 2026 •

edited by rustbot

Loading

RalfJung Sep 19, 2026 •

edited by rustbot

Loading

RalfJung Sep 19, 2026 •

edited by rustbot

Loading