SPOJ MOD Power Modulo Inverted

October 22, 20152015-10-22hahaschoolLeave a comment

题意

求最小的x使得A^x=B(mod C)。

思路

与HDU 2815 Mod Tree、POJ 3243 Clever Y完全一致，extended BSGS的模板题。

代码

#include <iostream>
#include <algorithm>
#include <queue>
#include <vector>
#include <cstdlib>
#include <cstdio>
#include <string>
#include <cstring>
#include <ctime>
#include <iomanip>
#include <cmath>
#include <set>
#include <stack>
#include <cmath>
#include <map>
#include <complex>

using namespace std;

#pragma mark - Extended Baby-Step-Gaint-Step Algorithm (exBSGS)

struct Hashtable{//Mod Based Hash Table
    static const int HASHMOD = 3894229;//A Prime
    long long top, hash[HASHMOD+100], value[HASHMOD+100], stk[1 << 16];
    int locate(long long &x){
        int h = x % HASHMOD;
        while(hash[h] != -1 && hash[h] != x){
            h++;
        }
        return h;
    }
    void insert(long long key,long long val){
        int h = locate(key);
        if(hash[h] == -1){
            hash[h] = key;
            value[h] = val;
            stk[++top] = h;
        }
    }
    long long find(long long &key){
        int h = locate(key);
        return (hash[h] == key)?value[h]:-1;
    }
    void clear(){
        while(top){
            hash[stk[top--]] = -1;
        }
    }
    Hashtable(){
        top = 0;
        memset(hash, -1, sizeof(hash));
    }
} hashtab;

struct Triple{
    long long x,y,z;
    Triple(long long _x, long long _y, long long _z){
        x = _x;
        y = _y;
        z = _z;
    }
    Triple(){
        x = y = z = 0;
    }
};

Triple exgcd(long long a,long long b){
    //for inv use(A/B % C):let a = B,b = C,return.x = inv,return.z = gcd
    //for solve use:return.x = specx,return.y = specy,returnz = gcd
    if (b == 0) {
        return Triple(1,0,a);
    }
    Triple last = exgcd(b, a%b);
    return Triple(last.y, last.x - a / b * last.y, last.z);
}

long long exBSGS(long long A,long long B,long long C){
    //A^x=B(mod C)
    //Ensure B is legal
    B %= C;
    if (B < 0) {
        B += C;
    }

    //Procceed A,B,C to normal BSGS if C is not prime
    long long tmp = 1, cnt = 0, D = 1;
    for (int i = 0; i < 64; i++) {
        if (tmp == B) {
            return i;
        }
        tmp = tmp * A % C;
    }

    for (Triple res; (res = exgcd(A, C)).z != 1; cnt++) {
        if (B % res.z) {
            return -1;
        }
        B /= res.z;
        C /= res.z;
        D = D * A / res.z % C;
    }
    //Normal BSGS
    long long sqrtn = (long long)(ceil(sqrt(C)));
    hashtab.clear();
    long long base = 1;
    for (int i = 0; i < sqrtn; i++) {
        hashtab.insert(base, i);
        base = base * A % C;
    }
    long long j = -1,i = 0;
    for (; i < sqrtn; i++) {
        Triple res = exgcd(D, C);
        long long c = C / res.z;
        res.x = (res.x * B/res.z % c + c) % c;
        j = hashtab.find(res.x);
        if (j != -1) {
            return i * sqrtn + j + cnt;
        }
        D = D * base % C;
    }
    return -1;
}


#pragma mark -

int main(){
    long long a,b,c;
    while (scanf(" %lld %lld %lld",&a,&c,&b) != EOF) {
        if (a == 0 && b == 0 && c == 0) {
            break;
        }
        long long ret = exBSGS(a, b, c);
        if (ret == -1) {
            puts("No Solution");
        }else{
            printf("%lld\n",ret);
        }
    }

    return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

#include <iostream>

#include <algorithm>

#include <queue>

#include <vector>

#include <cstdlib>

#include <cstdio>

#include <string>

#include <cstring>

#include <ctime>

#include <iomanip>

#include <cmath>

#include <set>

#include <stack>

#include <cmath>

#include <map>

#include <complex>

using namespace std;

#pragma mark - Extended Baby-Step-Gaint-Step Algorithm (exBSGS)

struct Hashtable{//Mod Based Hash Table

static const int HASHMOD = 3894229;//A Prime

long long top, hash[HASHMOD+100], value[HASHMOD+100], stk[1 << 16];

int locate(long long &x){

int h = x % HASHMOD;

while(hash[h] != -1 && hash[h] != x){

h++;

}

return h;

}

void insert(long long key,long long val){

int h = locate(key);

if(hash[h] == -1){

hash[h] = key;

value[h] = val;

stk[++top] = h;

}

long long find(long long &key){

int h = locate(key);

return (hash[h] == key)?value[h]:-1;

}

void clear(){

while(top){

hash[stk[top--]] = -1;

}

Hashtable(){

top = 0;

memset(hash, -1, sizeof(hash));

}

} hashtab;

struct Triple{

long long x,y,z;

Triple(long long _x, long long _y, long long _z){

x = _x;

y = _y;

z = _z;

}

Triple(){

x = y = z = 0;

}

};

Triple exgcd(long long a,long long b){

//for inv use(A/B % C):let a = B,b = C,return.x = inv,return.z = gcd

//for solve use:return.x = specx,return.y = specy,returnz = gcd

if (b == 0) {

return Triple(1,0,a);

}

Triple last = exgcd(b, a%b);

return Triple(last.y, last.x - a / b * last.y, last.z);

}

long long exBSGS(long long A,long long B,long long C){

//A^x=B(mod C)

//Ensure B is legal

B %= C;

if (B < 0) {

B += C;

}

//Procceed A,B,C to normal BSGS if C is not prime

long long tmp = 1, cnt = 0, D = 1;

for (int i = 0; i < 64; i++) {

if (tmp == B) {

return i;

}

tmp = tmp * A % C;

}

for (Triple res; (res = exgcd(A, C)).z != 1; cnt++) {

if (B % res.z) {

return -1;

}

B /= res.z;

C /= res.z;

D = D * A / res.z % C;

}

//Normal BSGS

long long sqrtn = (long long)(ceil(sqrt(C)));

hashtab.clear();

long long base = 1;

for (int i = 0; i < sqrtn; i++) {

hashtab.insert(base, i);

base = base * A % C;

}

long long j = -1,i = 0;

for (; i < sqrtn; i++) {

Triple res = exgcd(D, C);

long long c = C / res.z;

res.x = (res.x * B/res.z % c + c) % c;

j = hashtab.find(res.x);

if (j != -1) {

return i * sqrtn + j + cnt;

}

D = D * base % C;

}

return -1;

}

#pragma mark -

int main(){

long long a,b,c;

while (scanf(" %lld %lld %lld",&a,&c,&b) != EOF) {

if (a == 0 && b == 0 && c == 0) {

break;

}

long long ret = exBSGS(a, b, c);

if (ret == -1) {

puts("No Solution");

}else{

printf("%lld\n",ret);

}

return 0;

}

HDU 2815 Mod Tree

October 22, 20152015-10-22hahaschoolOne Comment

题意

题目翻译好之后实际上是求解最小的x使得A^x=B(mod C)成立，C不保证是素数。

思路

典型的Extended BSGS模板题。

需要注意的是，内存空间比较紧张，HASH空间不宜开得太大，需要牺牲一些时间换空间。

代码

#include <iostream>
#include <algorithm>
#include <queue>
#include <vector>
#include <cstdlib>
#include <cstdio>
#include <string>
#include <cstring>
#include <ctime>
#include <iomanip>
#include <cmath>
#include <set>
#include <stack>
#include <cmath>
#include <map>
#include <complex>

using namespace std;

#pragma mark - Extended Baby-Step-Gaint-Step Algorithm (exBSGS)

struct Hashtable{//Mod Based Hash Table
    static const int HASHMOD = 1200611;//A Prime
    long long top, hash[HASHMOD+100], value[HASHMOD+100], stk[1 << 16];
    int locate(long long &x){
        int h = x % HASHMOD;
        while(hash[h] != -1 && hash[h] != x){
            h++;
        }
        return h;
    }
    void insert(long long key,long long val){
        int h = locate(key);
        if(hash[h] == -1){
            hash[h] = key;
            value[h] = val;
            stk[++top] = h;
        }
    }
    long long find(long long &key){
        int h = locate(key);
        return (hash[h] == key)?value[h]:-1;
    }
    void clear(){
        while(top){
            hash[stk[top--]] = -1;
        }
    }
    Hashtable(){
        top = 0;
        memset(hash, -1, sizeof(hash));
    }
} hashtab;

struct Triple{
    long long x,y,z;
    Triple(long long _x, long long _y, long long _z){
        x = _x;
        y = _y;
        z = _z;
    }
    Triple(){
        x = y = z = 0;
    }
};

Triple exgcd(long long a,long long b){
    //for inv use(A/B % C):let a = B,b = C,return.x = inv,return.z = gcd
    //for solve use:return.x = specx,return.y = specy,returnz = gcd
    if (b == 0) {
        return Triple(1,0,a);
    }
    Triple last = exgcd(b, a%b);
    return Triple(last.y, last.x - a / b * last.y, last.z);
}

long long exBSGS(long long A,long long B,long long C){
    //A^x=B(mod C)
    //Ensure B is legal
    B %= C;
    if (B < 0) {
        B += C;
    }

    //Procceed A,B,C to normal BSGS if C is not prime
    long long tmp = 1, cnt = 0, D = 1;
    for (int i = 0; i < 64; i++) {
        if (tmp == B) {
            return i;
        }
        tmp = tmp * A % C;
    }

    for (Triple res; (res = exgcd(A, C)).z != 1; cnt++) {
        if (B % res.z) {
            return -1;
        }
        B /= res.z;
        C /= res.z;
        D = D * A / res.z % C;
    }
    //Normal BSGS
    long long sqrtn = (long long)(ceil(sqrt(C)));
    hashtab.clear();
    long long base = 1;
    for (int i = 0; i < sqrtn; i++) {
        hashtab.insert(base, i);
        base = base * A % C;
    }
    long long j = -1,i = 0;
    for (; i < sqrtn; i++) {
        Triple res = exgcd(D, C);
        long long c = C / res.z;
        res.x = (res.x * B/res.z % c + c) % c;
        j = hashtab.find(res.x);
        if (j != -1) {
            return i * sqrtn + j + cnt;
        }
        D = D * base % C;
    }
    return -1;
}


#pragma mark -

int main(){
    long long a,b,c;
    while (scanf(" %lld %lld %lld",&a,&c,&b) != EOF) {
        long long ret = exBSGS(a, b, c);
        if (b >= c || ret == -1) {
            puts("Orz,I can’t find D!");
        }else{
            printf("%lld\n",ret);
        }
    }

    return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

#include <iostream>

#include <algorithm>

#include <queue>

#include <vector>

#include <cstdlib>

#include <cstdio>

#include <string>

#include <cstring>

#include <ctime>

#include <iomanip>

#include <cmath>

#include <set>

#include <stack>

#include <cmath>

#include <map>

#include <complex>

using namespace std;

#pragma mark - Extended Baby-Step-Gaint-Step Algorithm (exBSGS)

struct Hashtable{//Mod Based Hash Table

static const int HASHMOD = 1200611;//A Prime

long long top, hash[HASHMOD+100], value[HASHMOD+100], stk[1 << 16];

int locate(long long &x){

int h = x % HASHMOD;

while(hash[h] != -1 && hash[h] != x){

h++;

}

return h;

}

void insert(long long key,long long val){

int h = locate(key);

if(hash[h] == -1){

hash[h] = key;

value[h] = val;

stk[++top] = h;

}

long long find(long long &key){

int h = locate(key);

return (hash[h] == key)?value[h]:-1;

}

void clear(){

while(top){

hash[stk[top--]] = -1;

}

Hashtable(){

top = 0;

memset(hash, -1, sizeof(hash));

}

} hashtab;

struct Triple{

long long x,y,z;

Triple(long long _x, long long _y, long long _z){

x = _x;

y = _y;

z = _z;

}

Triple(){

x = y = z = 0;

}

};

Triple exgcd(long long a,long long b){

//for inv use(A/B % C):let a = B,b = C,return.x = inv,return.z = gcd

//for solve use:return.x = specx,return.y = specy,returnz = gcd

if (b == 0) {

return Triple(1,0,a);

}

Triple last = exgcd(b, a%b);

return Triple(last.y, last.x - a / b * last.y, last.z);

}

long long exBSGS(long long A,long long B,long long C){

//A^x=B(mod C)

//Ensure B is legal

B %= C;

if (B < 0) {

B += C;

}

//Procceed A,B,C to normal BSGS if C is not prime

long long tmp = 1, cnt = 0, D = 1;

for (int i = 0; i < 64; i++) {

if (tmp == B) {

return i;

}

tmp = tmp * A % C;

}

for (Triple res; (res = exgcd(A, C)).z != 1; cnt++) {

if (B % res.z) {

return -1;

}

B /= res.z;

C /= res.z;

D = D * A / res.z % C;

}

//Normal BSGS

long long sqrtn = (long long)(ceil(sqrt(C)));

hashtab.clear();

long long base = 1;

for (int i = 0; i < sqrtn; i++) {

hashtab.insert(base, i);

base = base * A % C;

}

long long j = -1,i = 0;

for (; i < sqrtn; i++) {

Triple res = exgcd(D, C);

long long c = C / res.z;

res.x = (res.x * B/res.z % c + c) % c;

j = hashtab.find(res.x);

if (j != -1) {

return i * sqrtn + j + cnt;

}

D = D * base % C;

}

return -1;

}

#pragma mark -

int main(){

long long a,b,c;

while (scanf(" %lld %lld %lld",&a,&c,&b) != EOF) {

long long ret = exBSGS(a, b, c);

if (b >= c || ret == -1) {

puts("Orz,I can’t find D!");

}else{

printf("%lld\n",ret);

}

return 0;

}

POJ 3243 Clever Y

October 22, 20152015-10-22hahaschoolLeave a comment

题意

求解最小的x使得A^x=B(mod C)成立，C不保证是素数。

思路

C不保证素数的BSGS，这时使用扩展BSGS解决，扩展BSGS的实质是在BSGS之前利用这个定理：

Ad=Bd(mod C) <=> A=B(mod C/gcd(C,d))

将问题转化为一般BSGS能够解决的形式，再使用一般BSGS解决。

代码

POJ 2417 Discrete Logging这个问题中也使用了如下exBSGS的代码。

#include <iostream>
#include <algorithm>
#include <queue>
#include <vector>
#include <cstdlib>
#include <cstdio>
#include <string>
#include <cstring>
#include <ctime>
#include <iomanip>
#include <cmath>
#include <set>
#include <stack>
#include <cmath>
#include <map>
#include <complex>

using namespace std;

#pragma mark - Extended Baby-Step-Gaint-Step Algorithm (exBSGS)

struct Hashtable{//Mod Based Hash Table
    static const int HASHMOD = 3894229;//A Prime
    long long top, hash[HASHMOD+100], value[HASHMOD+100], stk[1 << 16];
    int locate(long long &x){
        int h = x % HASHMOD;
        while(hash[h] != -1 && hash[h] != x){
            h++;
        }
        return h;
    }
    void insert(long long key,long long val){
        int h = locate(key);
        if(hash[h] == -1){
            hash[h] = key;
            value[h] = val;
            stk[++top] = h;
        }
    }
    long long find(long long &key){
        int h = locate(key);
        return (hash[h] == key)?value[h]:-1;
    }
    void clear(){
        while(top){
            hash[stk[top--]] = -1;
        }
    }
    Hashtable(){
        top = 0;
        memset(hash, -1, sizeof(hash));
    }
} hashtab;

struct Triple{
    long long x,y,z;
    Triple(long long _x, long long _y, long long _z){
        x = _x;
        y = _y;
        z = _z;
    }
    Triple(){
        x = y = z = 0;
    }
};

Triple exgcd(long long a,long long b){
    //for inv use(A/B % C):let a = B,b = C,return.x = inv,return.z = gcd
    //for solve use:return.x = specx,return.y = specy,returnz = gcd
    if (b == 0) {
        return Triple(1,0,a);
    }
    Triple last = exgcd(b, a%b);
    return Triple(last.y, last.x - a / b * last.y, last.z);
}

long long exBSGS(long long A,long long B,long long C){
    //A^x=B(mod C)
    //Ensure B is legal
    B %= C;
    if (B < 0) {
        B += C;
    }

    //Procceed A,B,C to normal BSGS if C is not prime
    long long tmp = 1, cnt = 0, D = 1;
    for (int i = 0; i < 64; i++) {
        if (tmp == B) {
            return i;
        }
        tmp = tmp * A % C;
    }

    for (Triple res; (res = exgcd(A, C)).z != 1; cnt++) {
        if (B % res.z) {
            return -1;
        }
        B /= res.z;
        C /= res.z;
        D = D * A / res.z % C;
    }
    //Normal BSGS
    long long sqrtn = (long long)(ceil(sqrt(C)));
    hashtab.clear();
    long long base = 1;
    for (int i = 0; i < sqrtn; i++) {
        hashtab.insert(base, i);
        base = base * A % C;
    }
    long long j = -1,i = 0;
    for (; i < sqrtn; i++) {
        Triple res = exgcd(D, C);
        long long c = C / res.z;
        res.x = (res.x * B/res.z % c + c) % c;
        j = hashtab.find(res.x);
        if (j != -1) {
            return i * sqrtn + j + cnt;
        }
        D = D * base % C;
    }
    return -1;
}


#pragma mark -

int main(){
    long long a,b,c;
    while (scanf(" %lld %lld %lld",&a,&c,&b) != EOF) {
        if (a == 0 && b == 0 && c == 0) {
            break;
        }
        long long ret = exBSGS(a, b, c);
        if (ret == -1) {
            puts("No Solution");
        }else{
            printf("%lld\n",ret);
        }
    }

    return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

#include <iostream>

#include <algorithm>

#include <queue>

#include <vector>

#include <cstdlib>

#include <cstdio>

#include <string>

#include <cstring>

#include <ctime>

#include <iomanip>

#include <cmath>

#include <set>

#include <stack>

#include <cmath>

#include <map>

#include <complex>

using namespace std;

#pragma mark - Extended Baby-Step-Gaint-Step Algorithm (exBSGS)

struct Hashtable{//Mod Based Hash Table

static const int HASHMOD = 3894229;//A Prime

long long top, hash[HASHMOD+100], value[HASHMOD+100], stk[1 << 16];

int locate(long long &x){

int h = x % HASHMOD;

while(hash[h] != -1 && hash[h] != x){

h++;

}

return h;

}

void insert(long long key,long long val){

int h = locate(key);

if(hash[h] == -1){

hash[h] = key;

value[h] = val;

stk[++top] = h;

}

long long find(long long &key){

int h = locate(key);

return (hash[h] == key)?value[h]:-1;

}

void clear(){

while(top){

hash[stk[top--]] = -1;

}

Hashtable(){

top = 0;

memset(hash, -1, sizeof(hash));

}

} hashtab;

struct Triple{

long long x,y,z;

Triple(long long _x, long long _y, long long _z){

x = _x;

y = _y;

z = _z;

}

Triple(){

x = y = z = 0;

}

};

Triple exgcd(long long a,long long b){

//for inv use(A/B % C):let a = B,b = C,return.x = inv,return.z = gcd

//for solve use:return.x = specx,return.y = specy,returnz = gcd

if (b == 0) {

return Triple(1,0,a);

}

Triple last = exgcd(b, a%b);

return Triple(last.y, last.x - a / b * last.y, last.z);

}

long long exBSGS(long long A,long long B,long long C){

//A^x=B(mod C)

//Ensure B is legal

B %= C;

if (B < 0) {

B += C;

}

//Procceed A,B,C to normal BSGS if C is not prime

long long tmp = 1, cnt = 0, D = 1;

for (int i = 0; i < 64; i++) {

if (tmp == B) {

return i;

}

tmp = tmp * A % C;

}

for (Triple res; (res = exgcd(A, C)).z != 1; cnt++) {

if (B % res.z) {

return -1;

}

B /= res.z;

C /= res.z;

D = D * A / res.z % C;

}

//Normal BSGS

long long sqrtn = (long long)(ceil(sqrt(C)));

hashtab.clear();

long long base = 1;

for (int i = 0; i < sqrtn; i++) {

hashtab.insert(base, i);

base = base * A % C;

}

long long j = -1,i = 0;

for (; i < sqrtn; i++) {

Triple res = exgcd(D, C);

long long c = C / res.z;

res.x = (res.x * B/res.z % c + c) % c;

j = hashtab.find(res.x);

if (j != -1) {

return i * sqrtn + j + cnt;

}

D = D * base % C;

}

return -1;

}

#pragma mark -

int main(){

long long a,b,c;

while (scanf(" %lld %lld %lld",&a,&c,&b) != EOF) {

if (a == 0 && b == 0 && c == 0) {

break;

}

long long ret = exBSGS(a, b, c);

if (ret == -1) {

puts("No Solution");

}else{

printf("%lld\n",ret);

}

return 0;

}

POJ 2417 Discrete Logging

October 22, 20152015-10-22hahaschoolOne Comment

题意

求解最小的x使得A^x=B(mod C)成立。数据保证了C一定为素数。

思路

解决这种问题是有一种特定算法Baby Step Giant Step的。BSGS算法的具体思想这里不再赘述。直接套用就可以了。

代码

#include <iostream>
#include <algorithm>
#include <queue>
#include <vector>
#include <cstdlib>
#include <cstdio>
#include <string>
#include <cstring>
#include <ctime>
#include <iomanip>
#include <cmath>
#include <set>
#include <stack>
#include <cmath>
#include <map>
#include <complex>

using namespace std;

#pragma mark - Extended Baby-Step-Gaint-Step Algorithm (exBSGS)

struct Hashtable{//Mod Based Hash Table
    static const int HASHMOD = 3894229;//A Prime
    long long top, hash[HASHMOD+100], value[HASHMOD+100], stk[1 << 16];
    int locate(long long x){
        int h = x % HASHMOD;
        while(hash[h] != -1 && hash[h] != x){
            h++;
        }
        return h;
    }
    void insert(long long key,long long val){
        int h = locate(key);
        if(hash[h] == -1){
            hash[h] = key;
            value[h] = val;
            stk[++top] = h;
        }
    }
    long long find(long long key){
        int h = locate(key);
        return (hash[h] == key)?value[h]:-1;
    }
    void clear(){
        while(top){
            hash[stk[top--]] = -1;
        }
    }
    Hashtable(){
        top = 0;
        memset(hash, -1, sizeof(hash));
    }
} hashtab;

struct Triple{
    long long x,y,z;
    Triple(long long _x, long long _y, long long _z){
        x = _x;
        y = _y;
        z = _z;
    }
    Triple(){
        x = y = z = 0;
    }
};

Triple exgcd(long long a,long long b){
    //for inv use(A/B % C):let a = B,b = C,return.x = inv,return.z = gcd
    //for solve use:return.x = specx,return.y = specy,returnz = gcd
    if (b == 0) {
        return Triple(1,0,a);
    }
    Triple last = exgcd(b, a%b);
    return Triple(last.y, last.x - a / b * last.y, last.z);
}

long long exBSGS(long long A,long long B,long long C){
    //A^x=B(mod C)
    //Ensure B is legal
    B %= C;
    if (B < 0) {
        B += C;
    }
    //Procceed A,B,C to normal BSGS if C is not prime
    long long tmp = 1 % C, cnt = 0, D = 1 % C;
    for (int i = 0; i < 64; i++) {
        if (tmp == B) {
            return i;
        }
        tmp = tmp * A % C;
    }

    for (Triple res; (res = exgcd(A, C)).z != 1; cnt++) {
        if (B % res.z) {
            return -1;
        }
        B /= res.z;
        C /= res.z;
        D = D * A / res.z % C;
    }
    //Normal BSGS
    long long sqrtn = (long long)(ceil(sqrt(C)));
    hashtab.clear();
    long long base = 1 % C;
    for (int i = 0; i < sqrtn; i++) {
        hashtab.insert(base, i);
        base = base * A % C;
    }
    long long j = -1,i = 0;
    for (; i < sqrtn; i++) {
        Triple res = exgcd(D, C);
        long long c = C / res.z;
        res.x = (res.x * B/res.z % c + c) % c;
        j = hashtab.find(res.x);
        if (j != -1) {
            return i * sqrtn + j + cnt;
        }
        D = D * base % C;
    }
    return -1;
}


#pragma mark -

int main(){
    long long a,b,c;
    while (scanf(" %lld %lld %lld",&c,&a,&b) != EOF) {
        long long ret = exBSGS(a, b, c);
        if (ret == -1) {
            puts("no solution");
        }else{
            printf("%lld\n",ret);
        }
    }

    return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

#include <iostream>

#include <algorithm>

#include <queue>

#include <vector>

#include <cstdlib>

#include <cstdio>

#include <string>

#include <cstring>

#include <ctime>

#include <iomanip>

#include <cmath>

#include <set>

#include <stack>

#include <cmath>

#include <map>

#include <complex>

using namespace std;

#pragma mark - Extended Baby-Step-Gaint-Step Algorithm (exBSGS)

struct Hashtable{//Mod Based Hash Table

static const int HASHMOD = 3894229;//A Prime

long long top, hash[HASHMOD+100], value[HASHMOD+100], stk[1 << 16];

int locate(long long x){

int h = x % HASHMOD;

while(hash[h] != -1 && hash[h] != x){

h++;

}

return h;

}

void insert(long long key,long long val){

int h = locate(key);

if(hash[h] == -1){

hash[h] = key;

value[h] = val;

stk[++top] = h;

}

long long find(long long key){

int h = locate(key);

return (hash[h] == key)?value[h]:-1;

}

void clear(){

while(top){

hash[stk[top--]] = -1;

}

Hashtable(){

top = 0;

memset(hash, -1, sizeof(hash));

}

} hashtab;

struct Triple{

long long x,y,z;

Triple(long long _x, long long _y, long long _z){

x = _x;

y = _y;

z = _z;

}

Triple(){

x = y = z = 0;

}

};

Triple exgcd(long long a,long long b){

//for inv use(A/B % C):let a = B,b = C,return.x = inv,return.z = gcd

//for solve use:return.x = specx,return.y = specy,returnz = gcd

if (b == 0) {

return Triple(1,0,a);

}

Triple last = exgcd(b, a%b);

return Triple(last.y, last.x - a / b * last.y, last.z);

}

long long exBSGS(long long A,long long B,long long C){

//A^x=B(mod C)

//Ensure B is legal

B %= C;

if (B < 0) {

B += C;

}

//Procceed A,B,C to normal BSGS if C is not prime

long long tmp = 1 % C, cnt = 0, D = 1 % C;

for (int i = 0; i < 64; i++) {

if (tmp == B) {

return i;

}

tmp = tmp * A % C;

}

for (Triple res; (res = exgcd(A, C)).z != 1; cnt++) {

if (B % res.z) {

return -1;

}

B /= res.z;

C /= res.z;

D = D * A / res.z % C;

}

//Normal BSGS

long long sqrtn = (long long)(ceil(sqrt(C)));

hashtab.clear();

long long base = 1 % C;

for (int i = 0; i < sqrtn; i++) {

hashtab.insert(base, i);

base = base * A % C;

}

long long j = -1,i = 0;

for (; i < sqrtn; i++) {

Triple res = exgcd(D, C);

long long c = C / res.z;

res.x = (res.x * B/res.z % c + c) % c;

j = hashtab.find(res.x);

if (j != -1) {

return i * sqrtn + j + cnt;

}

D = D * base % C;

}

return -1;

}

#pragma mark -

int main(){

long long a,b,c;

while (scanf(" %lld %lld %lld",&c,&a,&b) != EOF) {

long long ret = exBSGS(a, b, c);

if (ret == -1) {

puts("no solution");

}else{

printf("%lld\n",ret);

}

return 0;

}

HDU 1402 A * B Problem Plus

October 22, 20152015-10-22hahaschoolLeave a comment

题意

实现两个大整数的加法

思路

整数加法可以看成多项式加法，然后就可以利用FFT通过进行卷积计算的方法来实现。

浮点数误差不会卡掉这个题。

代码

#include <iostream>
#include <algorithm>
#include <queue>
#include <vector>
#include <cstdlib>
#include <cstdio>
#include <string>
#include <cstring>
#include <ctime>
#include <iomanip>
#include <cmath>
#include <set>
#include <stack>
#include <cmath>
#include <map>
#include <complex>

using namespace std;

#pragma mark - Fast Fourier Transformation(FFT()) and Convolution Calaulation(Convolution())
const double PI = acos(-1.0);
typedef complex<double> Complex;

void FFT_Rader(Complex *F,int len){
    //BIT-REVERSE-COPY() 倒位序重新排列数组
    int j = len >> 1;
    for (int i = 1; i < len - 1; i++) {
        if (i < j) {
            swap(F[i], F[j]);
        }
        int k = len >> 1;
        while (j >= k) {
            j -= k;
            k >>= 1;
        }
        if (j < k) {
            j += k;
        }
    }
}

void FFT(Complex *F,int len,int on){
    //on是方向标记，当on为1时，执行FFT，时域->频域；当on为-1时，执行IFFT，频域->时域
    FFT_Rader(F, len);  //把F数组倒位序处理
    for (int h = 2; h <= len; h <<= 1) {    //分治后当前DFT的长度为h(不是使用阶段表示)
        Complex wn(cos(-on * 2 * PI / h),sin(-on * 2 * PI / h));    //使用欧拉定理展开复数单位根
        for (int j = 0; j < len; j += h) {  //向量树上节点对应开始区段
            Complex w(1,0);     //旋转因子
            for (int k = j; k < j + h / 2; k++) {
                Complex u = F[k];
                Complex t = w * F[k + h / 2];
                F[k] = u + t;   //蝶形操作合并
                F[k + h /2] = u - t;    //DFT的对称性
                w = w * wn;     //更新旋转因子
            }
        }
    }
    if (on == -1) {     //IFFT在sigma外面要除以N(DFT长度)
        for (int i = 0; i < len; i++) {
            F[i] /= len;
        }
    }
}

void Convolution(Complex *a,Complex *b,int len){
    //卷积定理:FFT(a卷b)=FFT(a)*FFT(b)
    //a和b的长度必须统合一致，不足的用0补齐
    //长度必须是2的整数次幂
    //将会直接用a向量来接受卷积结果，两个向量都会被破坏
    FFT(a, len, 1);
    FFT(b, len, 1);
    for (int i = 0; i < len; i++) {
        a[i] *= b[i];
    }
    FFT(a, len, -1);
}
#pragma mark -

const int MAXN = 500005;

Complex A[MAXN],B[MAXN];
char sa[MAXN],sb[MAXN];
int res[MAXN];

int prepare(char *a,char *b){
    int lena = (int)strlen(a);
    int lenb = (int)strlen(b);
    int len = 1;
    while (len < 2*lena || len < 2*lenb) {
        len <<= 1;
    }
    for (int i = 0; i < len; i++) {
        A[i].real(0);
        B[i].real(0);
        if (i < lena) {
            A[i].real(sa[lena-i-1]-'0');
        }
        if (i < lenb) {
            B[i].real(sb[lenb-i-1]-'0');
        }

        A[i].imag(0);
        B[i].imag(0);
    }
    return len;
}

void solve(){
    memset(res, 0, sizeof(res));
    int len = prepare(sa, sb);
    Convolution(A, B, len);
    for (int i = 0; i < len; i++) {
        A[i] += 0.5;    //规避浮点误差
        res[i] = A[i].real();
    }
    int bound = 0;
    for (int i = 0; i < len; i++) {
        res[i+1] += res[i] / 10;
        res[i] %= 10;
    }
    for (int i = len - 1; i >= 0; i--) {
        if (res[i]) {
            bound = i;
            break;
        }
    }
    for (int i = bound; i >= 0; i--) {
        printf("%d",res[i]);
    }
    putchar('\n');
}

int main(){
    while (scanf(" %s %s",sa,sb) != EOF) {
        solve();
    }
    return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

#include <iostream>

#include <algorithm>

#include <queue>

#include <vector>

#include <cstdlib>

#include <cstdio>

#include <string>

#include <cstring>

#include <ctime>

#include <iomanip>

#include <cmath>

#include <set>

#include <stack>

#include <cmath>

#include <map>

#include <complex>

using namespace std;

#pragma mark - Fast Fourier Transformation(FFT()) and Convolution Calaulation(Convolution())

const double PI = acos(-1.0);

typedef complex<double> Complex;

void FFT_Rader(Complex *F,int len){

//BIT-REVERSE-COPY() 倒位序重新排列数组

int j = len >> 1;

for (int i = 1; i < len - 1; i++) {

if (i < j) {

swap(F[i], F[j]);

}

int k = len >> 1;

while (j >= k) {

j -= k;

k >>= 1;

}

if (j < k) {

j += k;

}

void FFT(Complex *F,int len,int on){

//on是方向标记，当on为1时，执行FFT，时域->频域；当on为-1时，执行IFFT，频域->时域

FFT_Rader(F, len); //把F数组倒位序处理

for (int h = 2; h <= len; h <<= 1) { //分治后当前DFT的长度为h(不是使用阶段表示)

Complex wn(cos(-on * 2 * PI / h),sin(-on * 2 * PI / h)); //使用欧拉定理展开复数单位根

for (int j = 0; j < len; j += h) { //向量树上节点对应开始区段

Complex w(1,0); //旋转因子

for (int k = j; k < j + h / 2; k++) {

Complex u = F[k];

Complex t = w * F[k + h / 2];

F[k] = u + t; //蝶形操作合并

F[k + h /2] = u - t; //DFT的对称性

w = w * wn; //更新旋转因子

}

if (on == -1) { //IFFT在sigma外面要除以N(DFT长度)

for (int i = 0; i < len; i++) {

F[i] /= len;

}

void Convolution(Complex *a,Complex *b,int len){

//卷积定理:FFT(a卷b)=FFT(a)*FFT(b)

//a和b的长度必须统合一致，不足的用0补齐

//长度必须是2的整数次幂

//将会直接用a向量来接受卷积结果，两个向量都会被破坏

FFT(a, len, 1);

FFT(b, len, 1);

for (int i = 0; i < len; i++) {

a[i] *= b[i];

}

FFT(a, len, -1);

}

#pragma mark -

const int MAXN = 500005;

Complex A[MAXN],B[MAXN];

char sa[MAXN],sb[MAXN];

int res[MAXN];

int prepare(char *a,char *b){

int lena = (int)strlen(a);

int lenb = (int)strlen(b);

int len = 1;

while (len < 2*lena || len < 2*lenb) {

len <<= 1;

}

for (int i = 0; i < len; i++) {

A[i].real(0);

B[i].real(0);

if (i < lena) {

A[i].real(sa[lena-i-1]-'0');

}

if (i < lenb) {

B[i].real(sb[lenb-i-1]-'0');

}

A[i].imag(0);

B[i].imag(0);

}

return len;

}

void solve(){

memset(res, 0, sizeof(res));

int len = prepare(sa, sb);

Convolution(A, B, len);

for (int i = 0; i < len; i++) {

A[i] += 0.5; //规避浮点误差

res[i] = A[i].real();

}

int bound = 0;

for (int i = 0; i < len; i++) {

res[i+1] += res[i] / 10;

res[i] %= 10;

}

for (int i = len - 1; i >= 0; i--) {

if (res[i]) {

bound = i;

break;

}

for (int i = bound; i >= 0; i--) {

printf("%d",res[i]);

}

putchar('\n');

}

int main(){

while (scanf(" %s %s",sa,sb) != EOF) {

solve();

}

return 0;

}

HDU5446

正确的拆乘防爆，超大素数取模，Lucas定理，CRT合并。

#include <iostream>
#include <algorithm>
#include <queue>
#include <vector>
#include <cstdlib>
#include <cstdio>
#include <string>
#include <cstring>
#include <ctime>
#include <iomanip>
#include <cmath>
#include <set>
#include <stack>
#include <cmath>
#include <map>

using namespace std;


const int maxn = 100005;
long long mod = 1;
long long factor[12];

long long A[12];
long long fac[12][maxn];

inline long long mul(long long a,long long b,long long moder){
    long long ret = 0;
    if (a > b) {
        swap(a, b);
    }
    long long fac = a;
    while (b) {
        if (b & 1) {
            ret += fac;
            if (ret >= moder) {
                ret -= moder;
            }
        }
        fac += fac;
        if (fac >= moder) {
            fac -= moder;
        }
        b>>= 1;
    }
    return ret;
}


long long pow_mod(long long n, long long m, long long p) {
    n %= p;
    long long ret = 1;
    for(; m; m>>=1, n = mul(n, n, p))
        if(m & 1) ret = mul(ret, n, p);
    return ret;
}

long long inv(long long n, long long p) {
    return pow_mod(n, p-2, p);
}

long long C(long long n, long long m, int i,long long p) {
    if(n < m) return 0;
    return mul(fac[i][n], inv(mul(fac[i][m], fac[i][n-m], p), factor[i]),p);
}

long long lucas(long long n, long long m, int i) {
    if(m == 0) return 1;
    long long p = factor[i];
    return  mul(lucas(n/p, m/p, i), C(n%p, m%p, i,p), p);
}

long long CRT(int k) {
    long long m = mod, x = 0;
    for(int i = 0; i < k; i++) {
        long long w = m / factor[i];
        x += mul(mul(inv(w, factor[i]), A[i], m), w, m);
        x %= m;
    }
    return (x + m) % m;
}

int main()
{
    int caseCnt = 0;
    scanf(" %d",&caseCnt);
    for (int d = 1; d <= caseCnt; d++) {
        long long n,m,k;
        memset(fac, 0, sizeof(fac));
        mod = 1;
        scanf(" %lld %lld %lld",&n,&m,&k);
        long long maxfac = 0;
        for (int i = 0; i < k; i++) {
            cin >> factor[i];
            maxfac = max(factor[i],maxfac);
            mod *= factor[i];
        }
        for(int i = 0; i < k; i++) {
            fac[i][0] = 1;
            for(int j = 1; j <= maxfac; j++)
                fac[i][j] =   mul(fac[i][j-1], j, factor[i]);
        }
        memset(A, 0, sizeof(A));
        for(int i = 0; i < k; i++) {
            A[i] = (A[i] + lucas(n, m, i)) % factor[i];
        }
        printf("%lld\n",CRT(k));
    }
    return 0;
}

100

101

102

103

104

105

106

107

108

109

110

#include <iostream>

#include <algorithm>

#include <queue>

#include <vector>

#include <cstdlib>

#include <cstdio>

#include <string>

#include <cstring>

#include <ctime>

#include <iomanip>

#include <cmath>

#include <set>

#include <stack>

#include <cmath>

#include <map>

using namespace std;

const int maxn = 100005;

long long mod = 1;

long long factor[12];

long long A[12];

long long fac[12][maxn];

inline long long mul(long long a,long long b,long long moder){

long long ret = 0;

if (a > b) {

swap(a, b);

}

long long fac = a;

while (b) {

if (b & 1) {

ret += fac;

if (ret >= moder) {

ret -= moder;

}

fac += fac;

if (fac >= moder) {

fac -= moder;

}

b>>= 1;

}

return ret;

}

long long pow_mod(long long n, long long m, long long p) {

n %= p;

long long ret = 1;

for(; m; m>>=1, n = mul(n, n, p))

if(m & 1) ret = mul(ret, n, p);

return ret;

}

long long inv(long long n, long long p) {

return pow_mod(n, p-2, p);

}

long long C(long long n, long long m, int i,long long p) {

if(n < m) return 0;

return mul(fac[i][n], inv(mul(fac[i][m], fac[i][n-m], p), factor[i]),p);

}

long long lucas(long long n, long long m, int i) {

if(m == 0) return 1;

long long p = factor[i];

return mul(lucas(n/p, m/p, i), C(n%p, m%p, i,p), p);

}

long long CRT(int k) {

long long m = mod, x = 0;

for(int i = 0; i < k; i++) {

long long w = m / factor[i];

x += mul(mul(inv(w, factor[i]), A[i], m), w, m);

x %= m;

}

return (x + m) % m;

}

int main()

{

int caseCnt = 0;

scanf(" %d",&caseCnt);

for (int d = 1; d <= caseCnt; d++) {

long long n,m,k;

memset(fac, 0, sizeof(fac));

mod = 1;

scanf(" %lld %lld %lld",&n,&m,&k);

long long maxfac = 0;

for (int i = 0; i < k; i++) {

cin >> factor[i];

maxfac = max(factor[i],maxfac);

mod *= factor[i];

}

for(int i = 0; i < k; i++) {

fac[i][0] = 1;

for(int j = 1; j <= maxfac; j++)

fac[i][j] = mul(fac[i][j-1], j, factor[i]);

}

memset(A, 0, sizeof(A));

for(int i = 0; i < k; i++) {

A[i] = (A[i] + lucas(n, m, i)) % factor[i];

}

printf("%lld\n",CRT(k));

}

return 0;

}

(Extended) Baby Step Giant Step Algorithm

用于求出满足A^x=B(mod C)的最小的x。

参数配置：
- HASHMOD：Hashtable的Hash空间大小，爆内存要改小
调用方法：exBSGS(A,B,C);

#pragma mark - Extended Baby-Step-Gaint-Step Algorithm (exBSGS)

struct Hashtable{//Mod Based Hash Table
    static const int HASHMOD = 1200611;//A Prime
    long long top, hash[HASHMOD+100], value[HASHMOD+100], stk[1 << 16];
    int locate(const long long x){
        int h = x % HASHMOD;
        while(hash[h] != -1 && hash[h] != x){
            h++;
        }
        return h;
    }
    void insert(const long long key,const long long val){
        int h = locate(key);
        if(hash[h] == -1){
            hash[h] = key;
            value[h] = val;
            stk[++top] = h;
        }
    }
    long long find(const long long key){
        int h = locate(key);
        return (hash[h] == key)?value[h]:-1;
    }
    void clear(){
        while(top){
            hash[stk[top--]] = -1;
        }
    }
    Hashtable(){
        top = 0;
        memset(hash, -1, sizeof(hash));
    }
} hashtab;

struct Triple{
    long long x,y,z;
    Triple(const long long _x,const long long _y,const long long _z){
        x = _x;
        y = _y;
        z = _z;
    }
    Triple(){
        x = y = z = 0;
    }
};

Triple exgcd(const long long a,const long long b){
    //for inv use(A/B % C):let a = B,b = C,return.x = inv,return.z = gcd
    //for solve use:return.x = specx,return.y = specy,returnz = gcd
    if (b == 0) {
        return Triple(1,0,a);
    }
    Triple last = exgcd(b, a%b);
    return Triple(last.y, last.x - a / b * last.y, last.z);
}

long long exBSGS(long long A,long long B,long long C){
    //A^x=B(mod C)
    //Ensure B is legal
    B %= C;
    if (B < 0) {
        B += C;
    }
    //Procceed A,B,C to normal BSGS if C is not prime
    long long tmp = 1 % C, cnt = 0, D = 1 % C;
    for (int i = 0; i < 64; i++) {
        if (tmp == B) {
            return i;
        }
        tmp = tmp * A % C;
    }

    for (Triple res; (res = exgcd(A, C)).z != 1; cnt++) {
        if (B % res.z) {
            return -1;
        }
        B /= res.z;
        C /= res.z;
        D = D * A / res.z % C;
    }
    //Normal BSGS
    long long sqrtn = (long long)(ceil(sqrt(C)));
    hashtab.clear();
    long long base = 1 % C;
    for (int i = 0; i < sqrtn; i++) {
        hashtab.insert(base, i);
        base = base * A % C;
    }
    long long j = -1,i = 0;
    for (; i < sqrtn; i++) {
        Triple res = exgcd(D, C);
        long long c = C / res.z;
        res.x = (res.x * B/res.z % c + c) % c;
        j = hashtab.find(res.x);
        if (j != -1) {
            return i * sqrtn + j + cnt;
        }
        D = D * base % C;
    }
    return -1;
}


#pragma mark -

100

101

102

103

104

105

#pragma mark - Extended Baby-Step-Gaint-Step Algorithm (exBSGS)

struct Hashtable{//Mod Based Hash Table

static const int HASHMOD = 1200611;//A Prime

long long top, hash[HASHMOD+100], value[HASHMOD+100], stk[1 << 16];

int locate(const long long x){

int h = x % HASHMOD;

while(hash[h] != -1 && hash[h] != x){

h++;

}

return h;

}

void insert(const long long key,const long long val){

int h = locate(key);

if(hash[h] == -1){

hash[h] = key;

value[h] = val;

stk[++top] = h;

}

long long find(const long long key){

int h = locate(key);

return (hash[h] == key)?value[h]:-1;

}

void clear(){

while(top){

hash[stk[top--]] = -1;

}

Hashtable(){

top = 0;

memset(hash, -1, sizeof(hash));

}

} hashtab;

struct Triple{

long long x,y,z;

Triple(const long long _x,const long long _y,const long long _z){

x = _x;

y = _y;

z = _z;

}

Triple(){

x = y = z = 0;

}

};

Triple exgcd(const long long a,const long long b){

//for inv use(A/B % C):let a = B,b = C,return.x = inv,return.z = gcd

//for solve use:return.x = specx,return.y = specy,returnz = gcd

if (b == 0) {

return Triple(1,0,a);

}

Triple last = exgcd(b, a%b);

return Triple(last.y, last.x - a / b * last.y, last.z);

}

long long exBSGS(long long A,long long B,long long C){

//A^x=B(mod C)

//Ensure B is legal

B %= C;

if (B < 0) {

B += C;

}

//Procceed A,B,C to normal BSGS if C is not prime

long long tmp = 1 % C, cnt = 0, D = 1 % C;

for (int i = 0; i < 64; i++) {

if (tmp == B) {

return i;

}

tmp = tmp * A % C;

}

for (Triple res; (res = exgcd(A, C)).z != 1; cnt++) {

if (B % res.z) {

return -1;

}

B /= res.z;

C /= res.z;

D = D * A / res.z % C;

}

//Normal BSGS

long long sqrtn = (long long)(ceil(sqrt(C)));

hashtab.clear();

long long base = 1 % C;

for (int i = 0; i < sqrtn; i++) {

hashtab.insert(base, i);

base = base * A % C;

}

long long j = -1,i = 0;

for (; i < sqrtn; i++) {

Triple res = exgcd(D, C);

long long c = C / res.z;

res.x = (res.x * B/res.z % c + c) % c;

j = hashtab.find(res.x);

if (j != -1) {

return i * sqrtn + j + cnt;

}

D = D * base % C;

}

return -1;

}

#pragma mark -

指数循环节

若Ａ^start = B (mod P)求出最小的len使得Ａ^(k*len) = B (mod P)，原理是这种对于一个值的循环节长度一定是Phi(P)的约数，所以不断约掉因子并尝试。

long long fastpow(long long a,long long p,long long m){
    long long ret = 1;
    long long fac = a;
    while(p){
        if(p&1){
            ret *= fac;
            ret %= m;
        }
        p >>= 1;
        fac *= fac;
        fac %= m;
    }
    return ret;
}

long long euler_phi(long long x){
    long long ret = 1;
    for(long long i = 2; i * i <= x; i++){
        if(x % i == 0){
            x /= i;
            ret *= i-1;
            while(x % i == 0){
                x /= i;
                ret *= i;
            }
        }
    }
    if(x > 1){
        ret *= x-1;
    }
    return ret;
}
long long loop_len(long long start,long long A,long long D,long long P){
    long long phi = euler_phi(P),ret = phi;
    for(long long i = 2; i * i <= phi; i++){
        while(phi % i == 0){
            phi /= i;
        }
        while(ret % i == 0 && fastpow(A,start + ret/i,P) == D){
            ret /= i;
        }
    }
    if (phi > 1) {
        while(ret % phi == 0 && fastpow(A,start + ret/phi,P) == D){
            ret /= phi;
        }
    }
    return ret;

}

long long fastpow(long long a,long long p,long long m){

long long ret = 1;

long long fac = a;

while(p){

if(p&1){

ret *= fac;

ret %= m;

}

p >>= 1;

fac *= fac;

fac %= m;

}

return ret;

}

long long euler_phi(long long x){

long long ret = 1;

for(long long i = 2; i * i <= x; i++){

if(x % i == 0){

x /= i;

ret *= i-1;

while(x % i == 0){

x /= i;

ret *= i;

}

if(x > 1){

ret *= x-1;

}

return ret;

}

long long loop_len(long long start,long long A,long long D,long long P){

long long phi = euler_phi(P),ret = phi;

for(long long i = 2; i * i <= phi; i++){

while(phi % i == 0){

phi /= i;

}

while(ret % i == 0 && fastpow(A,start + ret/i,P) == D){

ret /= i;

}

if (phi > 1) {

while(ret % phi == 0 && fastpow(A,start + ret/phi,P) == D){

ret /= phi;

}

return ret;

}

字符串算法小总结

October 22, 20152015-10-22hahaschoolLeave a comment

Extended Knuth-Morris-Pratt Algorithm

模板

参数调整：MAXN是字符串的理论最大长度
调用：exKMP_getextend(char* 母串,char* 模式串)，会同时处理好next数组和extend数组，数据间无需初始化。
exKMP_next[i]表示模式串mode[i~(len-1)]与整个模式串的最长公共前缀长度
exKMP_extend[i]表示母串str[i~(len-1)]于整个模式串的最长公共前缀长度

#pragma mark - Extended Knuth-Morris-Pratt Algorithm (exKMP)
const int MAXN = 100010;
int exKMP_next[MAXN],exKMP_extend[MAXN];

void exKMP_getnext(char* str){
    //Preproccess mode string
    //next array indicates the public prefix of str[0~n-1] and str[i~n-1]
    //相当于自己和自己匹配的exKMP
    int len = strlen(str);
    exKMP_next[0] = len;
    //get next[1] specially
    int pos = 0;
    while(pos+1 < len && str[pos] == str[pos+1]){
        pos++;
    }
    exKMP_next[1] = pos;
    pos = 1;
    for(int i = 2; i < len; i++){
        if(exKMP_next[i - pos] + i < exKMP_next[pos] + pos){
            //situ. 1
            exKMP_next[i] = exKMP_next[i - pos];
        }else{
            //situ. 2
            int j = exKMP_next[pos] + pos - i;
            if(j < 0){
                j = 0;
            }
            while(i + j < len && str[i+j] == str[j]){
                j++;
            }
            exKMP_next[i] = j;
            pos = i;
        }
    }
}

void exKMP_getextend(char *str,char *mode){
    int len = strlen(str),lenmode = strlen(mode);
    exKMP_getnext(mode);
    int pos = 0;
    while(pos < len && pos < lenmode && str[pos] == mode[pos]){
        pos++;
    }
    exKMP_extend[0] = pos;
    pos = 0;//pos是最近一次成功匹配的开始位置
    for(int i = 1; i < len; i++){
        if(exKMP_next[i - pos] + i < exKMP_extend[pos] + pos){
            //situ.1 模式串的等价前置状态中，等价区没有超过当前匹配范围
            exKMP_extend[i] = exKMP_next[i - pos];
        }else{
            //situ.2 模式串的等价前置状态中等价去超出了当前已经确认的范围
            int j = exKMP_extend[pos] + pos - i;
            if(j < 0){
                j = 0;
            }
            while(i + j < len && j < lenmode && str[j+i] == mode[j]){
                j++;
            }
            exKMP_extend[i] = j;
            pos = i;
        }
    }
}

#pragma mark -

#pragma mark - Extended Knuth-Morris-Pratt Algorithm (exKMP)

const int MAXN = 100010;

int exKMP_next[MAXN],exKMP_extend[MAXN];

void exKMP_getnext(char* str){

//Preproccess mode string

//next array indicates the public prefix of str[0~n-1] and str[i~n-1]

//相当于自己和自己匹配的exKMP

int len = strlen(str);

exKMP_next[0] = len;

//get next[1] specially

int pos = 0;

while(pos+1 < len && str[pos] == str[pos+1]){

pos++;

}

exKMP_next[1] = pos;

pos = 1;

for(int i = 2; i < len; i++){

if(exKMP_next[i - pos] + i < exKMP_next[pos] + pos){

//situ. 1

exKMP_next[i] = exKMP_next[i - pos];

}else{

//situ. 2

int j = exKMP_next[pos] + pos - i;

if(j < 0){

j = 0;

}

while(i + j < len && str[i+j] == str[j]){

j++;

}

exKMP_next[i] = j;

pos = i;

}

void exKMP_getextend(char *str,char *mode){

int len = strlen(str),lenmode = strlen(mode);

exKMP_getnext(mode);

int pos = 0;

while(pos < len && pos < lenmode && str[pos] == mode[pos]){

pos++;

}

exKMP_extend[0] = pos;

pos = 0;//pos是最近一次成功匹配的开始位置

for(int i = 1; i < len; i++){

if(exKMP_next[i - pos] + i < exKMP_extend[pos] + pos){

//situ.1 模式串的等价前置状态中，等价区没有超过当前匹配范围

exKMP_extend[i] = exKMP_next[i - pos];

}else{

//situ.2 模式串的等价前置状态中等价去超出了当前已经确认的范围

int j = exKMP_extend[pos] + pos - i;

if(j < 0){

j = 0;

}

while(i + j < len && j < lenmode && str[j+i] == mode[j]){

j++;

}

exKMP_extend[i] = j;

pos = i;

}

#pragma mark -

题目

HDU 4300 Clairewd’s message

Aho-Corasick Automaton

模板

参数调整：
- SP：字符集大小
- MAXN：理论上Tire树的节点上限
- mark()：用来将字符映射到SP规定的字符集大小内的函数
- ACA_match()：需要按照题目要求更改匹配时行为
调用：
- 初始化：init()
- 插入模式串：insert()
- 构建fail指针：ACA_buildfail()
- 匹配母串：ACA_match()
活用注记：
- root表示没有匹配到任何字符，如果超出字符集，也要转移到root
- 每一个节点的fail指针对应节点是等价的更早状态（包括root状态）
- 每个节点指向失配的出边实际上指向的是其fail指针对应状态中对应出边指向的状态

#pragma mark - AC Automaton
const int SP = 26;
const int MAXN = 600005;

struct Trie{
    int nxt[MAXN][SP];
    int fail[MAXN];
    int freq[MAXN];
    int last[MAXN];//优化，last指向的是当前状态的前缀里面能够配上字符串的状态，而不是所有前缀
    int root,tot;
    inline int mark(char a){
        return a - 'a';
    }
    void init(){
        tot = 0;
        root = newnode();
    }
    int newnode(){
        for(int i = 0; i < SP; i++){
            nxt[tot][i] = -1;
        }
        freq[tot] = 0;
        fail[tot] = -1;
        return tot++;
    }
    void insert(char* str){
        int len = (int)strlen(str);
        int now = root;
        for(int i = 0; i < len; i++){
            if(nxt[now][mark(str[i])] == -1){
                nxt[now][mark(str[i])] = newnode();
            }
            now = nxt[now][mark(str[i])];
        }
        freq[now]++;
    }
    void ACA_buildfail(){
        queue<int> que;
        fail[root] = root;
        last[root] = root;
        for(int i = 0; i < SP; i++){
            if(nxt[root][i] == -1){
                nxt[root][i] = root;
            }else{
                fail[nxt[root][i]] = root;
                last[nxt[root][i]] = root;
                que.push(nxt[root][i]);
            }
        }
        while(!que.empty()){
            int now = que.front();
            que.pop();
            for(int i = 0; i < SP; i++){
                if(nxt[now][i] == -1){
                    nxt[now][i] = nxt[fail[now]][i];
                }else{
                    fail[nxt[now][i]] = nxt[fail[now]][i];
                    last[nxt[now][i]] = freq[nxt[fail[now]][i]]?nxt[fail[now]][i]:last[nxt[fail[now]][i]];
                    que.push(nxt[now][i]);
                }
            }
        }
    }
    int ACA_match(char* str){
        int len = (int)strlen(str);
        int now = root;
        int res = 0;
        for(int i = 0; i < len; i++){
            now = nxt[now][mark(str[i])];
            int temp = now;
            while(temp != root){
                //Proccess matching here
                //res += freq[temp];
                //freq[temp] = 0;//防止重复,but this will damage aca,try alternative
                //temp = last[temp];
            }
        }
        return res;
    }
};
#pragma mark -

#pragma mark - AC Automaton

const int SP = 26;

const int MAXN = 600005;

struct Trie{

int nxt[MAXN][SP];

int fail[MAXN];

int freq[MAXN];

int last[MAXN];//优化，last指向的是当前状态的前缀里面能够配上字符串的状态，而不是所有前缀

int root,tot;

inline int mark(char a){

return a - 'a';

}

void init(){

tot = 0;

root = newnode();

}

int newnode(){

for(int i = 0; i < SP; i++){

nxt[tot][i] = -1;

}

freq[tot] = 0;

fail[tot] = -1;

return tot++;

}

void insert(char* str){

int len = (int)strlen(str);

int now = root;

for(int i = 0; i < len; i++){

if(nxt[now][mark(str[i])] == -1){

nxt[now][mark(str[i])] = newnode();

}

now = nxt[now][mark(str[i])];

}

freq[now]++;

}

void ACA_buildfail(){

queue<int> que;

fail[root] = root;

last[root] = root;

for(int i = 0; i < SP; i++){

if(nxt[root][i] == -1){

nxt[root][i] = root;

}else{

fail[nxt[root][i]] = root;

last[nxt[root][i]] = root;

que.push(nxt[root][i]);

}

while(!que.empty()){

int now = que.front();

que.pop();

for(int i = 0; i < SP; i++){

if(nxt[now][i] == -1){

nxt[now][i] = nxt[fail[now]][i];

}else{

fail[nxt[now][i]] = nxt[fail[now]][i];

last[nxt[now][i]] = freq[nxt[fail[now]][i]]?nxt[fail[now]][i]:last[nxt[fail[now]][i]];

que.push(nxt[now][i]);

}

int ACA_match(char* str){

int len = (int)strlen(str);

int now = root;

int res = 0;

for(int i = 0; i < len; i++){

now = nxt[now][mark(str[i])];

int temp = now;

while(temp != root){

//Proccess matching here

//res += freq[temp];

//freq[temp] = 0;//防止重复,but this will damage aca,try alternative

//temp = last[temp];

}

return res;

}

};

#pragma mark -

题目

[Pending]

Suffix Automaton

模板１～指针实现，动态分配内存

struct SAM_State{
    int mxn;
    SAM_State *to[26];
    SAM_State *parent;
    SAM_State(int _mxn){
        mxn = _mxn;
        memset(to, 0, sizeof(to));
        parent = NULL;
    }
};

SAM_State* SAM_Append(SAM_State *last,SAM_State *root,int ch){
    SAM_State *p = last;
    SAM_State *np = new SAM_State(p -> mxn + 1);
    while (p && p -> to[ch] == NULL) {
        p -> to[ch] = np;
        p = p -> parent;
    }
    if (p != NULL) {
        SAM_State *q = p -> to[ch];
        if (q -> mxn == p -> mxn + 1) {
            np -> parent = q;
        }else{
            SAM_State *nq = new SAM_State(p -> mxn + 1);
            nq -> parent = q -> parent;
            q -> parent = nq;
            np -> parent = nq;
            memcpy(nq -> to, q -> to, sizeof(q -> to));
            while (p && p -> to[ch] == q) {
                p -> to[ch] = nq;
                p = p -> parent;
            }
        }
    }else{
        np -> parent = root;
    }
    last = np;
    return last;
}

SAM_State* SAM_Build(string str){
    SAM_State *root = new SAM_State(0);
    SAM_State *last = root;
    for (int i = 0; i < str.size(); i++) {
        last = SAM_Append(last, root, str[i] - 'a');
    }
    return root;
}

struct SAM_State{

int mxn;

SAM_State *to[26];

SAM_State *parent;

SAM_State(int _mxn){

mxn = _mxn;

memset(to, 0, sizeof(to));

parent = NULL;

}

};

SAM_State* SAM_Append(SAM_State *last,SAM_State *root,int ch){

SAM_State *p = last;

SAM_State *np = new SAM_State(p -> mxn + 1);

while (p && p -> to[ch] == NULL) {

p -> to[ch] = np;

p = p -> parent;

}

if (p != NULL) {

SAM_State *q = p -> to[ch];

if (q -> mxn == p -> mxn + 1) {

np -> parent = q;

}else{

SAM_State *nq = new SAM_State(p -> mxn + 1);

nq -> parent = q -> parent;

q -> parent = nq;

np -> parent = nq;

memcpy(nq -> to, q -> to, sizeof(q -> to));

while (p && p -> to[ch] == q) {

p -> to[ch] = nq;

p = p -> parent;

}

}else{

np -> parent = root;

}

last = np;

return last;

}

SAM_State* SAM_Build(string str){

SAM_State *root = new SAM_State(0);

SAM_State *last = root;

for (int i = 0; i < str.size(); i++) {

last = SAM_Append(last, root, str[i] - 'a');

}

return root;

}

模板２～数组实现，支持删除后缀

const int MAXN = 510005;

struct SAM_Node{
    int mxn;
    int to[26];
    int parent;

    int del;
    int pos;
    int endflg;

    void init(){
        mxn = 0;
        parent = -1;
        del = 0;
        pos = -1;
        endflg = 0;
        memset(to,-1,sizeof(to));
    }
} pool[MAXN*2];

int stempos[MAXN];
bool delflag[MAXN];

int SAM_tot;
int del_tot;
int dfs_tot;

void SAM_Init(){
    pool[0].init();
    SAM_tot = 1;
    del_tot = 1;
    dfs_tot = 0;
    delflag[0] = false;
}

int SAM_Append(int last,int ch){
    int p = last;
    int np = SAM_tot++;
    pool[np].init();
    pool[np].mxn = pool[p].mxn + 1;

    pool[np].del = del_tot++;
    delflag[pool[np].del] = false;
    pool[np].pos = pool[np].mxn - 1;
    stempos[pool[np].pos] = np;

    while(p != -1 && (pool[p].to[ch] == -1 || delflag[pool[pool[p].to[ch]].del])){
        pool[p].to[ch] = np;
        p = pool[p].parent;
    }
    if(p != -1){
        int q = pool[p].to[ch];
        if(pool[q].mxn == pool[p].mxn + 1){
            pool[np].parent = q;
        }else{
            int nq = SAM_tot++;
            pool[nq] = pool[q];
            pool[nq].mxn = pool[p].mxn + 1;
            pool[q].parent = nq;
            pool[np].parent = nq;
            while(p != -1 && pool[p].to[ch] == q){
                pool[p].to[ch] = nq;
                p = pool[p].parent;
            }
        }
    }else{
        pool[np].parent = 0;
    }
    return np;
}

int len,last,qlen;
char str[MAXN];

void suffix_delete(int siz){
    int i,j;
    for(i = len-1,j = 1; j <= siz; i--,j++){
        delflag[pool[stempos[i]].del] = true;
    }
    len -= siz;
    last = stempos[i];
}

void suffix_append(){
    int apdlen = strlen(str);
    for(int i = 0; i < apdlen; i++){
        last = SAM_Append(last,str[i]-'a');
    }
    len += apdlen;
}

const int MAXN = 510005;

struct SAM_Node{

int mxn;

int to[26];

int parent;

int del;

int pos;

int endflg;

void init(){

mxn = 0;

parent = -1;

del = 0;

pos = -1;

endflg = 0;

memset(to,-1,sizeof(to));

}

} pool[MAXN*2];

int stempos[MAXN];

bool delflag[MAXN];

int SAM_tot;

int del_tot;

int dfs_tot;

void SAM_Init(){

pool[0].init();

SAM_tot = 1;

del_tot = 1;

dfs_tot = 0;

delflag[0] = false;

}

int SAM_Append(int last,int ch){

int p = last;

int np = SAM_tot++;

pool[np].init();

pool[np].mxn = pool[p].mxn + 1;

pool[np].del = del_tot++;

delflag[pool[np].del] = false;

pool[np].pos = pool[np].mxn - 1;

stempos[pool[np].pos] = np;

while(p != -1 && (pool[p].to[ch] == -1 || delflag[pool[pool[p].to[ch]].del])){

pool[p].to[ch] = np;

p = pool[p].parent;

}

if(p != -1){

int q = pool[p].to[ch];

if(pool[q].mxn == pool[p].mxn + 1){

pool[np].parent = q;

}else{

int nq = SAM_tot++;

pool[nq] = pool[q];

pool[nq].mxn = pool[p].mxn + 1;

pool[q].parent = nq;

pool[np].parent = nq;

while(p != -1 && pool[p].to[ch] == q){

pool[p].to[ch] = nq;

p = pool[p].parent;

}

}else{

pool[np].parent = 0;

}

return np;

}

int len,last,qlen;

char str[MAXN];

void suffix_delete(int siz){

int i,j;

for(i = len-1,j = 1; j <= siz; i--,j++){

delflag[pool[stempos[i]].del] = true;

}

len -= siz;

last = stempos[i];

}

void suffix_append(){

int apdlen = strlen(str);

for(int i = 0; i < apdlen; i++){

last = SAM_Append(last,str[i]-'a');

}

len += apdlen;

}

欧拉图和哈密顿图小总结

October 18, 20152018-10-16hahaschoolOne Comment

欧拉图

概念

定义

欧拉回路：图中经过每条边一次且仅一次的环；
欧拉路径：图中经过每条边一次且仅一次的路径；
欧拉图：有至少一个欧拉回路的图；
半欧拉图：没有欧拉环，但有至少一条欧拉路径的图。

无向图判定条件

图是联通的（度数为0的点不用考虑，因为是孤立的点，没影响）～先决条件
所有点的度数都是偶数～存在欧拉回路，欧拉图
有且仅有两个点的度数是奇数～存在欧拉路径，半欧拉图
两个奇数度数的点一个是路径起点一个是路径终点

证明：因为任意一个点，欧拉环（或欧拉路径）从它这里进去多少次就要出来多少次，故(进去的次数+出来的次数)为偶数，又因为(进去的次数+出来的次数)=该点的度数（根据定义），所以该点的度数为偶数。

有向图判定条件

把所有的边强行改为无向边之后（这个图叫原图的基图），联通～先决条件
所有点出度等于入度～存在欧拉回路，欧拉图
有且仅有两个点出度不等于入度，一个出度比入度大1，另一个入度比出度大1～欧拉路径，半欧拉图
出度比入度大1的点是路径起点，入度比出度大1的点是路径终点

证明：与无向图证明类似，一个点进去多少次就要出来多少次。

模板

struct Graph{
    int head[MAXN];
    int nxt[MAXE];
    int ind[MAXN];
    int outd[MAXN];
    int from[MAXE];
    int end[MAXE];
    int vis[MAXE];
    int q;
    int V,E;
    inline void addEdge(int _from,int _to){
        from[q] = _from;
        end[q] = _to;
        outd[_from]++;
        ind[_to]++;
        nxt[q] = head[_from];
        head[_from] = q++;
    }
    inline void clear(){
        q = 1;
        memset(head, 0, sizeof(head));
        memset(nxt, 0, sizeof(nxt));
        memset(end, 0, sizeof(end));
        memset(vis, false, sizeof(vis));
        memset(from, 0, sizeof(from));
        memset(ind, 0, sizeof(ind));
        memset(outd, 0, sizeof(outd));
    }
};

stack<int> res;
void euler_dfs(Graph &g,int cur){
    for (int i = g.head[cur]; i; i = g.nxt[i]) {
        if (!g.vis[i]) {
            g.vis[i] = true;
            //g.vis[i^1] = true; 无向图的时候，取消这句注释
            euler_dfs(g,g.end[i]);
            res.push(i);
        }
    }
}

vector<int> singular;
int euler_judge(Graph &g){
    singular.clear();
    for (int i = 1; i <= g.V; i++) {
        if (g.ind[i] != g.outd[i]) {
            singular.push_back(i);
        }
    }
    if (singular.size() == 2) {
        if(g.ind[singular[0]] == g.outd[singular[0]] + 1){
            swap(singular[0],singular[1]);
        }
        if(g.ind[singular[0]] + 1 == g.outd[singular[0]] && g.outd[singular[1]] + 1 == g.ind[singular[1]]){
            return 1;
        }
        return 0;//euler path
        //判断必须严格，差一点也不行，跑出合适个数的结果都可能是错的！
    }else if(singular.size() == 0){
        for (int i = 1; i<= g.V; i++) {
            if (g.ind[i]) {
                singular.push_back(i);
            }
        }
        return 2;//euler loop
    }else{
        return 0;//failed
    }
}

//singular[0]是起点

struct Graph{

int head[MAXN];

int nxt[MAXE];

int ind[MAXN];

int outd[MAXN];

int from[MAXE];

int end[MAXE];

int vis[MAXE];

int q;

int V,E;

inline void addEdge(int _from,int _to){

from[q] = _from;

end[q] = _to;

outd[_from]++;

ind[_to]++;

nxt[q] = head[_from];

head[_from] = q++;

}

inline void clear(){

q = 1;

memset(head, 0, sizeof(head));

memset(nxt, 0, sizeof(nxt));

memset(end, 0, sizeof(end));

memset(vis, false, sizeof(vis));

memset(from, 0, sizeof(from));

memset(ind, 0, sizeof(ind));

memset(outd, 0, sizeof(outd));

}

};

stack<int> res;

void euler_dfs(Graph &g,int cur){

for (int i = g.head[cur]; i; i = g.nxt[i]) {

if (!g.vis[i]) {

g.vis[i] = true;

//g.vis[i^1] = true; 无向图的时候，取消这句注释

euler_dfs(g,g.end[i]);

res.push(i);

}

vector<int> singular;

int euler_judge(Graph &g){

singular.clear();

for (int i = 1; i <= g.V; i++) {

if (g.ind[i] != g.outd[i]) {

singular.push_back(i);

}

if (singular.size() == 2) {

if(g.ind[singular[0]] == g.outd[singular[0]] + 1){

swap(singular[0],singular[1]);

}

if(g.ind[singular[0]] + 1 == g.outd[singular[0]] && g.outd[singular[1]] + 1 == g.ind[singular[1]]){

return 1;

}

return 0;//euler path

//判断必须严格，差一点也不行，跑出合适个数的结果都可能是错的！

}else if(singular.size() == 0){

for (int i = 1; i<= g.V; i++) {

if (g.ind[i]) {

singular.push_back(i);

}

return 2;//euler loop

}else{

return 0;//failed

}

//singular[0]是起点

ACM Course 2 中文配套题面

October 08, 20152015-10-08hahaschoolLeave a comment

Problem B（Original：POJ 3624）音乐游戏

PROBLEM DESCRIPTION

hahaschool非常喜欢玩音乐游戏，虽然他的水平差到不行，但他依然非常享受。众所周知，音乐游戏有好多种，不过作为音游狗，hahaschool显然不会去玩粗制滥造的音乐游戏（比如节操大师），也不会去摸某些打着音乐游戏旗号的音乐播放器（比如deemo）。因为hahaschool比较胖，腿部的力量和灵活性有限，所以也无法玩主要依靠腿部游玩的游戏（比如Dance Dance Revolution），同时因为hahaschool是一个害羞的男孩子，所以不会去玩要求模仿游戏里的动作的游戏（比如Dance Evolution）。即使面临诸多限制，音乐游戏千千万，在音乐游戏的世界里hahaschool总是能找到偷税(划掉)愉悦的好办法，不过因为hahaschool有选择困难症，又是一个理性派，所以每当他走进机厅，总会陷入不知道玩什么的迷茫，浪费宝贵的时间。

欢乐的时光总是短暂的，为了更好地享受音乐游戏，刚刚到达广州市区某机厅的hahaschool拜托你帮他作出选择。假设这家机厅里面有N种音乐游戏（因为机厅非常有钱，所以游戏的种类可能会多达3402种），对于第i种游戏有愉悦程度Di(1<=Di<=100)，第i种游戏要投币Wi个(1<=Wi<=400)，对于每种游戏，hahaschool只会玩一盘，因为这样才能享受到不同的愉悦，玩完一盘之后hahaschool就能得到游戏对应的愉悦程度。现在hahaschool的手里有M枚游戏币，因为不明原因，M可能多达12880，hahaschool并不需要花完全部的游戏币，请你算出hahaschool最多能得到多少愉悦？

INPUT

第一行：两个空格分开的整数N和M

第二行一直到第N+1行：对于第i行，会有两个有空格隔开的数Wi和Di

N,M,Wi,Di的含义及约定范围在题面中已经给出

请注意，测试数据有多组，所有的数据都会放在一个文件里，请处理到文件结束。

OUTPUT

请输出一个整数，是hahaschool在题目的约束条件之下最高能获得的愉悦程度总和，独占一行。

SAMPLE INPUT

4 6
1 4
2 6
3 12
2 7

SAMPLE OUTPUT

Problem C（Original：HDU 1114）壕的魔法财布

PROBLEM DESCRIPTION

众所周知gy非常壕，那么他那么多钱都是从哪里来的呢？

根据《刑法》，我们有：

第三百九十五条第一款国家工作人员的财产、支出明显超过合法收入，差额巨大的，可以责令该国家工作人员说明来源，不能说明来源的，差额部分以非法所得论，处五年以下有期徒刑或者拘役；差额特别巨大的，处五年以上十年以下有期徒刑。财产的差额部分予以追缴。

第九十三条本法所称国家工作人员，是指国家机关中从事公务的人员。国有公司、企业、事业单位、人民团体中从事公务的人员和国家机关、国有公司、企业、事业单位委派到非国有公司、企业、事业单位、社会团体从事公务的人员，以及其他依照法律从事公务的人员，以国家工作人员论。

显然这两条并没有什么卵用，既然我们无法通过巨额财产来源不明罪威胁gy说出他如此壕的原因，我们只能雇佣专业的情报人员。

很快，我们就得到了情报，gy如此壕的秘密在于它有一个宝具“魔法财布”，每当gy需要钱的时候，他只要打开财布，里面就会掉出来金币银币各种币，然后财布会暂时性的变空。

因为gy是壕，所以跟gy跟级花约会的时候一定会买非常奢华的礼物，然而财布掉出来的钱币未必能够支付这件东西的价格，这样就会让gy显得非常没面子，所以gy必须要事先估算出财布里面最少装了多少钱，选择支付得起的最奢华的礼物，这样才能突显他的高贵。

INPUT

在评测过程中，所有的数据都是在一个文件里，请处理到文件结束。

在输入文件中包括了T组数据，在文件的开头，给出了T，独占一行。

之后，对于每组数据的第一行，是两个用空格分隔开来的整数E和F。E表示一个空的财布的重量。F表示gy拿出财布时，他感受到的财布的重量。重量的单位是克，最多10000克，毕竟gy有神力。(1<=E<=F<=10000)

在每组数据的第二行，给出了一个整数N，独占一行。表示财布可能会生成N种钱币.(1<=N<=500)

接下来是N行，每一行会有两个整数P和W，用空格分隔。P表示当前行所表示钱币的价值，W表示当前行所表示钱币的重量。(1<=P<=50000,1<=W<=10000,重量以克给出)。

OUTPUT

请输出财布中所含钱币的最小价值，独占一行。请注意，如果通过给出的钱币无法通过某种方式使财布的重量正好等于gy感受到的财布的重量，请输出-1.

SAMPLE INPUT

3
10 110
2
1 1
30 50
10 110
2
1 1
50 30
1 6
2
10 3
20 4

SAMPLE OUTPUT

60
100
-1