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
|