2014.12.29 - intege.diff.hpp

        // Extended Euclidean Algorithm.
        struct eea_result{
            DoubleInt result, s, t;
        };

        static eea_result eea(DoubleInt a, DoubleInt b){
            eea_result p;
            if(b == 0){
                p.result = a;
                p.s = 1;
                p.t = 0;
                return p;
            }
            p = eea(b, std::abs(a % b));
            eea_result q;
            q.result =  p.result;
            q.s = p.t;
            q.t = p.s - (a / b) * p.t;
            return q;
        }

        // FFT.
        static data_type fft(const data_type &arg_data, int coe){
            // arg_data の内部表現データのサイズを 2 の冪乗に繰り上げて得る.
            std::size_t n = ceil_pow2_single(static_cast<UInt>(arg_data.size()));

            // log_2(n) を得る.
            std::size_t lg = bit_num(n) - 1;

            // arg_data 内部表現を並べ替えて a にコピーする.
            std::vector<DoubleInt> a(n), b(n, 1);
            for(std::size_t i = 0; i < arg_data.size(); ++i){
                a[bit_reverse(i, lg)] = static_cast<DoubleInt>(arg_data[i]);
            }

            // n ごとに素数 p = k * 2^n + 1 の最小の k を示すテーブル.
            // 実用上は 2^64 以下までで十分.
            // enumeration_prime_k 関数によって計算された.
            static const UInt k_table[] = {
                1,      1,      2,      1,      3,      3,      2,      1,
                15,     12,     6,      3,      5,      4,      2,      1,
                6,      3,      11,     7,      33,     25,     20,     10,
                5,      7,      15,     12,     6,      3,      35,     18,
                9,      12,     6,      3,      20,     10,     5,      6,
                3,      9,      9,      15,     35,     19,     27,     15,
                14,     7,      14,     7,      20,     10,     5,      27,
                29,     54,     27,     31
            };

            for(std::size_t i = 1; i <= lg; ++i){
                std::size_t m = 1 << i;
                UInt q = k_table[i - 1], p = q * static_cast<UInt>(m) + 1;
                DoubleInt zeta = eea(coe * static_cast<DoubleInt>(q), p).s, omega = 1;
                for(std::size_t j = 0; j < m / 2; ++j){
                    for(std::size_t k = j; k < n; k += m){
                        std::size_t l = k + m / 2;
                        DoubleInt t = mod(omega * b[l], p), u = b[k];
                        b[k] = mod(t + u, p);
                        if(b[k] == 0){
                            a[k] = a[l] + a[k];
                        }
                        b[l] = mod(t - u, p);
                        if(b[l] == 0){
                            a[l] = a[l] - a[k];
                        }
                    }
                    omega = mod(omega * zeta, p);
                }
            }

            data_type result;
            result.resize(a.size());
            for(std::size_t i = 0; i < a.size(); ++i){
                result[i] = static_cast<UInt>(a[i]);
            }
            return result;
        }

        // Bit Reverse.
        static std::size_t bit_reverse(UInt k, UInt n){
            std::size_t l = 0u;
            for(UInt i = 0u; i < n; ++i){
                l += ((k & (1 << i)) ? 1 : 0) << (n - i - 1);
            }
            return l;
        }

        // prime = k * n + 1 になる k を探す.
        static int search_prime_k(const integer &n){
            int k = 1;
            while(!miller_rabin_test(k * n + 1, 100)){
                k += 1;
            }
            return k;
        }

        // prime = k * n + 1 になる k を列挙する.
        template<class OStream>
        static void enumeration_prime_k(OStream &os){
            integer bound = 1;
            for(int i = 0; i < 64; ++i){
                bound <<= 1;
            }
            integer p;
            for(int a = 1; ; ++a){
                p = 1;
                p <<= a;
                if(p > bound){
                    break;
                }
                int k = search_prime_k(p);
                os << k << (a % 8 == 0 ? ",\n" : ",\t");
                if(k * p + 1 > bound){
                    break;
                }
            }
        }

        // Primality Test (2).
        static bool miller_rabin_test(const integer &n, std::size_t s){
            static std::mt19937 g(0xABABABABu);
            integer a, m = n - 2;
            for(std::size_t i = 0; i < s; ++i){
                if(witness_test(m.random(g) + 1, n)){
                    return false;
                }
            }
            return true;
        }

        // Primality Test (1).
        static bool witness_test(const integer &a, const integer &n){
            integer b = n - 1, d = 1, x, t;
            for(int i = static_cast<int>(b.bit_num()) - 1; i >= 0; --i){
                x = d;
                d = (d * d) % n;
                if(d == 1 && x != 1 && x != b){
                    return true;
                }
                t = b >> static_cast<UInt>(i);
                if(t > 0 && t.data[0] % 2 == 1){
                    d = (d * a) % n;
                }
            }
            return d != 1;
        }

        // Random Number
        template<class Generator>
        integer random(Generator &g) const{
            std::size_t s = g() % data.size() + 1;
            integer x;
            x.data.resize(s);
            for(std::size_t i = 0; i < s - 1; ++i){
                x.data[i] = g() & base_mask;
            }
            if(s == data.size()){
                x.data.back() = g() % data.back();
            }
            x.normalize_data_size();
            return x;
        }

        static DoubleInt mod(DoubleInt x, DoubleInt y){
            x = x % y;
            return x >= 0 ? x : (y + x);
        }