mulmod failing

[ATH]

I thought I used this well known mulmod code in quite a few programs, so I would have noticed this before: It seems to fail when a and b are large and c is small in "a*b (mod c)". By failing I mean the program is stuck in this function or it takes a very long time compared to what it should.

For example a=245673636173; b=245673636173; and c=2047; and it also "fail" when a!=b.


I thought I did a lot of tests previously with many different values a,b,c < 2^64;

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t c) {
  uint64_t d; /* to hold the result of a*b mod c */
  /* calculates a*b mod c, stores result in d */
  asm ("mov %1, %%rax;"        /* put a into rax */
       "mul %2;"               /* mul a*b -> rdx:rax */
       "div %3;"               /* (a*b)/c -> quot in rax remainder in rdx */
       "mov %%rdx, %0;"        /* store result in d */
       :"=r"(d)                /* output */
       :"r"(a), "r"(b), "r"(c) /* input */
       :"%rax", "%rdx"         /* clobbered registers */
      );
  return d;
}

[danaj]

In my code I have this comment:

  /* GCC on a 64-bit Intel x86, help from WraithX and Wojciech Izykowski */
  /* Beware: if (a*b)/c > 2^64, there will be an FP exception */
  static inline uint64_t _mulmod(uint64_t a, uint64_t b, uint64_t n) {
    uint64_t d, dummy;                    /* d will get a*b mod c */
    asm ("mulq %3\n\t"              /* mul a*b -> rdx:rax */
         "divq %4\n\t"              /* (a*b)/c -> quot in rax remainder in rdx */
         :"=a"(dummy), "=&d"(d)     /* output */
         :"a"(a), "rm"(b), "rm"(n)  /* input */
         :"cc"                      /* mulq and divq can set conditions */
        );
    return d;
  }

It seems the exact behavior varies -- my Mac hangs, while some other machines I have exit with an FP exception.

[rogue]

It is safest to compute a%c and b%c first. Yes, it is slower, but in the code I work with I typically need the value of a%c and b%c.

[ATH]

Thank you. You are both correct. I found my example was working when I lowered a down below the limit a*a/2047 < 2^64.

What got me confused was, that I was sure I had tested mulmod with random values up to 2^64 for all 3 variables. But now I found that old test code, and I actually did a%c and b%c before mulmod back then, and I completely forgot

[retina]

With the values you gave your code should crash with a divide-by-zero* exception. Unless the calling code has installed an exception handler or something.

Any time the high order word is larger or equal than the divisor the result from div can't fit into a single register.

245,673,636,173 x 245,673,636,173 = 60,355,535,510,463,574,085,929

60,355,535,510,463,574,085,929 / 2,047 = 29,484,873,234,227,442,152 & 785 remainder

Note that 29,484,873,234,227,442,152 cannot fit into a single 64-bit register.

This can also be detected by observing that the high order word of 60,355,535,510,463,574,085,929 is 3,271. So even before the divide is done we can know that it will fail because 3,271 >= 2,047.

* Even though the divisor is not actually zero you will still get the divide-by-zero exception whenever the result is too large.