A^x = D (mod P)
0 <= x <= M, here M is a given integer.
1 <= A, P < 2^31, 0 <= D < P, 1 <= M < 2^63
----------------------------------------------------------
裸拓展baby step giant step
先转成非拓展版，然后如果在转的过程中就出解了，则return 1,否则就找出=D的起点和循环节长度，然后直接求解。
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
const int Hash_MOD = 1000003;
typedef long long LL;
typedef unsigned long long ULL;
LL times;
struct link
{
LL link;
int val,next;
}es[100000];
struct triple
{
LL x,y,g;
triple(const LL _x = 0,const LL _y = 0,const LL _g = 0):
x(_x),y(_y),g(_g){}
};
int H[Hash_MOD * 2],ec;
triple exgcd(const LL a,const LL b)
{
if (!b) return triple(1,0,a);
const triple last(exgcd(b,a%b));
return triple(last.y,last.x - a / b * last.y,last.g);
}
LL rem_equ(const LL a,const LL b,const LL c)
{
//ax == c (mod b)
//ax mod b == c
const triple tmp(exgcd(a,b));
const LL MOD = (b / tmp.g);
return ((tmp.x * (c / tmp.g) % MOD) + MOD) % MOD;
}
LL find(const LL x)
{
for (int i = H[x%Hash_MOD]; i; i = es[i].next)
if (es[i].link == x) return es[i].val;
return -1;
}
void insert(LL x,const LL v)
{
if (find(x) != -1) return;
const LL key = x % Hash_MOD;
es[++ec].next = H[key];
H[key] = ec;
es[ec].link = x;
es[ec].val = v;
}
LL BSGS(LL A,LL P,LL D)
{
LL AA = 1 % P,x = 1 % P,MOD = P,DD = D,m = static_cast<LL>(std::ceil(std::sqrt((double)P)));
times = 0;
for (triple tmp; (tmp = exgcd(A,P)).g != 1; ) {
if (x == DD) return times++;
if (D % tmp.g) return -1;
P /= tmp.g;
D /= tmp.g;
(AA *= A / tmp.g) %= P;
++times;
(x *= A) %= MOD;
}
A %= P;
ec = 0;
memset(H,0,sizeof H);
LL tmp = 1 % P;
for (int i = 0; i < m; ++i,(tmp *= A) %= P)
insert(tmp,i);
x = 1 % P;
for (LL i = 0; i < m; ++i,(x *= tmp) %= P) {
const int j = find(rem_equ(AA*x%P,P,D));
if (j != -1) return i * m + j + times;
}
return -1;
}
LL getphi(LL x)
{
LL res = 1;
for (LL i = 2; i * i <= x; ++i)
if (x % i == 0) {
x /= i;
res *= i - 1;
while (x % i == 0) x /= i,res *= i;
}
if (x > 1) res *= x - 1;
return res;
}
LL power(LL a,LL k,LL MOD)
{
//a ^ k % MOD
LL res = 1 % MOD;
for (; k; k/=2,(a *= a) %= MOD)
if (k & 1) (res *= a) %= MOD;
return res;
}
LL calc_len(const LL start,const LL A,const LL P,const LL D)
{
//A ^ (start + *this) == A ^ start == D (mod P)
LL phi = getphi(P),res = phi;
for (LL i = 2; i * i <= phi; ++i) {
for (; phi % i == 0; phi /= i);
for (; res % i == 0 && power(A,start + res/i,P) == D; res /= i);
}
if (phi > 1)
for (; res % phi == 0 && power(A,start + res/phi,P) == D; res /= phi);
return res;
}
ULL work(const LL A,const LL P,const LL D,const ULL M)
{
LL start = BSGS(A,P,D);
//printf("%lld
",start);
if (start == -1 || start > M) return 0;
else if (start < times) return 1;
// LL phi=getphi(P);
// if (power(A,start + phi,P) != D) return 1;
ULL len = calc_len(start,A,P,D);
return (M - start) / len + 1;
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("3254.in","r",stdin);
freopen("3254.out","w",stdout);
#endif
LL A,P,D;
ULL M;
while (~scanf("%lld%lld%lld%llu",&A,&P,&D,&M))
printf("%llu
",work(A%P,P,D,M));
}