### balbit's blog

By balbit, history, 16 months ago,

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?

• +7

 » 16 months ago, # |   +18 Try using __int128 if compiler supports it.
 » 16 months ago, # |   +56 $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.
•  » » 16 months ago, # ^ |   0 Why can't it differ by 2? Can you prove that the difference is not bigger than 1?
•  » » » 13 months ago, # ^ | ← Rev. 2 →   +28 Assumes $a  » 16 months ago, # | +18 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.  » 13 months ago, # | -15 without using bigint That's a mistake. Bigint is less bug-prone than obscure tricks in this case. •  » » 13 months ago, # ^ | +10 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... •  » » » 13 months ago, # ^ | ← Rev. 2 → -14 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.  » 13 months ago, # | 0 what this algo should do? on n=21, 22 and 25 (and others examples) it returns -1  » 13 months ago, # | 0 Try checking out the source code for function "multiplyHigh" in java.lang.Math •  » » 13 months ago, # ^ | 0 Even if you have the 128-bit representation of the product, it's not trivial to compute that value modulo mod.  » 13 months ago, # | ← Rev. 5 → 0 Another method: In case the modulus is fixed and a perfect square (such as$10^{18}\$), it's possible to be reasonably fast using integer arithmetic only. (It's easy to derive the function yourself. See function splitdigit in code below) About performance:Some notes about Pollard Rho being slow: https://codeforces.com/blog/entry/62287#comment-462676Benchmark: code//#pragma GCC optimize ("fast-math") #include #include #include #include #include #include using namespace std; uint64_t constexpr MOD=1000000000000000000ULL; uint64_t subquotient1(uint64_t a,uint64_t b){ auto q=uint64_t((long double)a*b/MOD); // error: -1, 0 or 1 auto out=a*b-q*MOD+MOD; // error: 0, MOD or 2*MOD if(out>=MOD){ out-=MOD; if(out>=MOD)out-=MOD; } return out; } uint64_t subquotient2(uint64_t a,uint64_t b){ auto q=uint64_t((long double)a*b/MOD-0.5); // error: 0 or 1 auto out=a*b-q*MOD; // error: 0 or MOD if(out>=MOD)out-=MOD; return out; } uint64_t subquotient3(uint64_t a,uint64_t b){ auto q=uint64_t(std::round((long double)a*b/MOD))-1; // error: 0 or 1 auto out=a*b-q*MOD; // error: 0 or MOD if(out>=MOD)out-=MOD; return out; } uint64_t subquotient4(uint64_t a,uint64_t b){ auto q=uint64_t((long double)a*b/MOD); auto out=a*b-q*MOD+MOD; if(out>=MOD)out-=MOD; if(out>=MOD)out-=MOD; return out; } uint64_t subquotient5(uint64_t a,uint64_t b){ auto q=uint64_t((long double)a*b/MOD); auto out=a*b-q*MOD+MOD; return out%MOD; } #ifdef __SIZEOF_INT128__ uint64_t useint128(uint64_t a,uint64_t b){ return uint64_t((unsigned __int128)a*b%MOD); } #endif uint64_t splitdigit(uint64_t ax,uint64_t bx){ // specific to MOD 1000000000000000000 (or any perfect square mod) auto ah=ax/1000000000,al=ax%1000000000; auto bh=bx/1000000000,bl=bx%1000000000; auto out=(ah*bl+al*bh)%1000000000*1000000000 + al*bl; if(out>=MOD)out-=MOD; return out; }  Benchmark codeint constexpr N_ITER=10000000; // compute 1234**(2**N_ITER) for all possible base values with N_ITER applications of f vector> benchresult; void bench(char const* name,auto f){ auto start=clock(); uint64_t result=1234; for(int i=0; i
 » 4 weeks ago, # |   -30 How is this function multiplying two numbers??? Any article explaining this algo ?? , i found some factorisation algos with this name.Btw ,what is most optimal algo to multiply two large float numbers in c++? Any idea