21 Days Training // Day 12

T4. leap

Description(Simplified) :
Find the mininum positive interger x that satisfies \frac{x(x+1)}{2} \mid N , i.e. {x(x+1)} \mid 2N .

Wrong idea :

Enumerate the factors of 2N, and find 2 factors p_1, p_2 that satisfies p_2 \equiv \pm 1 \pmod{p_1} .
If, p_2 \equiv 1 \pmod{p_1} , x=p_2-1=kp_1 .
If, p_2 \equiv -1 \pmod{p_1} , x=p_2+1=kp_1 .

For all possiblities, find the mininum one.

Bugs :
This idea did not consider the condition that x=ap_1, x+1=bp_2 .

(Right) idea :
When enumerating all factors, we need to find the mininum positive inteeger a, b that satisfies ap_1+1=bp_2 , i.e. bp_2-ap_1=1 .

For this equation use exgcd to solve it.
After that, calculate x and find the mininum one.
Specially, when p_1=1, current x equals to N-1.

Code :

#include <bits/stdc++.h>
using namespace std;
const long long LIM=1000000;
bitset<LIM+10> isprime, vis;
long long plist[80000], dpf[30];
int pcount=0, dpfcount=0;
set<long long> fac;
map<long long, int> pfac;

void exgcd(long long a, long long b, long long& x, long long& y) {
    if (b == 0) {
        x = 1, y = 0;
        return;
    }
    exgcd(b, a % b, y, x);
    y -= a / b * x;
}

void preprocess() {
    isprime[2]=1, vis[2]=1;
    plist[pcount++]=2;
    for (long long j=4;j<=LIM+10;j+=2) vis[j]=1;
    for (long long i=3;i<=LIM+10;i++) {
        if (vis[i]) continue;
        plist[pcount++]=i;
        isprime[i]=1, vis[i]=1;
        for (long long j=2*i;j<=LIM+10;j+=i) vis[j]=1;
    }
}

void dfs(int L, long long cur) {
    if (L==dpfcount) {
        fac.insert(cur);
        return ;
    }
    long long bk=cur;
    for (int i=0;i<=pfac[dpf[L]];i++) {
        dfs(L+1, bk);
        bk*=dpf[L];
    }
}

void leap() {
    long long N;cin >> N;
    N*=2;
    long long bk=N;
    pfac.clear();
    dpfcount=0;
    for (int i=0;i<pcount;i++) {
        if (!(bk%plist[i])) {
            dpf[dpfcount++]=plist[i];
            pfac[plist[i]]=0;
            while (!(bk%plist[i])) pfac[plist[i]]++, bk/=plist[i];
        }
        if (bk==1) break;
    }
    if (bk!=1) pfac[bk]=1, dpf[dpfcount++]=bk;
    fac.clear();
    dfs(0,1);
    long long ans=LONG_LONG_MAX;
    for (set<long long>::iterator it=fac.begin();it!=fac.end();it++) {
        long long minfac=*it, maxfac;
        if (minfac>sqrt(N)) break;
        maxfac = N/minfac;
        if (minfac==1) {
            ans = min(ans, N-1);
            continue;
        }
        if (__gcd(maxfac, minfac)!=1) continue;
        //minfac*x-maxfac*y=1
        //maxfac*x-minfac*y=1
        long long x1, y1, x2, y2;
        exgcd(minfac, -maxfac, x1, y1);
        exgcd(maxfac, -minfac, x2, y2);
        ans = min(ans, minfac*x1);
        ans = min(ans, maxfac*x2);
    }
    cout << ans << endl;
}

int main() {
    preprocess();
    int C;cin >> C;
    for (int i=0;i<C;i++) leap();
    return 0;
}

Optimize :

1. When \gcd(p_1, p_2) \ne 1, we can directly skip it. (Line 70).
2. Preprocess : Calculate all primes between 1 and 1e6. (Line 20~30).