foj 1971 A math problem

AC大牛的神题。

比赛时候没写出来。

首先，是思维问题。

一开始大部分队伍想的都是算出N，M然后找公式求，也确实找到一些，但是都不可行。赛后看了解题报告才明白将 ANS 写成

ANS = wi mod (pi ^ ci)

的形式，通过求 wi 。 再套用中国剩余定理。确定 ANS 。

首先， 直接求出 N 。

然后 求 wi = 1^A + 2^A + 3^A + ...+ N ^A mod (pi^ci)

化成

wi = t * (1^A + 2^A + 3^A + 4^A + 5^A +....(pi^ci-1)^A ) + 1^A + 2^A + 3^A +....+bi^A.

其中 t = N / pi^ci.

这个应该好推。

1...bi 部分直接暴力

然后关键变成求前一部分。

也就是 1 .2. 3 .. pi^ci-1 这些数字的A 次方和。

这里，将这些数字划分成两个集合 R.S.

R为 pi^ci 的简化剩余系。 S 为 pi^ci 的非简化剩余系。

由于 A >> cI . 所以集合 S 中的所有元素 比满足 k^A = 0 mod(pi ^ ci)。

重点在与求R中所有元素。

(AC大牛说这里有个结论可以 直接上。但是我硬是没证明出来。求解释。这里是用第二种方法写的)

g为 pi^ci 的原根 (PS: 原根存在条件，2,4,p^l,2p^l,p为奇素数。当pi=2&&ci>2时没原根，要单独处理。)

那么R中的元素可以表示成 g^i. 其中i 遍历 1到 euler(i) 。

那么R中元素的A次方和相当与一个等比数列求和取模。

感觉比赛时候如果开头思路对了应该可以推到最后，但是肯定会被pi=2的情况活活卡死..唉。。就这样吧

#include <stdio.h>
#include <iostream>
#include <math.h>
#include <sstream>
#include <algorithm>
using namespace std;
typedef __int64 u64;
template<class T>inline string toStr(const T& v){ostringstream os;os<<v;return os.str();}
#define FOR(i,s,e) for (int (i)=(s);(i)<(e);i++)
#define DEBUG(a)     printf("%s = %s\n", #a, toStr(a).c_str())
const double eps=1e-8;
template<class T>inline T isqr(T v) {return v*v;}
template<class T>inline int isgn(T v) {return (v>eps) - (v<-eps);}
const double pi=acos(-1.0);
#define FIN(s) freopen(s,"r",stdin);
#define FOUT(s) freopen(s,"w",stdout);
u64 b[30];
u64 p[30];
u64 c[30];
u64 ph[30];
u64 m[30];
bool is[130000];//用来求素数[1,130000]
int prm[30000];//用来保存素数
int totleprm;//记录素数总个数
int getprm(u64 n)
{
     u64 i, j, k = 0;
    int s, e = (int)(sqrt(0.0 + n) + 1);
    memset(is, 1, sizeof(is));
    prm[k++] = 2;
    is[0] = is[1] = 0;
    for (i = 4; i < n; i +=2) is[i] = 0;
    for (i = 3; i < e; i +=2) if(is[i])
    { prm[k++] = i;
        for (s = i * 2, j = i * i ; j < n; j +=s)
        is[j] = 0;
        }
    for (; i < n; i +=2)
    if (is[i])     prm[k++] = i;
    return k;
}

u64 gcd(u64 a,u64 b)
{
    if (a < b)
    return gcd(b,a);
    u64 c;
    while (b)
    {c=a%b;a=b;b=c;}
    return a;
}
u64 product_mod(u64 a,u64 b,u64 c)
{
    u64 ret=0,tmp=a%c;
    while(b)
    {
        if(b&0x1)
            if((ret+=tmp)>=c)
                ret-=c;
            if((tmp<<=1)>=c)
                tmp-=c;
            b>>=1;
    }
    return ret;
}
u64 power_mod(u64 A, u64 B, u64 C)
{
    u64 R = 1, D = A;
    while (B )
    {
        if (B&1) R = product_mod(R, D, C);
        D = product_mod(D, D, C);
        B >>=1;
    }
    return R;
}
u64 power(u64 A, u64 B)
{
    u64 R = 1, D = A;
    while (B )
    {
        if (B&1) R =R*D;
        D = D*D;
        B >>=1;
    }
    return R;
}
u64 extgcd(u64 a,u64 b,u64 &x,u64 &y)
{
    if (b == 0) {x=1;y=0;return a;}
    u64 d = extgcd(b, a%b,x,y);
    u64 t = x; x = y; y = t - a/b * y;
    return d;
}
u64 getPrmRoot(u64 m,u64 phi)
{
    if (m == 2)
    return 1;
    u64 fac[1000];
    u64 n = phi;
    int e = (int)(sqrt(0.0 + n) + 1);
    int length = 0;
    for (int i=0; i<totleprm && prm[i]*prm[i] <= phi; i++)
    if (n % prm[i] == 0)
    {
        fac[length++] = phi/prm[i];
        while ( n % prm[i] == 0) n/=prm[i];
    }

    if (n!=1)
    fac[length++] = phi/n;
    u64 g = 2;
    while (g < m)
    {
        if (gcd(g,m) != 1)
        {g++;continue;}
        bool sign = true;
        FOR(i,0,length)
        if (power_mod(g,fac[i],m) == 1 )
        {
            sign = false;
            break;
        }
        if (sign)
        break;
        g++;
    }
//表示原根不存在
    if (g == m)
    return -1;

    return g;
}


u64 china(int k, u64 a[],u64 m[], u64& lcm)
{
    bool flag = false;
    u64 e, x, y, i,d;
    u64 result;
    u64 a1,m1;
    u64 a2,m2;
    m1 = m[0]; a1 = a[0];
    FOR(i,1,k)
    {
        m2 = m[i]; a2 = a[i];
        d = extgcd( m1, m2, x, y );
        if ( (a2-a1) % d !=   0 )
        flag = 1;
        result = ( x * ((a2-a1) / d ) % m2 + m2 ) % m2;
        a1 = a1 + m1 * result; //对于求多个方程
        m1 = (m1 * m2) / d;     //lcm(m1,m2)最小公倍数；
        a1 = (a1 % m1 + m1) % m1;
    }
    lcm = m1;
    if (flag)
    return -1;
    else
    return a1;
}

u64 getSum(u64 p, u64 n, u64 m)
{
    if (n == 1)
    return (p%m);
    u64 k = getSum(p,(n)/2,m);
    if (n % 2 == 0)
    {
        u64 k_help = power_mod(p,n/2,m);
        return (k+product_mod(k_help,k,m) )%m;
    }
    else
    if (n % 2 == 1)
    {
        u64 k_help = power_mod(p,n/2+1,m);
       return (k+product_mod(k_help,k,m)+k_help )%m;
    }

}
int main()
{
// FIN("in.txt");
   // FOUT("out.txt");
    totleprm = getprm(100000);
    int t ;
    u64 A,k;
    u64 sum;
    u64 sum_help;
    u64 time;
    u64 N,lcm;
    u64 ans;
    u64 g;
    int num=1;
    scanf("%d",&t);
    while (t--)
    {
        scanf("%I64d%I64d",&A,&k);
        k++;
        FOR(i,0,k)
        {
            scanf("%I64d%I64d%I64d",&b[i],&p[i],&c[i]);
            m[i]=power(p[i],c[i]);
            ph[i]=m[i] / p[i]* (p[i]-1);
        }
        N = china(k,b,m,lcm);
        if (N == 0)
        N = lcm;

//暴力求解规模较小的情况可以避免处理一些边缘数据，应该算是数论中的一个小技巧吧
        if (N < 50)
        {
            sum = 0;
            FOR(i,1,N+1)
            sum = (sum + power_mod(i,A,lcm)) % lcm;
            printf("Case %d: %I64d\n",num++,sum);
            continue;
        }
        else
        FOR(i,0,k)
        {
            sum = 0;

            FOR(j,1,b[i]+1)
            sum=(sum + power_mod(j,A,m[i]) ) %m[i];

            time = N / m[i];

//这里是单独处理pi=2的情况，但是我还是想知道那个结论是咋推出来的
            if (p[i] == 2)
            {
                b[i] = (sum + time*ph[i])%m[i];
                continue;
            }

//求原根
            g = getPrmRoot(m[i],ph[i]);
            g = power_mod(g,A,m[i]);

//等比数列求和

             sum_help = getSum(g,ph[i],m[i]);
             sum = (sum+time*sum_help)%m[i];
             b[i] = sum%m[i];
        }
        ans = china(k,b,m,lcm);
        printf("Case %d: %I64d\n",num++,ans);
    }
    return 0;
}