balbit's blog

By balbit, history, 5 years ago, In English

Here is a normal implementation of Pollard's Rho algorithm.

ll c = 1;
ll g(ll x, ll n){
    return (x*x+c)%n;
}

ll po(ll n){
    ll x = 2, y = 2, d = 1;
    while (d==1){
        x = g(x,n); y = g(g(y,n),n);
        d = __gcd(llabs(x-y),n);
    }
    if (d==n) return -1;
    return d;
}

However, the product x*x will overflow if n is too big (outside of int range). Is there a way to do the multiplication quickly (without using bigint or adding another log factor), or replacing the polynomial with some other function?

  • Vote: I like it
  • +7
  • Vote: I do not like it

| Write comment?
»
5 years ago, # |
  Vote: I like it +18 Vote: I do not like it

Try using __int128 if compiler supports it.

»
5 years ago, # |
  Vote: I like it +56 Vote: I do not like it

$$$ab \mod n = ab - \lfloor ab/n \rfloor n$$$. Calculate $$$ab/n$$$ in long double, then cast it to long long, it will differ from the real value of $$$\lfloor ab/n \rfloor$$$ by $$$\pm 1$$$. Let the value you got be $$$d$$$. Then calculate $$$ab - dn$$$ in long long. You will get the answer $$$\pm n$$$. Of course, technically long long overflow is undefined behavior, but in most cases it behaves just like multiplying modulo $$$2^{64}$$$ so it's not a big deal. You can always use unsigned long long to be extra sure, but then you need to be careful with negative values.

  • »
    »
    5 years ago, # ^ |
      Vote: I like it 0 Vote: I do not like it

    Why can't it differ by 2? Can you prove that the difference is not bigger than 1?

    • »
      »
      »
      5 years ago, # ^ |
      Rev. 2   Vote: I like it +28 Vote: I do not like it

      Assumes $$$a<n$$$ and $$$b<n$$$, $$$\lfloor ab/n \rfloor <n$$$, and long double has 64 bits for mantissa, so it would be (almost) enough precision.

      Proof sketch: ($$$c_1$$$ and $$$c_2$$$ are some "small enough" constants)

      Assumes the computer will always round to the nearest representable long double value. Note that conversion from long long to long double does not lose information (because long double mantissa is 64 bits long).

      Therefore the relative error of $$$ab$$$ does not exceed $$$2^{-64} \times c_1$$$, and the relative error of $$$ab/n$$$ doesn't exceed $$$2^{-64} \times c_2$$$, and (because $$$ab/n < n$$$) the absolute error doesn't exceed $$$n \times 2^{-64} \times c_2 < 1$$$. (QED)

»
5 years ago, # |
  Vote: I like it +18 Vote: I do not like it

You can compute this function modulo n in O(log x) using something like binary powering. You need to use that x * a = (x * (a — 1) + x) % n if a is odd and x * a = (2 * (x * (a / 2)) % n if a is even.

»
5 years ago, # |
  Vote: I like it -15 Vote: I do not like it

without using bigint

That's a mistake. Bigint is less bug-prone than obscure tricks in this case.

  • »
    »
    5 years ago, # ^ |
      Vote: I like it +10 Vote: I do not like it

    If you can use prepared template, then both are not bug-prone. If you can't, big int may be more bug prone because it's more complex. (In this case you need to implement modulus)

    Still, the tricks described here are pretty simple, and in the unlikely case you can't get it correct in about 5 minutes, you can switch to big int without too much time wasted...

    • »
      »
      »
      5 years ago, # ^ |
      Rev. 2   Vote: I like it -14 Vote: I do not like it

      You won't always have a prepared template for a specific thing, or you can be in a situation where you need to modify your template a bit. If you don't know how, you're screwed.

      You should know how to code a simple bigint — bruteforce sum, difference, product — by hand quickly, always without exception. If your problem can be solved by using five 64-bit ints to store a 128-bit number, that's what you should do because it wastes as little time as possible.

»
5 years ago, # |
  Vote: I like it 0 Vote: I do not like it

what this algo should do? on n=21, 22 and 25 (and others examples) it returns -1

»
5 years ago, # |
  Vote: I like it 0 Vote: I do not like it

Try checking out the source code for function "multiplyHigh" in java.lang.Math