Skip to content

Complex remainder doesn't find the nearest multiplicator of the divider due to approximations #445

Description

@Astrac

I noticed that if I calculate cetain remainders I get results that do not agree with what I calculate using other tools; for example, 10 % (5 + 7i) gives 10 as a result with this library but other tools (e.g. wolfram alpha) give me -2 - 2i. To corroborate wolfram's result just notice that |10| = 10 and |5 + 7i| ~ 8.6, while the modulus of the reminder should be smaller than the one of the divisor.

I found that the issue causing this is the way the remainder is computed using only integer types: the implementation always find the biggest multiplier of the divisor that is smaller than dividend, not the nearest mutiplier. If I rewrite the algorithm as follows I get the very same results as wolfram and other code of mine that depends on this calculation behaves correctly:

fn calc_mod(n: Complex<i32>, modulus: Complex<i32>) -> Complex<i32> {
    // Transform to f32
    let nf = Complex::new(n.re as f32, n.im as f32);
    let modf = Complex::new(modulus.re as f32, modulus.im as f32);
    
    // The usual complex division formula but applied to Complex<f32> instead of <Complex<i32>
    let div = (nf * modf.conj()) / modf.norm_sqr();

    // use round() to ensure closest multiplier is found
    let rounded = Complex::new(div.re.round() as i32, div.im.round() as i32); 
    let nearest_mult = rounded * modulus;

    return n - nearest_mult;
}

Activity

  1. cuviper commented on Dec 15, 2025

    @cuviper
    Member

    The definition we've chosen is consistent (for integers) such that (a / b) * b + a % b == a.

    I acknowledge that this isn't necessarily the most useful definition, but I think that also depends on your domain and how you're going to use the result -- similar to the tension between plain integer Rem vs rem_euclid.

    Maybe we need a similar specialized method for Complex -- if that can be written in a generic way.

  2. cuviper commented on Dec 15, 2025

    @cuviper
    Member

    I found some history in #327, which did try to be more clever at first.

    One of the points I made then was, "I think it's a fundamental property that operations on Complex values with zero imaginary should match the equivalent operation on the real part alone." I still think that's desirable, at least for the normal operators -- it definitely doesn't hold for other methods like cbrt of negatives.

  3. Astrac commented on Dec 16, 2025

    @Astrac
    Author

    Thank you for getting back to me - what you say makes sense in particular with the rem vs rem_euclid thing that I didn't notice before.

    It seems to me that a point could be made that what kind of modulus you get should depend on what number set your using in the first place; in particular, the current implementation is perfect for positive integers but when we extend it to the relative numbers we have the possibility to count back from the nearest multiplier of the divisor. If we reason like this then the non-complex remainder operation should give the same result as the complex one. For example, 10 % 13 is 10 restricted to natural numbers but it is also (and in some cases more correctly) -3 if we are dealing with relative numbers.

    I'm not sure that making this depending on whether we are using e.g. i32 or u32 makes it even more confusing though, so maybe having the two rem and rem_euclid functions and explicitly choosing is a good compromise.

  4. cuviper commented on Dec 16, 2025

    @cuviper
    Member

    If we reason like this then the non-complex remainder operation should give the same result as the complex one. For example, 10 % 13 is 10 restricted to natural numbers but it is also (and in some cases more correctly) -3 if we are dealing with relative numbers.

    Interesting -- I don't know of any scalar implementation that gives that result, so at the very least we need to give it a new name. Maybe rem_relative or rem_nearest, or something like rem_minimal to indicate that it's the remainder with the least magnitude. (and maybe the least [0, 2π) argument in case of tied magnitude?)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions