hdu 3930 Broot 二次剩余

模素数p的原根g的优美体现在每个模p的非零数以g的幂次出现。所以，对任何数1 <= a < p，我们可选择幂

 　　　　　　　　g,g^2,g^2,```````,g^(p-2),g^(p-1)

中恰好一个与a模p同余。相应的指数被称为以g为底的a模p的指标。假设p与g已给定，则记指标为I(a)。

以下以2模13的所有幂的形式：

I　　　　　　　   1　  2　　3　   4　  5　　 6　　 7　　 8　   9　　10　  11　12

2^I(mod 13)　  2　　4　　8　　3　　6　　12　　11　　9　　5　　10　　7　　1

例如，为求I(11)，我们搜寻表的第二行直到找到数11，则指标I(11)=7可从第一行得到。

指标法则

　　a).　　I(ab)=I(a)+I(b) (mod p-1)

　　b).　　I(a^k)=kI(a) (mod p-1)

g^I(ab)=ab=g^I(a)g^I(b)=g^(I(a)+I(b)) (mod p)

故有g^I(x) = x (mod p-1)

也可以这么理解：g^I(value) = id (mod p-1)

例题：3*x^30=4(mod 37)

I(3*x^30)=I(4)

I(3)+30*I(x)=I(4) (mod 36)

26+30*I(x)=2 (mod 36)

30*I(x)=-24=12 (mod 36)

对于这里I(4)=2解释一下：

其实就是g^2 = 4 (mod 36)

其实这个可以根据g^x =b (mod p)来求

对于本题，g就是37的原根2,b=I(x) p =36

也就是求满足2^x = I(x) (mod 36)的x。用Baby Step Giant Step算法即可哦 

提醒：两边不要除以6以得到5*I(x)+2 (mod 36),否则会丢失一些解。

ax = c (mod m)

由扩展欧几里德，知：

I(x)=4,10,16,22,28,34

最后，有指标表（书论书上有哦，自己不放写程序看看），得到x的对应值

I(16)=4,I(25)=10,I(9)=16

 I(21)=22,I(12)=28,I(28)=34

其实这里，也可以根据指标的原式公式g^I(x) = x (mod p-1)

故，同余式3*x^30=4 (mod 37)有6个解，即

x=16,25,9,21,12,28 (mod 37)

/*
 * hdu3930.c
 *
 *  Created on: 2011-10-12
 *      Author: bjfuwangzhu
 */

#include<math.h>
#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#define LL long long
#define nmax 1000010
typedef struct num {
	int ii;
	LL value;
} num;
num Num[nmax];
int flag[nmax], prime[nmax], cpfactor[100];
LL pfactor[100];
int plen, len_pfactor;
LL k, n, p, proot, x, y;
void mkprime() {
	int i, j;
	memset(flag, -1, sizeof(flag));
	for (i = 2, plen = 0; i < nmax; i++) {
		if (flag[i]) {
			prime[plen++] = i;
		}
		for (j = 0; (j < plen) && (i * prime[j] < nmax); j++) {
			flag[i * prime[j]] = 0;
			if (i % prime[j] == 0) {
				break;
			}
		}
	}
}
void findpFactor(LL n) {
	int i, te, cnt;
	te = (int) sqrt(n * 1.0);
	for (i = 0, len_pfactor = 0; (i < plen) && (prime[i] <= te); i++) {
		if (n % prime[i] == 0) {
			cnt = 0;
			while (n % prime[i] == 0) {
				cnt++;
				n /= prime[i];
			}
			pfactor[len_pfactor] = (LL) prime[i];
			cpfactor[len_pfactor++] = cnt;
		}
	}
	if (n > 1) {
		pfactor[len_pfactor] = n;
		cpfactor[len_pfactor++] = 1;
	}
}

LL modular_multi(LL a, LL b, LL c) {
	LL res, temp;
	res = 0, temp = a % c;
	while (b) {
		if (b & 1) {
			res += temp;
			if (res >= c) {
				res -= c;
			}
		}
		temp <<= 1;
		if (temp >= c) {
			temp -= c;
		}
		b >>= 1;
	}
	return res;
}
/*快速幂取余a^b%c*/LL modular_exp(LL a, LL b, LL c) {
	LL res, temp;
	res = 1 % c, temp = a % c;
	while (b) {
		if (b & 1) {
			res = modular_multi(res, temp, c);
		}
		temp = modular_multi(temp, temp, c);
		b >>= 1;
	}
	return res;
}
int dfs(int depth, LL now) {
	int i;
	LL res, temp;
	if (depth == len_pfactor) {
		res = modular_exp(proot, now, p);
		if ((res == 1) && (now != (p - 1))) {
			return 0;
		}
		return 1;
	}
	for (i = 0, temp = 1; i <= cpfactor[depth]; i++) {
		if (!dfs(depth + 1, now * temp)) {
			return 0;
		}
		temp = temp * pfactor[depth];
	}
	return 1;
}
void primitive() {
	findpFactor(p - 1);
	for (proot = 2;; proot++) {
		if (dfs(0, 1)) {
			return;
		}
	}
}
LL extend_gcd(LL a, LL b) {
	LL d, xx;
	if (b == 0) {
		x = 1, y = 0;
		return a;
	}
	d = extend_gcd(b, a % b);
	xx = x;
	x = y, y = xx - a / b * y;
	return d;
}

int bfindNum(LL key, int n) {
	int left, right, mid;
	left = 0, right = n;
	while (left <= right) {
		mid = (left + right) >> 1;
		if (Num[mid].value == key) {
			return Num[mid].ii;
		} else if (Num[mid].value > key) {
			right = mid - 1;
		} else {
			left = mid + 1;
		}
	}
	return -1;
}
int cmp(const void *a, const void *b) {
	num n = *(num *) a;
	num m = *(num *) b;
	LL temp = n.value - m.value;
	if (temp > 0) {
		return 1;
	} else if (temp < 0) {
		return -1;
	}
	return 0;
}
/* a^x = b (mod c)*/LL baby_step_giant_step(LL a, LL b, LL c) {
	int i, j, te;
	LL temp, xx, aa;
	te = (int) (sqrt(c * 1.0) + 0.5);
	for (i = 0, temp = 1 % c; i <= te; i++) {
		Num[i].ii = i;
		Num[i].value = temp;
		temp = temp * a % c;
	}
	aa = Num[te].value;
	qsort(Num, te + 1, sizeof(Num[0]), cmp);
	for (i = 0, temp = 1; i <= te; i++) {
		extend_gcd((int) (temp), c);
		xx = x;
		xx = xx * b;
		xx = xx % c + c;
		x = xx % c;
		j = bfindNum(x, te + 1);
		if (j != -1) {
			return (LL) (i) * te + j;
		}
		temp = temp * aa % c;
	}
	return -1;
}
int rcmp(const void *a, const void *b) {
	LL temp = *(LL *) a - *(LL *) b;
	if (temp > 0) {
		return 1;
	} else if (temp < 0) {
		return -1;
	}
	return 0;

}
/* ax = b (mod c)*/LL result[1001];
void solve(LL a, LL b, LL c) {
	int i;
	LL d;
	d = extend_gcd(a, c);
	if (b % d) {
		puts("-1");
		return;
	}
	b /= d, c /= d;
	result[0] = ((LL) x * b % c + c) % c;
	for (i = 1; i < d; i++) {
		result[i] = result[i - 1] + c;
	}
	for (i = 0; i < d; i++) {
		result[i] = modular_exp(proot, result[i], p);
	}
	qsort(result, d, sizeof(result[0]), rcmp);
	for (i = 0; i < d; i++) {
		printf("%I64d\n", result[i]);
	}
}
int main() {
#ifndef ONLINE_JUDGE
	freopen("data.in", "r", stdin);
#endif
	int cas;
	LL a, b, c;
	mkprime();
	cas = 0;
	/*x^k = n (mod p) */
	while (~scanf("%I64d %I64d %I64d", &k, &p, &n)) {
		primitive();
		b = baby_step_giant_step(proot, n, p);
		a = k, c = p - 1;
		printf("case%d:\n", ++cas);
		solve(a, b, c);
	}
	return 0;
}