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, 0sizeof( 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    forint 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    forint 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( 01 ) ) 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    forint i = 1; i < z; i++ ) {
127        e = mul_mod( e, a, p );
128        if!x.count( e ) ) {
129            x[ e ] = i;
130        }

131    }

132    forint 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        forint i = 1; i < d; i++ ) {
153            sol[ i ] = sol[ i - 1 ] + n;
154        }

155        forint i = 0; i < d; i++ ) {
156            sol[ i ] = pow_mod( g, sol[ i ], m );
157        }

158        sort( sol, sol + d );
159        forint 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