http://acm.hdu.edu.cn/showproblem.php?pid=3930題目意思很簡單 對于方程 x^k = b mod p，給出k，b，p，求所有的x ( 0<= x < p )，題目的數據范圍很惡心，其實沒有那么大，只有10^12那么大，可以打素數表了親，再大的話，只能rho來分解了。。。 1) 先暴力求p的原根g 2) 大步小步求g^t1 = b mod p 3) 則g^(t1+n*t2) = b mod p, t2 = p - 1 4) 設x = g^y mod p, x^k = (g^y)^k = g^(yk) = g^(t1+n*t2) 那么就是求同余方程 yk = t1( mod t2 ) 求出y之后帶到x = g^y mod p,解出x
1  /**//*
2 author: AmazingCaddy
3 time: 2011/8/13 13:41
4 */
5 #include <cstdio>
6#include <complex>
7#include <cstdlib>
8#include <iostream>
9#include <cstring>
10#include <string>
11#include <algorithm>
12#include <cmath>
13#include <ctime>
14#include <vector>
15#include <map>
16#include <queue>
17using namespace std;
18
19//typedef __int64 ll;
20typedef long long ll;
21
22const int maxn = 1000004;
23int vis[ maxn ], p[ maxn ];
24int plen, flen;
25ll a[ 64 ], b[ 64 ];
26ll k, m, newx, g;
27
28void prime( )  {
29 memset( vis, 0, sizeof( vis ) );
30 plen = 0;
31  for( ll i = 2, k = 4; i < maxn; ++i, k += i + i - 1 ) {
32 if( !vis[ i ] ) {
33 p[ plen++ ] = i;
34 if( k < maxn ) {
35 for( ll j = k; j < maxn; j += i ) {
36 vis[ j ] = 1;
37 }
38 }
39 }
40 }
41}
42
43ll mul_mod( ll a, ll b, ll mod ) {
44 ll ans = 0, d = a % mod;
45 while( b ) {
46 if( b & 1 ) {
47 ans += d;
48 if( ans > mod ) {
49 ans -= mod;
50 }
51 }
52 d += d;
53 if( d > mod ) {
54 d -= mod;
55 }
56 b >>= 1;
57 }
58 return ans;
59}
60
61ll pow_mod( ll a, ll b, ll mod ) {
62 ll ans = 1, d = a % mod;
63 while( b ) {
64 if( b & 1 ) {
65 ans = mul_mod( ans, d, mod );
66 }
67 d = mul_mod( d, d, mod );
68 b >>= 1;
69 }
70 return ans;
71}
72
73void factor( ll n ) {
74 flen = 0;
75 for( int i = 0;(ll) p[ i ] * p[ i ] <= n; i++ ) {
76 if( n % p[ i ] == 0 ) {
77 for( b[ flen ] = 0; n % p[ i ] == 0; ++b[ flen ], n /= p[ i ] );
78 a[ flen++ ] = p[ i ];
79 }
80 }
81 if( n > 1 ) a[ flen ] = n, b[ flen++ ] = 1;
82}
83
84void ex_gcd( ll a, ll b, ll & d, ll &x, ll &y ) {
85 if( !b ) {
86 x = 1, y = 0, d = a;
87 }
88 else {
89 ex_gcd( b, a % b, d, y, x );
90 y -= x * ( a / b );
91 }
92}
93
94ll Inv( ll n, ll p ) {
95 return pow_mod( n, p - 2, p );
96}
97
98bool dfs( int dep, ll t ) {
99 if( dep == flen ) {
100 ll ans = pow_mod( g, t, m );
101 if( ans == 1 && t != m - 1 ) return false;
102 return true;
103 }
104 ll tmp = 1;
105 for( int i = 0; i <= b[ dep ]; i++ ) {
106 if( !dfs( dep + 1, t * tmp ) ) return false;
107 tmp *= a[ dep ];
108 }
109 return true;
110}
111
112void find_g( ) {
113 factor( m - 1 );
114 for( g = 2; ; g++ ) {
115 if( dfs( 0, 1 ) ) break;
116 }
117}
118
119// a^x = b ( mod p ) find x, p is prime
120ll log_x( ll a, ll b, ll p ) {
121 map< ll, int > x;
122 ll z = (ll)ceil( sqrt( p * 1.0 ) );
123 ll v = Inv( pow_mod( a, z, p ), p );
124 ll e = 1;
125 x[ 1 ] = 0;
126 for( int i = 1; i < z; i++ ) {
127 e = mul_mod( e, a, p );
128 if( !x.count( e ) ) {
129 x[ e ] = i;
130 }
131 }
132 for( int i = 0; i < z; i++ ) {
133 if( x.count( b ) ) {
134 return i * z + x[ b ];
135 }
136 b = mul_mod( b, v, p );
137 }
138 return -1;
139}
140
141ll sol[ 1000 ];
142void solve( ll a, ll b, ll n ) {
143 ll d, x, y;
144 ex_gcd( a, n, d, x, y );
145 if( b % d ) {
146 printf( "-1\n" );
147 }
148 else {
149 //ax+ny=d, so a'x+n'y=1, x is the inverse of a' mod n'
150 n /= d, b /= d;
151 sol[ 0 ] = ( x * b % n + n ) % n; // first soltion
152 for( int i = 1; i < d; i++ ) {
153 sol[ i ] = sol[ i - 1 ] + n;
154 }
155 for( int i = 0; i < d; i++ ) {
156 sol[ i ] = pow_mod( g, sol[ i ], m );
157 }
158 sort( sol, sol + d );
159 for( int i = 0; i < d; i++ ) {
160 printf( "%I64d\n", sol[ i ] );
161 }
162 }
163}
164
165int main(int argc, char *argv[])
166{
167// freopen( "in.in", "r", stdin );
168// freopen( "out", "w", stdout );
169 int cas = 1;
170 prime( );
171 while( scanf( "%I64d%I64d%I64d", &k, &m, &newx ) != EOF ) {
172 find_g( );
173 ll t1 = log_x( g, newx, m );
174 ll t2 = m - 1;
175 printf( "case%d:\n", cas++ );
176 solve( k, t1, t2 );
177 }
178 return 0;
179}
180
|