Language: C++

#ifndef __INTMOD_H__
#define __INTMOD_H__

/* updated: 2020-06-24 */

#include <iostream>
#include <cassert>
#include <vector>
#include <cstdint>
#include <functional>
#include <unordered_map>

template <uint32_t Mod>
class IntMod {
	static_assert(Mod != 0, "Mod must not be 0");
	static_assert(Mod < 0x80000000, "Mod must be less than 0x80000000");

private:
	uint32_t value_m;

public:
	IntMod();
	IntMod(int32_t value);
	IntMod(uint32_t value);
	IntMod(int64_t value);
	IntMod(uint64_t value);

	const IntMod<Mod>& operator+() const;
	IntMod<Mod> operator-() const;
	IntMod<Mod>& operator++();
	IntMod<Mod>& operator--();
	IntMod<Mod> operator++(int dummy);
	IntMod<Mod> operator--(int dummy);
	IntMod<Mod>& operator+=(const IntMod<Mod>& right);
	IntMod<Mod>& operator-=(const IntMod<Mod>& right);
	IntMod<Mod>& operator*=(const IntMod<Mod>& right);
	IntMod<Mod>& operator/=(const IntMod<Mod>& right);
	IntMod<Mod> operator[](uint64_t exp) const;

	static IntMod<Mod> add(const IntMod<Mod>& a, const IntMod<Mod>& b);
	static IntMod<Mod> subtract(const IntMod<Mod>& a, const IntMod<Mod>& b);
	static IntMod<Mod> multiply(const IntMod<Mod>& a, const IntMod<Mod>& b);
	static IntMod<Mod> divide(const IntMod<Mod>& a, const IntMod<Mod>& b);
	static int compare(const IntMod<Mod>& a, const IntMod<Mod>& b);

	IntMod<Mod> inverse() const;
	IntMod<Mod> power(uint64_t exp) const;
	uint32_t get() const;

	static IntMod<Mod> fact(uint32_t num);
	static IntMod<Mod> invFact(uint32_t num);
	static IntMod<Mod> combi(uint32_t n, uint32_t r);
	static IntMod<Mod> homcombi(uint32_t n, uint32_t r);
	static IntMod<Mod> perm(uint32_t n, uint32_t r);
	static std::vector<IntMod<Mod>> calcInverseList(int max);
};

template <uint32_t Mod>
IntMod<Mod> operator+(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
IntMod<Mod> operator-(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
IntMod<Mod> operator*(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
IntMod<Mod> operator/(const IntMod<Mod>& a, const IntMod<Mod>& b);

template <uint32_t Mod, class Int>
IntMod<Mod> operator+(const IntMod<Mod>& a, const Int& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator-(const IntMod<Mod>& a, const Int& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator*(const IntMod<Mod>& a, const Int& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator/(const IntMod<Mod>& a, const Int& b);

template <uint32_t Mod, class Int>
IntMod<Mod> operator+(const Int& a, const IntMod<Mod>& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator-(const Int& a, const IntMod<Mod>& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator*(const Int& a, const IntMod<Mod>& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator/(const Int& a, const IntMod<Mod>& b);

template <uint32_t Mod>
bool operator==(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
bool operator!=(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
bool operator<(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
bool operator<=(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
bool operator>(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
bool operator>=(const IntMod<Mod>& a, const IntMod<Mod>& b);

template <uint32_t Mod, class Int>
bool operator==(const IntMod<Mod>& a, const Int& b);
template <uint32_t Mod, class Int>
bool operator!=(const IntMod<Mod>& a, const Int& b);

template <uint32_t Mod>
std::istream& operator>>(std::istream& ist, IntMod<Mod>& object);
template <uint32_t Mod>
std::ostream& operator<<(std::ostream& ost, const IntMod<Mod>& object);

using MInt = IntMod<1000000007>;
using MF = IntMod<998244353>;

// for unordered_set/map
namespace std {
	template <uint32_t Mod>
	struct hash<IntMod<Mod>> {
		size_t operator() (const IntMod<Mod>& object) const {
			return std::hash<uint32_t>()(object.get());
		}
	};
}

/*------------------------------------*/
/*-------   Implementations   --------*/
/*------------------------------------*/

template <uint32_t Mod>
IntMod<Mod>::IntMod() 
	: value_m(0) {
}

template <uint32_t Mod>
IntMod<Mod>::IntMod(int32_t value) {
	int32_t tmp = value % int32_t(Mod);
	value_m = tmp >= 0 ? tmp : Mod - uint32_t(-tmp);
}

template <uint32_t Mod>
IntMod<Mod>::IntMod(uint32_t value)
	: value_m(value % Mod) {
}

template <uint32_t Mod>
IntMod<Mod>::IntMod(int64_t value) {
	int32_t tmp = value % int64_t(Mod);
	value_m = tmp >= 0 ? tmp : Mod - uint32_t(-tmp);
}

template <uint32_t Mod>
IntMod<Mod>::IntMod(uint64_t value)
	: value_m(value % Mod) {
}

template <uint32_t Mod>
const IntMod<Mod>& IntMod<Mod>::operator+() const {
	return *this;
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::operator-() const {
	IntMod<Mod> ret;
	if (value_m == 0) {
		ret.value_m = 0;
	} else {
		ret.value_m = Mod - value_m;
	}
	return ret;
}

template <uint32_t Mod>
IntMod<Mod>& IntMod<Mod>::operator++() {
	++value_m;
	if (value_m == Mod) {
		value_m = 0;
	}
	return *this;
}

template <uint32_t Mod>
IntMod<Mod>& IntMod<Mod>::operator--() {
	if (value_m == 0) {
		value_m = Mod;
	}
	--value_m;
	return *this;
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::operator++(int dummy) {
	IntMod ret(*this);
	++*this;
	return ret;
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::operator--(int dummy) {
	IntMod ret(*this);
	--*this;
	return ret;
}

template <uint32_t Mod>
IntMod<Mod>& IntMod<Mod>::operator+=(const IntMod<Mod>& right) {
	value_m += right.value_m;		// note: value_m < 0x80000000
	if (value_m >= Mod) {
		value_m -= Mod;
	}
	return *this;
}

template <uint32_t Mod>
IntMod<Mod>& IntMod<Mod>::operator-=(const IntMod<Mod>& right) {
	if (value_m < right.value_m) {
		value_m += Mod;				// note: value_m < 0x80000000
	}
	value_m -= right.value_m;
	return *this;
}

template <uint32_t Mod>
IntMod<Mod>& IntMod<Mod>::operator*=(const IntMod<Mod>& right) {
	value_m = (uint64_t(value_m) * right.value_m) % Mod;
	return *this;
}

template <uint32_t Mod>
IntMod<Mod>& IntMod<Mod>::operator/=(const IntMod<Mod>& right) {
	*this *= right.inverse();
	return *this;
}

// for power
template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::operator[](uint64_t exp) const {
	return power(exp);
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::add(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	IntMod<Mod> ret(a);
	ret += b;
	return ret;
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::subtract(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	IntMod<Mod> ret(a);
	ret -= b;
	return ret;
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::multiply(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	IntMod<Mod> ret(a);
	ret *= b;
	return ret;
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::divide(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	IntMod<Mod> ret(a);
	ret /= b;
	return ret;
}

template <uint32_t Mod>
int IntMod<Mod>::compare(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	if (a.value_m < b.value_m) return -1;
	if (a.value_m > b.value_m) return 1;
	return 0;
}

/* Mod must be a prime. */
template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::inverse() const {
	assert(value_m != 0);
	return power(Mod - 2);
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::power(uint64_t exp) const {
	IntMod<Mod> product(1);
	IntMod<Mod> factor(*this);
	while (exp > 0) {
		if (exp & 1) {
			product *= factor;
		}
		factor *= factor;
		exp >>= 1;
	}
	return product;
}

template <uint32_t Mod>
uint32_t IntMod<Mod>::get() const {
	return value_m;
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::fact(uint32_t num) {
	assert(num <= 0x80000000);

	static std::vector<IntMod<Mod>> table(1, 1);
	if (table.size() > num) return table[num];

	uint32_t old_size = table.size();
	table.resize(num + 1);
	for (uint32_t i = old_size; i <= num; ++i) {
		table[i] = table[i - 1] * i;
	}
	return table[num];
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::invFact(uint32_t num) {
	assert(num <= 0x80000000);

	static std::vector<IntMod<Mod>> table(1, 1);
	if (table.size() > num) return table[num];

	uint32_t old_size = table.size();
	table.resize(num + 1);
	table[num] = fact(num).inverse();
	for (uint32_t i = num; i >= old_size; --i) {
		table[i - 1] = table[i] * i;
	}
	return table[num];
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::combi(uint32_t n, uint32_t r) {
	assert(n >= r);
	return fact(n) * invFact(n - r) * invFact(r);
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::homcombi(uint32_t n, uint32_t r) {
	assert(n != 0 || r == 0);

	if (n == 0 && r == 0) return IntMod<Mod>(1);
	return combi(n + r - 1, r);
}

template <uint32_t Mod>
IntMod<Mod> IntMod<Mod>::perm(uint32_t n, uint32_t r) {
	assert(n >= r);
	return fact(n) * invFact(n - r);
}

template <uint32_t Mod>
std::vector<IntMod<Mod>> IntMod<Mod>::calcInverseList(int max) {
	assert(max < Mod);

	std::vector<IntMod<Mod>> table(max + 1);
	table[1] = 1;
	for (int i = 2; i <= max; ++i) {
		table[i] = (uint64_t(Mod - Mod / i) * table[Mod % i].get()) % Mod;
	}
	return table;
}

template <uint32_t Mod>
IntMod<Mod> operator+(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::add(a, b);
}

template <uint32_t Mod>
IntMod<Mod> operator-(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::subtract(a, b);
}

template <uint32_t Mod>
IntMod<Mod> operator*(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::multiply(a, b);
}

template <uint32_t Mod>
IntMod<Mod> operator/(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::divide(a, b);
}

template <uint32_t Mod, class Int>
IntMod<Mod> operator+(const IntMod<Mod>& a, const Int& b) {
	return a + IntMod<Mod>(b);
}

template <uint32_t Mod, class Int>
IntMod<Mod> operator-(const IntMod<Mod>& a, const Int& b) {
	return a - IntMod<Mod>(b);
}

template <uint32_t Mod, class Int>
IntMod<Mod> operator*(const IntMod<Mod>& a, const Int& b) {
	return a * IntMod<Mod>(b);
}

template <uint32_t Mod, class Int>
IntMod<Mod> operator/(const IntMod<Mod>& a, const Int& b) {
	return a / IntMod<Mod>(b);
}

template <uint32_t Mod, class Int>
IntMod<Mod> operator+(const Int& a, const IntMod<Mod>& b) {
	return IntMod<Mod>(a) + b;
}

template <uint32_t Mod, class Int>
IntMod<Mod> operator-(const Int& a, const IntMod<Mod>& b) {
	return IntMod<Mod>(a) - b;
}

template <uint32_t Mod, class Int>
IntMod<Mod> operator*(const Int& a, const IntMod<Mod>& b) {
	return IntMod<Mod>(a) * b;
}

template <uint32_t Mod, class Int>
IntMod<Mod> operator/(const Int& a, const IntMod<Mod>& b) {
	return IntMod<Mod>(a) / b;
}

template <uint32_t Mod>
bool operator==(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) == 0;
}

template <uint32_t Mod>
bool operator!=(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) != 0;
}

template <uint32_t Mod>
bool operator<(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) < 0;
}

template <uint32_t Mod>
bool operator<=(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) <= 0;
}

template <uint32_t Mod>
bool operator>(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) > 0;
}

template <uint32_t Mod>
bool operator>=(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) >= 0;
}

template <uint32_t Mod, class Int>
bool operator==(const IntMod<Mod>& a, const Int& b) {
	return a == IntMod<Mod>(b);
}

template <uint32_t Mod, class Int>
bool operator!=(const IntMod<Mod>& a, const Int& b) {
	return a != IntMod<Mod>(b);
}

template <uint32_t Mod>
std::istream& operator>>(std::istream& ist, IntMod<Mod>& object) {
	uint64_t value;
	ist >> value;
	object = IntMod<Mod>(value);
	return ist;
}

template <uint32_t Mod>
std::ostream& operator<<(std::ostream& ost, const IntMod<Mod>& object) {
	ost << object.get();
	return ost;
}

#endif
#ifndef __SETS_H__
#define __SETS_H__

/* updated: 2020-03-28 */

// 集合系の関数などを置きたい

#include <cstdint>
#include <vector>
#include <cassert>

#ifdef _MSC_VER
#include <intrin.h>
#endif

// ビット演算系

/* 立っているビットの数を数える */
inline int popcount(uint64_t value) {
#ifdef _MSC_VER
	return __popcnt64(value);
#elif defined(__GNUC__)
	return __builtin_popcountll(value);
#else
	value = (value & 0x5555555555555555) + ((value >>  1) & 0x5555555555555555);
	value = (value & 0x3333333333333333) + ((value >>  2) & 0x3333333333333333);
	value = (value & 0x0f0f0f0f0f0f0f0f) + ((value >>  4) & 0x0f0f0f0f0f0f0f0f);
	value = (value & 0x00ff00ff00ff00ff) + ((value >>  8) & 0x00ff00ff00ff00ff);
	value = (value & 0x0000ffff0000ffff) + ((value >> 16) & 0x0000ffff0000ffff);
	value = (value & 0x00000000ffffffff) + ((value >> 32) & 0x00000000ffffffff);
	return value;
#endif
}

/* 最上位ビットの位置を返す (0-based, value == 0 のときは -1 を返す) */
/* floor(log2(value)) と一致する */
inline int msbPosition(uint64_t value) {
	value |= value >> 1;
	value |= value >> 2;
	value |= value >> 4;
	value |= value >> 8;
	value |= value >> 16;
	value |= value >> 32;
	return popcount(value) - 1;
}

/* 最下位ビットの位置を返す (0-based, value == 0 のときは 64 を返す) */
inline int lsbPosition(uint64_t value) {
	value |= value << 1;
	value |= value << 2;
	value |= value << 4;
	value |= value << 8;
	value |= value << 16;
	value |= value << 32;
	return 64 - popcount(value);
}

/* 最上位ビットを取得 (value == 0 のときは 0 を返す) */
inline uint64_t msb(uint64_t value) {
	if (value == 0) return uint64_t();
	return uint64_t(1) << msbPosition(value);
}

/* 最下位ビットを取得 (value == 0 のときは 0 を返す) */
/* TODO: これに関しては明らかに value & -value のほうが速いと思うけどどうしようか */
inline uint64_t lsb(uint64_t value) {
	if (value == 0) return uint64_t();
	return uint64_t(1) << lsbPosition(value);
}

/* ceil(log2(value)) を返す (value == 0 のときは -1 を返す) */
/* 戻り値の範囲は -1 <= ret <= 64 */
inline int ceilLog2(uint64_t value) {
	if (value == 0) return -1;
	return msbPosition(value - 1) + 1;
}

/* ビットの並びを逆転させる */
inline uint64_t bitReverse(uint64_t val) {
	val = (val & 0x5555555555555555) <<  1 | (val & 0xAAAAAAAAAAAAAAAA) >>  1;
	val = (val & 0x3333333333333333) <<  2 | (val & 0xCCCCCCCCCCCCCCCC) >>  2;
	val = (val & 0x0F0F0F0F0F0F0F0F) <<  4 | (val & 0xF0F0F0F0F0F0F0F0) >>  4;
	val = (val & 0x00FF00FF00FF00FF) <<  8 | (val & 0xFF00FF00FF00FF00) >>  8;
	val = (val & 0x0000FFFF0000FFFF) << 16 | (val & 0xFFFF0000FFFF0000) >> 16;
	val = (val & 0x00000000FFFFFFFF) << 32 | (val & 0xFFFFFFFF00000000) >> 32;
	return val;
}

/* 下位 digits ビットについてビットの並びを逆転させる (上位ビットは無視) */
inline uint64_t bitReverse(uint64_t val, int digits) {
	assert(digits > 0 || digits <= 64);
	return bitReverse(val) >> (64 - digits);
}

/* 1 << quantity を返す */
inline uint64_t shift(int quantity) {
	return uint64_t(1) << quantity;
}

/* 整数 sets を集合としてみたとき、element 番目の要素を持つかどうか bool(s & (1 << element)) */
inline bool has(uint64_t sets, int element) {
	return sets & shift(element);
}

template <class T>
void bitReverseArray(std::vector<T>& arr);

/* 配列の要素を添字のビット逆順に並べ替える (配列長は2冪) */
template <class T>
void bitReverseArray(std::vector<T>& arr) {
	const int rank = msbPosition(arr.size());
	const int length = 1 << rank;
	assert(arr.size() == length);

	for (int i = 0; i < length; ++i) {
		int j = bitReverse(i, rank);
		if (j > i) {
			std::swap(arr[i], arr[j]);
		}
	}
}


// 畳み込み演算系

template <class T>
void transformSubsetZeta(std::vector<T>& arr);
template <class T>
void transformInverseSubsetZeta(std::vector<T>& arr);
template <class T>
void transformSupersetZeta(std::vector<T>& arr);
template <class T>
void transformInverseSupersetZeta(std::vector<T>& arr);
template <class T>
void transformHadamard(std::vector<T>& vec);
template <class T>
std::vector<T> convolveOnXor(const std::vector<T>& a, const std::vector<T>& b);

/* 逆変換とまとめられないか？ */
/* 高速ゼータ変換 O(n log n), n は2冪 */
/* b[s] = sum[t is a subset of s] a[t] */
template <class T>
void transformSubsetZeta(std::vector<T>& arr) {
	const uint64_t length = arr.size();
	const int level = msbPosition(length);
	assert(length == shift(level));

	for (int k = 0; k < level; ++k) {
		for (uint64_t s = 0; s < length; ++s) {
			if (has(s, k)) continue;
			arr[s] += arr[s | shift(k)];
		}
	}
}

/* 高速ゼータ変換の逆変換 O(n log n), n は2冪 */
template <class T>
void transformInverseSubsetZeta(std::vector<T>& arr) {
	const uint64_t length = arr.size();
	const int level = msbPosition(length);
	assert(length == (uint64_t(1) << level));

	for (int k = 0; k < level; ++k) {
		for (uint64_t s = 0; s < length; ++s) {
			if (has(s, k)) continue;
			arr[s | shift(k)] -= arr[s];
		}
	}
}

/* 高速(上位集合)ゼータ変換 O(n log n), n は2冪 */
/* b[s] = sum[t is a superset of s] a[t] */
template <class T>
void transformSupersetZeta(std::vector<T>& arr) {
	const uint64_t length = arr.size();
	const int level = msbPosition(length);
	assert(length == (uint64_t(1) << level));

	for (int k = 0; k < level; ++k) {
		for (uint64_t s = 0; s < length; ++s) {
			if (has(s, k)) continue;
			arr[s | shift(k)] += arr[s];
		}
	}
}

/* 高速(上位集合)ゼータ変換の逆変換 O(n log n), n は2冪 */
template <class T>
void transformInverseSupersetZeta(std::vector<T>& arr) {
	const uint64_t length = arr.size();
	const int level = msbPosition(length);
	assert(length == (uint64_t(1) << level));

	for (int k = 0; k < level; ++k) {
		for (uint64_t s = 0; s < length; ++s) {
			if (has(s, k)) continue;
			arr[s] -= arr[s | shift(k)];
		}
	}
}

/* 高速アダマール変換 O(n log n), n は2冪 (オーバーフロー注意) */
/* オーバーフロー注意 */
template <class T>
void transformHadamard(std::vector<T>& vec) {
	const int rank = msbPosition(vec.size());
	const int length = 1 << rank;
	assert(vec.size() == length);

	/* バタフライ演算 */
	auto butterfly = [](T& a, T& b) -> void {
		T tmp(a);
		a += b;
		b = tmp - b;
	};

	/* FFT 同様に各ビットについてバタフライ演算を回す */
	for (int cnt = 0; cnt < rank; ++cnt) {
		const int blk_max = length >> (cnt + 1);
		const int blk_size = 1 << (cnt + 1);

		for (int blk = 0; blk < blk_max; ++blk) {
			const int idx_max = 1 << cnt;

			for (int idx = 0; idx < idx_max; ++idx) {
				const int left = (blk * blk_size + idx);
				const int right = left + blk_size / 2;

				butterfly(vec[left], vec[right]);
			}
		}
	}
}

/* xor 畳み込み O(n log n), n は2冪 (オーバーフロー注意) */
/* ret[n] = sum [i xor j == n] (a[i] * b[j]) */
template <class T>
std::vector<T> convolveOnXor(const std::vector<T>& a, const std::vector<T>& b) {
	const int rank = msbPosition(a.size());
	const int size = 1 << rank;
	assert(a.size() == b.size());
	assert(a.size() == size);

	std::vector<T> a_copied = a;
	std::vector<T> b_copied = b;
	transformHadamard(a_copied);
	transformHadamard(b_copied);

	std::vector<T> ret(size);
	for (int i = 0; i < size; ++i) {
		ret[i] = a_copied[i] * b_copied[i];
	}
	transformHadamard(ret);

	/* 関数 transformHadamard はユニタリ変換ではないので係数を調整 */
	/* IntMod など -> 逆元を前計算することも可能 */
	for (int i = 0; i < size; ++i) {
		ret[i] /= size;
	}
	return ret;
}



#endif
#ifndef __FFT_H__
#define __FFT_H__

/* updated: 2020-03-28 */

// required: IntMod, Sets

#ifdef ONLY_MY_ENVIR
#include "IntMod.h"
#include "Sets.h"
#endif

#include <cstdint>
#include <iostream>
#include <vector>
#include <utility>
#include <cassert>
#include <algorithm>
#include <complex>

class MyComplex {
	using Real = double;
	Real re;
	Real im;

public:
	MyComplex();
	explicit MyComplex(Real re);
	explicit MyComplex(Real re, Real im);
	explicit MyComplex(const std::complex<double>& other);
	MyComplex(const MyComplex&) = default;
	MyComplex& operator=(const MyComplex&) = default;

	MyComplex& operator+=(const MyComplex& right);
	MyComplex& operator-=(const MyComplex& right);
	MyComplex& operator*=(const MyComplex& right);
	MyComplex& operator/=(const MyComplex& right);

	static MyComplex multiply(const MyComplex& left, const MyComplex& right);
	static MyComplex divide(const MyComplex& left, const MyComplex& right);

	Real real() const;
	Real image() const;
};

MyComplex operator+(const MyComplex& left, const MyComplex& right);
MyComplex operator-(const MyComplex& left, const MyComplex& right);
MyComplex operator*(const MyComplex& left, const MyComplex& right);
MyComplex operator/(const MyComplex& left, const MyComplex& right);

void uMul64(uint64_t a, uint64_t b, uint64_t* ret_lo, uint64_t* ret_hi);

// 逆元は高速化できる？
class IntModSpecial {
	// value_m must be less than Mod
	uint64_t value_m;
public:
	static const uint64_t Mod = 0xFFFFFFFF00000001;

	IntModSpecial();
	explicit IntModSpecial(uint64_t value);
	IntModSpecial(const IntModSpecial&) = default;
	IntModSpecial& operator=(const IntModSpecial&) = default;

	IntModSpecial& operator+=(const IntModSpecial& right);
	IntModSpecial& operator-=(const IntModSpecial& right);
	IntModSpecial& operator*=(const IntModSpecial& right);
	IntModSpecial& operator/=(const IntModSpecial& right);

	uint64_t get() const;
	IntModSpecial power(uint64_t exp) const;
	IntModSpecial inverse() const;
	static IntModSpecial primitiveRoot(int rank);
};

IntModSpecial operator+(const IntModSpecial& left, const IntModSpecial& right);
IntModSpecial operator-(const IntModSpecial& left, const IntModSpecial& right);
IntModSpecial operator*(const IntModSpecial& left, const IntModSpecial& right);
IntModSpecial operator/(const IntModSpecial& left, const IntModSpecial& right);

std::istream& operator>>(std::istream& ist, IntModSpecial& object);
std::ostream& operator<<(std::ostream& ost, const IntModSpecial& object);

//uint32_t bitReverse(uint32_t val);
//uint32_t bitReverse(uint32_t val, int digits);

/*------------------------------------------*/

/* もう少し良い方法はないか */
template <uint32_t Mod>
struct Divisors {
	static constexpr int calcSize() {
		uint32_t n = Mod;
		int count = 0;
		for (uint32_t i = 2; i * i <= n; ++i) {
			if (n % i == 0) {
				++count;
			}

			while (n % i == 0) {
				n /= i;
			}
		}

		if (n > 1) {
			++count;
		}

		return count;
	}

	constexpr static int size = calcSize();
	uint32_t list[size];

	constexpr Divisors() : list() {
		uint32_t n = Mod;
		int index = 0;
		for (uint32_t i = 2; i * i <= n; ++i) {
			if (n % i == 0) {
				list[index] = i;
				++index;
			}

			while (n % i == 0) {
				n /= i;
			}
		}

		if (n > 1) {
			list[index] = n;
			++index;
		}
	}
};

constexpr uint32_t powMod(uint32_t base, uint64_t exp, uint32_t modulus) {
	uint64_t product = 1;
	uint64_t multiplier = base;
	while (exp > 0) {
		if (exp % 2 == 1) {
			product = product * multiplier % modulus;
		}
		multiplier = multiplier * multiplier % modulus;
		exp /= 2;
	}
	return uint32_t(product);
}

/* 素数 Mod に対して原始根を一つ求める */
template <uint32_t Mod>
constexpr uint32_t primitiveRoot() {
	const uint32_t order = Mod - 1;
	Divisors<order> divisors;
	for (int i = 2; i < Mod; ++i) {
		bool is_primitive = true;
		for (uint32_t factor : divisors.list) {
			if (powMod(i, order / factor, Mod) == 1) {
				is_primitive = false;
				break;
			}
		}

		if (is_primitive) {
			return i;
		}
	}
}

template <class T, uint32_t Rank>
struct FFTRing {
	static_assert(Rank == 0 && Rank != 0, "Type T not supported");
	// T : 値型(複素数や有限体など)
	// root() : 1 の length 乗根
};

template <uint32_t Mod, uint32_t Rank>
struct FFTRing<IntMod<Mod>, Rank> {
	using T = IntMod<Mod>;
	static constexpr uint32_t rank = Rank;
	static constexpr uint32_t length = 1u << Rank;
	static constexpr uint32_t order = (Mod - 1) & (uint32_t)(-int32_t(Mod - 1));
	static constexpr uint32_t primitive = primitiveRoot<Mod>();
	static_assert(Rank < 32 && length <= order, "array too large");

	static T root() {
		const uint32_t coeff = (Mod - 1) / order;
		return T(primitive).power(uint64_t(coeff) * (order / length));
	}
};

template <uint32_t Rank>
struct FFTRing<IntModSpecial, Rank> {
	using T = IntModSpecial;
	static constexpr uint32_t rank = Rank;
	static constexpr uint32_t length = 1u << Rank;
	static_assert(Rank < 32, "array too long");

	static T root() {
		return T::primitiveRoot(Rank);
	}
};

template <uint32_t Rank>
struct FFTRing<MyComplex, Rank> {
	using T = MyComplex;
	static constexpr uint32_t rank = Rank;
	static constexpr uint32_t length = 1u << Rank;

	static T root() {
		const double PI = std::acos(-1);
		return T(std::polar(1.0, 2 * PI / length));
	}
};

/*----------------------------------------*/
/*--------   FourierTransformer   --------*/
/*----------------------------------------*/

template <class T, uint32_t Rank>
class FourierTransformer {
	using Ring = FFTRing<T, Rank>;

	static constexpr uint32_t length = 1u << Rank;
	static_assert(Rank < 32, "array too long");

	const T length_inv;

	// 2つ分配列をとるのは無駄かもしれない
	const std::vector<T> roots;
	const std::vector<T> roots_inv;

public:
	FourierTransformer();

private:
	std::vector<T> calcTable(bool inverse);

public:
	std::vector<T> generateArray() const;
	void execute(std::vector<T>& array) const;
	void inverse(std::vector<T>& array) const;
	std::vector<T> convolve(const std::vector<T>& a, const std::vector<T>& b) const;

private:
	//void bitReverseArray(std::vector<T>& arr) const;
	void butterfly(T& a, T& b, const T& root) const;
	void transform(std::vector<T>& arr, const std::vector<T>& roots) const;
};

template <class T, uint32_t Rank>
FourierTransformer<T, Rank>::FourierTransformer()
	: length_inv(T(1) / T(length)), 
	roots(calcTable(false)),
	roots_inv(calcTable(true)) {
}

template <class T, uint32_t Rank>
std::vector<T> FourierTransformer<T, Rank>::calcTable(bool inverse) {
	std::vector<T> table(length + 1);
	const T root = Ring::root();

	table[0] = T(1);
	for (int i = 1; i <= length; ++i) {
		table[i] = table[i - 1] * root;
	}

	if (inverse) {
		std::reverse(table.begin(), table.end());
	}
	return table;
}

template <class T, uint32_t Rank>
std::vector<T> FourierTransformer<T, Rank>::generateArray() const {
	return std::vector<T>(length);
}

template <class T, uint32_t Rank>
void FourierTransformer<T, Rank>::execute(std::vector<T>& array) const {
	assert(array.size() == length);
	transform(array, roots);
}

template <class T, uint32_t Rank>
void FourierTransformer<T, Rank>::inverse(std::vector<T>& array) const {
	assert(array.size() == length);
	transform(array, roots_inv);

	for (T& element : array) {
		element *= length_inv;
	}
}

template <class T, uint32_t Rank>
std::vector<T> FourierTransformer<T, Rank>::convolve(const std::vector<T>& a, const std::vector<T>& b) const {
	assert(a.size() == length && b.size() == length);

	std::vector<T> a_copied = a;
	std::vector<T> b_copied = b;

	execute(a_copied);
	execute(b_copied);

	for (int i = 0; i < length; ++i) {
		a_copied[i] *= b_copied[i];
	}

	inverse(a_copied);
	return a_copied;
}

//template <class T, uint32_t Rank>
//void FourierTransformer<T, Rank>::bitReverseArray(std::vector<T>& arr) const {
//	for (uint32_t i = 0; i < length; ++i) {
//		uint32_t j = bitReverse(i, Rank);
//		if (j > i) {
//			std::swap(arr[i], arr[j]);
//		}
//	}
//}

template <class T, uint32_t Rank>
void FourierTransformer<T, Rank>::butterfly(T& a, T& b, const T& root) const {
	T temp = b * root;
	b = a - temp;
	a += temp;
}

/* FFT の本体 */
template <class T, uint32_t Rank>
void FourierTransformer<T, Rank>::transform(std::vector<T>& arr, const std::vector<T>& roots) const {
	bitReverseArray(arr);
	for (int cnt = 0; cnt < Rank; ++cnt) {
		const int blk_max = length >> (cnt + 1);
		const int blk_size = 1 << (cnt + 1);
		const int root_interval = blk_max;

		for (int blk = 0; blk < blk_max; ++blk) {
			const int idx_max = 1 << cnt;

			for (int idx = 0; idx < idx_max; ++idx) {
				const int left = (blk * blk_size + idx);
				const int right = left + blk_size / 2;

				butterfly(arr[left], arr[right], roots[idx * root_interval]);
			}
		}
	}
}

#endif

#define _CRT_SECURE_NO_WARNINGS
#define _SCL_SECURE_NO_WARNINGS
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cassert>
#include <iostream>
#include <string>
#include <vector>
#include <list>
#include <utility>
#include <algorithm>
#include <functional>
#include <cmath>
#include <stack>
#include <queue>
#include <set>
#include <map>
#include <iomanip>
#include <sstream>
#include <bitset>
#include <limits>
#include <numeric>
#include <valarray>
#include <fstream>
#include <array>
#include <random>
#include <unordered_set>
#include <unordered_map>

using namespace std;
using uint = uint32_t;
using LL = int64_t;
using ULL = uint64_t;
using PP = pair<LL, LL>;
template <typename T> using PriorityQ = priority_queue<T, vector<T>, greater<T> >;
#define REP(i, a, n) for(LL i = (a), i##_max_ = (n); i < i##_max_; ++i)
#define REM(i, a, n) for(LL i = (LL)(n) - 1, i##_min_ = (a); i >= i##_min_; --i)
#define FLOAT fixed << setprecision(16)
#define SPEEDUP { cin.tie(NULL); ios::sync_with_stdio(false); }
const int INF = 0x3FFFFFFF;
const LL INFLL = 0x3FFFFFFF3FFFFFFF;
const double INFD = 1.0e+308;
const double EPS = 1.0e-9;

void YesNo(bool b) { cout << (b ? "Yes" : "No") << endl; }
void YESNO(bool b) { cout << (b ? "YES" : "NO") << endl; }
template <class T, class U> istream& operator>>(istream& ist, pair<T, U>& right) { return ist >> right.first >> right.second; }
template <class T, class U> ostream& operator<<(ostream& ost, const pair<T, U>& right) { return ost << right.first << ' ' << right.second; }
template <class T, class TCompatible, size_t N> void Fill(T(&dest)[N], const TCompatible& val) { fill(dest, dest + N, val); }
template <class T, class TCompatible, size_t M, size_t N> void Fill(T(&dest)[M][N], const TCompatible& val) { for (int i = 0; i < M; ++i) Fill(dest[i], val); }
template <class T> T Next() { T buf; cin >> buf; return buf; }
istream& Ignore(istream& ist) { string s; ist >> s; return ist; }
bool Inside(int i, int j, int h, int w) { return i >= 0 && i < h&& j >= 0 && j < w; }

#ifdef ONLY_MY_ENVIR
#include "Accumulator.h"
#include "Algebraic.h"
#include "BinaryMatrix.h"
#include "BinaryTree.h"
#include "Bipartite.h"
#include "BIT.h"
#include "Compressor.h"
#include "Decompositions.h"
#include "DiscreteLog.h"
#include "DynamicMod.h"
#include "Exponential.h"
#include "Factorization.h"
#include "FFT.h"
#include "FlowSolver.h"
#include "Geometric2D.h"
#include "Geometric2DFloat.h"
#include "Geometric3D.h"
#include "Geometric3DFloat.h"
#include "Graph.h"
#include "GraphUtil.h"
#include "Interpolation.h"
#include "IntMod.h"
#include "KDTree.h"
#include "LazySegmentTree.h"
#include "LIS.h"
#include "List.h"
#include "Math.h"
#include "MathUtil.h"
#include "Matrix.h"
#include "MinCostFlowSolver.h"
#include "MinMax.h"
#include "Numbers.h"
#include "Optimize.h"
#include "Permutation.h"
#include "Polynomial.h"
#include "Position.h"
#include "Random.h"
#include "Range.h"
#include "Rational.h"
#include "RuntimeMod.h"
#include "SegmentTree.h"
#include "SegmentTree2D.h"
#include "Sets.h"
#include "Shortest.h"
#include "SlidingWindow.h"
#include "SpanningTree.h"
#include "StringSearching.h"
#include "SuffixArray.h"
#include "SwitchList.h"
#include "Timer.h"
#include "Tree.h"
#include "TreeUtil.h"
#include "UnionFind.h"
#include "Util.h"
#include "VectorUtil.h"
#endif

#ifdef __GNUC__
typedef __int128 LLL;
istream& operator>>(istream& ist, __int128& val) { LL tmp; ist >> tmp; val = tmp; return ist; }
ostream& operator<<(ostream& ost, __int128 val) { LL tmp = val; ost << tmp; return ost; }
#endif

const int P = 200003;
using MP = IntMod<P>;

vector<LL> explist() {
	vector<LL> ret(P);
	MP prod = 1;
	MP g = 2;
	REP(i, 0, P - 1) {
		ret[i] = prod.get();
		prod *= g;
	}
	return ret;
}

vector<LL> loglist() {
	vector<LL> ret(P);
	MP prod = 1;
	MP g = 2;
	REP(i, 0, P - 1) {
		ret[prod.get()] = i;
		prod *= g;
	}
	return ret;
}

int main() {
	int n = Next<int>();
	vector<LL> v(n);
	REP(i, 0, n) v[i] = Next<LL>();
	vector<LL> lg = loglist();
	vector<LL> ex = explist();

	FourierTransformer<IntModSpecial, 19> ft;
	auto a = ft.generateArray();

	REP(i, 0, n) {
		if (v[i] == 0) continue;
		int pos = lg[v[i]];
		a[pos] += IntModSpecial(1);
	}
	auto res = ft.convolve(a, a);

	LL sum = 0;
	REP(i, 0, res.size()) {
		sum += res[i].get() * ex[i % (P - 1)];
	}

	REP(i, 0, n) {
		sum -= v[i] * v[i] % P;
	}

	cout << sum / 2 << endl;
	return 0;
}








Compiler output:
/usr/bin/ld: /tmp/ccDuPxJZ.o: in function `main':
a.cc:(.text+0x4ad): undefined reference to `IntModSpecial::IntModSpecial(unsigned long)'
/usr/bin/ld: a.cc:(.text+0x4d7): undefined reference to `IntModSpecial::operator+=(IntModSpecial const&)'
/usr/bin/ld: a.cc:(.text+0x557): undefined reference to `IntModSpecial::get() const'
/usr/bin/ld: /tmp/ccDuPxJZ.o: in function `FourierTransformer<IntModSpecial, 19u>::FourierTransformer()':
a.cc:(.text._ZN18FourierTransformerI13IntModSpecialLj19EEC2Ev[_ZN18FourierTransformerI13IntModSpecialLj19EEC5Ev]+0x1a): undefined reference to `IntModSpecial::IntModSpecial(unsigned long)'
/usr/bin/ld: a.cc:(.text._ZN18FourierTransformerI13IntModSpecialLj19EEC2Ev[_ZN18FourierTransformerI13IntModSpecialLj19EEC5Ev]+0x2b): undefined reference to `IntModSpecial::IntModSpecial(unsigned long)'
/usr/bin/ld: a.cc:(.text._ZN18FourierTransformerI13IntModSpecialLj19EEC2Ev[_ZN18FourierTransformerI13IntModSpecialLj19EEC5Ev]+0x42): undefined reference to `operator/(IntModSpecial const&, IntModSpecial const&)'
/usr/bin/ld: /tmp/ccDuPxJZ.o: in function `FourierTransformer<IntModSpecial, 19u>::convolve(std::vector<IntModSpecial, std::allocator<IntModSpecial> > const&, std::vector<IntModSpecial, std::allocator<IntModSpecial> > const&) const':
a.cc:(.text._ZNK18FourierTransformerI13IntModSpecialLj19EE8convolveERKSt6vectorIS0_SaIS0_EES6_[_ZNK18FourierTransformerI13IntModSpecialLj19EE8convolveERKSt6vectorIS0_SaIS0_EES6_]+0xf2): undefined reference to `IntModSpecial::operator*=(IntModSpecial const&)'
/usr/bin/ld: /tmp/ccDuPxJZ.o: in function `FourierTransformer<IntModSpecial, 19u>::calcTable(bool)':
a.cc:(.text._ZN18FourierTransformerI13IntModSpecialLj19EE9calcTableEb[_ZN18FourierTransformerI13IntModSpecialLj19EE9calcTableEb]+0x58): undefined reference to `IntModSpecial::IntModSpecial(unsigned long)'
/usr/bin/ld: a.cc:(.text._ZN18FourierTransformerI13IntModSpecialLj19EE9calcTableEb[_ZN18FourierTransformerI13IntModSpecialLj19EE9calcTableEb]+0xa3): undefined reference to `operator*(IntModSpecial const&, IntModSpecial const&)'
/usr/bin/ld: /tmp/ccDuPxJZ.o: in function `FourierTransformer<IntModSpecial, 19u>::inverse(std::vector<IntModSpecial, std::allocator<IntModSpecial> >&) const':
a.cc:(.text._ZNK18FourierTransformerI13IntModSpecialLj19EE7inverseERSt6vectorIS0_SaIS0_EE[_ZNK18FourierTransformerI13IntModSpecialLj19EE7inverseERSt6vectorIS0_SaIS0_EE]+0xb0): undefined reference to `IntModSpecial::operator*=(IntModSpecial const&)'
/usr/bin/ld: /tmp/ccDuPxJZ.o: in function `FFTRing<IntModSpecial, 19u>::root()':
a.cc:(.text._ZN7FFTRingI13IntModSpecialLj19EE4rootEv[_ZN7FFTRingI13IntModSpecialLj19EE4rootEv]+0xa): undefined reference to `IntModSpecial::primitiveRoot(int)'
/usr/bin/ld: /tmp/ccDuPxJZ.o: in function `FourierTransformer<IntModSpecial, 19u>::butterfly(IntModSpecial&, IntModSpecial&, IntModSpecial const&) const':
a.cc:(.text._ZNK18FourierTransformerI13IntModSpecialLj19EE9butterflyERS0_S2_RKS0_[_ZNK18FourierTransformerI13IntModSpecialLj19EE9butterflyERS0_S2_RKS0_]+0x27): undefined reference to `operator*(IntModSpecial const&, IntModSpecial const&)'
/usr/bin/ld: a.cc:(.text._ZNK18FourierTransformerI13IntModSpecialLj19EE9butterflyERS0_S2_RKS0_[_ZNK18FourierTransformerI13IntModSpecialLj19EE9butterflyERS0_S2_RKS0_]+0x3e): undefined reference to `operator-(IntModSpecial const&, IntModSpecial const&)'
/usr/bin/ld: a.cc:(.text._ZNK18FourierTransformerI13IntModSpecialLj19EE9butterflyERS0_S2_RKS0_[_ZNK18FourierTransformerI13IntModSpecialLj19EE9butterflyERS0_S2_RKS0_]+0x5b): undefined reference to `IntModSpecial::operator+=(IntModSpecial const&)'
/usr/bin/ld: /tmp/ccDuPxJZ.o: in function `void std::_Construct<IntModSpecial>(IntModSpecial*)':
a.cc:(.text._ZSt10_ConstructI13IntModSpecialJEEvPT_DpOT0_[_ZSt10_ConstructI13IntModSpecialJEEvPT_DpOT0_]+0x31): undefined reference to `IntModSpecial::IntModSpecial()'
collect2: error: ld returned 1 exit status