【快速筛可能是大费马素数的数】

为什么写这样一篇博客呢……，原因是今天的多校有道题用到了NTT，然后用双模数然后中国剩余定理我被卡到死，直到我中间运算数组换成了int才卡过。之后我发现，竟然有一个O(1)的大数快速乘法（进行运算的数都在long long内）。根据这个，我发现NTT完全可以通过寻找一个超级大的费马素数来实现（比运算后的系数还大~），代替双模数的中国剩余定理的NTT写法，因为这个常数太大了……

通过我无聊到极点的耐心，我把这个方法写的更加无脑了，鉴于本人低能，以及方便的前提下，只枚举2的幂次，以及另外的两个素数，代码如下，就是纯暴力~

这样随便选一组可行的费马素数，就可以了……前提是这一组确实是素数，233……

my code:my~~code:

[code]#include <bits/stdc++.h>
using namespace std ;

typedef long long LL ;
typedef unsigned long long ULL ;

const int MAXN = 5005 ;

LL mod ;
ULL L_bound = 3e18 , R_bound = 4.5e18 ;//筛取费马素数的上下界
bool prime[MAXN] ;

LL mul ( LL x , LL y ) {//O(1)快速乘法，进行运算的数都在long long内
    return ( x * y - ( long long ) ( x / ( long double ) mod * y + 1e-3 ) * mod + mod ) % mod ;
}

LL power ( LL a , LL b ) {//快速幂
    LL res = 1 , tmp = a ;
    while ( b ) {
        if ( b & 1 ) res = mul ( res , tmp ) ;
        tmp = mul ( tmp , tmp ) ;
        b >>= 1 ;
    }
    return res ;
}

void preprocess () {//预处理素数
    memset ( prime , 1 , sizeof prime ) ;
    prime[0] = prime[1] = 0 ;
    for ( LL i = 2 ; i < MAXN ; ++ i ) if ( prime[i] ) {
        for ( LL j = i * i ; j < MAXN ; j += i ) prime[j] = 0 ;
    }
}

int main () {
    int cas = 0 ;
    preprocess () ;
    for ( LL P0 = 23 ; ( 1LL << P0 ) <= R_bound ; ++ P0 ) if ( prime[P0] ) {
        ULL tmp = ( 1LL << P0 ) ;
        for ( LL P1 = 3 ; tmp * P1 <= R_bound && P1 < MAXN ; ++ P1 ) if ( prime[P1] ) {
            for ( LL P2 = 3 ; tmp * P1 * P2 <= R_bound && P2 < MAXN ; ++ P2 ) if ( P2 > P1 && prime[P2] ) {
                mod = ( 1LL << P0 ) * P1 * P2 + 1 ;
                if ( mod < L_bound ) continue ;
                LL t = mod - 1 ;
                LL x0 = t / 2 ;
                LL x1 = t / P1 ;
                LL x2 = t / P2 ;
                int ok = 0 ;
                for ( LL i = 1 ; i <= 10 ; ++ i ) {
                    if ( power ( i , x0 ) != 1LL && power ( i , x1 ) != 1LL && power ( i , x2 ) != 1LL ) {
                        if ( i != 2 ) {
                            /*
                                原根得到是2的基本可以确定是非素数……
                                而不是2的应该大部分都是素数吧……
                                我猜的……
                                这个可以人工验算，或者去套用Miller_Rabbin（米勒-拉宾）算法判断下就好……
                            */
                            LL tmp = mod , x = 0 ;
                            for ( ; tmp ; ++ x , tmp /= 10 ) ;
                            printf ( "Case #%d:\n" , ++ cas ) ;
                            printf ( "P0 = %I64d , P1 = %I64d , P2 = %I64d\n" , P0 , P1 , P2 ) ;
                            printf ( "mod = %I64d , g = %I64d\n" , mod , i ) ;
                            printf ( "%I64d\n\n" , x ) ;
                            ok = 1 ;
                        }
                        break ;
                    }
                    if ( ok ) break ;
                }
            }
        }
    }
    return 0 ;
}