Cyclotomic Polynomial
(math/cyclotomic_polynomial.hpp)
- View this file on GitHub
- Last update: 2026-07-03 15:39:11+09:00
- Include:
#include "math/cyclotomic_polynomial.hpp"
Overview
The $n$-th cyclotomic polynomial $\Phi_n(x)$ is the monic polynomial whose roots are the primitive $n$-th roots of unity. It is characterized by
\[x^n - 1 = \prod_{d \mid n} \Phi_d(x).\]This header constructs its coefficients using the Möbius product
\[\Phi_n(x) = \prod_{d \mid n} (1-x^d)^{\mu(n/d)}\]as a truncated formal power series. Only sparse multiplication and exact division by factors of the form $1-x^d$ are needed.
Function
template <class T = long long>
std::vector<T> cyclotomic_polynomial(std::uint64_t index);
index must be positive. The result stores coefficients in ascending degree
order: element i is the coefficient of $x^i$. Its size is
$\varphi(\mathtt{index})+1$.
T must be constructible from -1, 0, and 1, and support addition and
subtraction. The coefficients are integers, so no multiplicative inverse or
prime modulus is required. With the default long long, every coefficient
must fit in that type.
| Function | Description | Complexity |
|---|---|---|
cyclotomic_polynomial<T>(index) |
Returns the coefficients of $\Phi_{\mathtt{index}}(x)$. | $O(2^{\omega(n)}\varphi(n))$ time and $O(\varphi(n))$ memory, plus 64-bit factorization |
Here $\omega(n)$ is the number of distinct prime factors. Factorization uses the library’s deterministic Miller–Rabin and Pollard–Rho implementation.
For repunits, the defining identity immediately gives
\[R_n(x) = \frac{x^n-1}{x-1} = \prod_{\substack{d \mid n \\ d>1}} \Phi_d(x).\]Example
#include "math/cyclotomic_polynomial.hpp"
#include <iostream>
#include <vector>
int main() {
std::vector<long long> polynomial =
m1une::math::cyclotomic_polynomial(12);
for (long long coefficient : polynomial) {
std::cout << coefficient << ' ';
}
std::cout << '\n'; // 1 0 -1 0 1
}
Depends on
Required by
Verified with
Code
#ifndef M1UNE_MATH_CYCLOTOMIC_POLYNOMIAL_HPP
#define M1UNE_MATH_CYCLOTOMIC_POLYNOMIAL_HPP 1
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <vector>
#include "prime_factorization.hpp"
namespace m1une {
namespace math {
template <class T = long long>
std::vector<T> cyclotomic_polynomial(std::uint64_t index) {
assert(index >= 1);
if (index == 1) return {T(-1), T(1)};
const std::vector<std::pair<std::uint64_t, int>> factors =
prime_factorize(index);
std::uint64_t degree = index;
for (const auto& factor : factors) {
degree = degree / factor.first * (factor.first - 1);
}
assert(degree < std::numeric_limits<std::size_t>::max());
std::vector<T> result(static_cast<std::size_t>(degree) + 1, T(0));
result[0] = T(1);
const std::size_t subset_count = std::size_t(1) << factors.size();
for (std::size_t mask = 0; mask < subset_count; mask++) {
std::uint64_t exponent = index;
bool negative_mobius = false;
for (std::size_t i = 0; i < factors.size(); i++) {
if ((mask >> i) & 1) {
exponent /= factors[i].first;
negative_mobius = !negative_mobius;
}
}
if (exponent > degree) continue;
const std::size_t shift = static_cast<std::size_t>(exponent);
if (negative_mobius) {
// Divide by 1 - x^shift as a truncated formal power series.
for (std::size_t i = shift; i <= degree; i++) {
result[i] += result[i - shift];
}
} else {
// Multiply by 1 - x^shift.
for (std::size_t i = static_cast<std::size_t>(degree);
i >= shift;
i--) {
result[i] -= result[i - shift];
if (i == shift) break;
}
}
}
return result;
}
} // namespace math
} // namespace m1une
#endif // M1UNE_MATH_CYCLOTOMIC_POLYNOMIAL_HPP#line 1 "math/cyclotomic_polynomial.hpp"
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <vector>
#line 1 "math/prime_factorization.hpp"
#include <algorithm>
#line 7 "math/prime_factorization.hpp"
#include <numeric>
#include <utility>
#line 10 "math/prime_factorization.hpp"
namespace m1une {
namespace math {
namespace internal {
inline uint64_t multiply_mod(uint64_t a, uint64_t b, uint64_t mod) {
return static_cast<uint64_t>(static_cast<unsigned __int128>(a) * b % mod);
}
inline uint64_t power_mod(uint64_t base, uint64_t exponent, uint64_t mod) {
uint64_t result = 1;
while (exponent > 0) {
if (exponent & 1) result = multiply_mod(result, base, mod);
base = multiply_mod(base, base, mod);
exponent >>= 1;
}
return result;
}
inline uint64_t pollard_random() {
static uint64_t state = 0x123456789abcdef0ULL;
state += 0x9e3779b97f4a7c15ULL;
uint64_t value = state;
value = (value ^ (value >> 30)) * 0xbf58476d1ce4e5b9ULL;
value = (value ^ (value >> 27)) * 0x94d049bb133111ebULL;
return value ^ (value >> 31);
}
} // namespace internal
inline bool is_prime(uint64_t value) {
if (value < 2) return false;
for (uint64_t prime : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
if (value % prime == 0) return value == prime;
}
uint64_t odd_part = value - 1;
int power_of_two = 0;
while ((odd_part & 1) == 0) {
odd_part >>= 1;
power_of_two++;
}
for (uint64_t base : {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}) {
if (base % value == 0) continue;
uint64_t x = internal::power_mod(base % value, odd_part, value);
if (x == 1 || x == value - 1) continue;
bool composite = true;
for (int i = 1; i < power_of_two; i++) {
x = internal::multiply_mod(x, x, value);
if (x == value - 1) {
composite = false;
break;
}
}
if (composite) return false;
}
return true;
}
namespace internal {
inline uint64_t pollard_rho(uint64_t value) {
for (uint64_t prime : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
if (value % prime == 0) return prime;
}
while (true) {
const uint64_t constant = pollard_random() % (value - 1) + 1;
uint64_t y = pollard_random() % (value - 1) + 1;
uint64_t x = 0;
uint64_t saved_y = 0;
uint64_t gcd = 1;
uint64_t segment_length = 1;
auto advance = [&](uint64_t current) {
return static_cast<uint64_t>(
(static_cast<unsigned __int128>(multiply_mod(current, current, value)) + constant) % value);
};
while (gcd == 1) {
x = y;
for (uint64_t i = 0; i < segment_length; i++) y = advance(y);
for (uint64_t offset = 0; offset < segment_length && gcd == 1; offset += 128) {
saved_y = y;
uint64_t product = 1;
const uint64_t block = std::min<uint64_t>(128, segment_length - offset);
for (uint64_t i = 0; i < block; i++) {
y = advance(y);
const uint64_t difference = x > y ? x - y : y - x;
product = multiply_mod(product, difference, value);
}
gcd = std::gcd(product, value);
}
segment_length <<= 1;
}
if (gcd == value) {
do {
saved_y = advance(saved_y);
const uint64_t difference = x > saved_y ? x - saved_y : saved_y - x;
gcd = std::gcd(difference, value);
} while (gcd == 1);
}
if (gcd != value) return gcd;
}
}
inline void factor_recursively(uint64_t value, std::vector<uint64_t>& factors) {
if (value == 1) return;
if (is_prime(value)) {
factors.push_back(value);
return;
}
const uint64_t divisor = pollard_rho(value);
factor_recursively(divisor, factors);
factor_recursively(value / divisor, factors);
}
} // namespace internal
inline std::vector<uint64_t> prime_factors(uint64_t value) {
assert(value >= 1);
std::vector<uint64_t> result;
internal::factor_recursively(value, result);
std::sort(result.begin(), result.end());
return result;
}
inline std::vector<std::pair<uint64_t, int>> prime_factorize(uint64_t value) {
std::vector<uint64_t> factors = prime_factors(value);
std::vector<std::pair<uint64_t, int>> result;
for (uint64_t prime : factors) {
if (result.empty() || result.back().first != prime) {
result.emplace_back(prime, 1);
} else {
result.back().second++;
}
}
return result;
}
inline std::vector<uint64_t> divisors(uint64_t value) {
std::vector<uint64_t> result = {1};
for (const auto& factor : prime_factorize(value)) {
const int current_size = int(result.size());
uint64_t power = 1;
for (int exponent = 1; exponent <= factor.second; exponent++) {
power *= factor.first;
for (int i = 0; i < current_size; i++) {
result.push_back(result[i] * power);
}
}
}
std::sort(result.begin(), result.end());
return result;
}
inline uint64_t euler_phi(uint64_t value) {
assert(value >= 1);
uint64_t result = value;
for (const auto& factor : prime_factorize(value)) {
result = result / factor.first * (factor.first - 1);
}
return result;
}
inline int mobius(uint64_t value) {
assert(value >= 1);
int result = 1;
for (const auto& factor : prime_factorize(value)) {
if (factor.second >= 2) return 0;
result = -result;
}
return result;
}
} // namespace math
} // namespace m1une
#line 11 "math/cyclotomic_polynomial.hpp"
namespace m1une {
namespace math {
template <class T = long long>
std::vector<T> cyclotomic_polynomial(std::uint64_t index) {
assert(index >= 1);
if (index == 1) return {T(-1), T(1)};
const std::vector<std::pair<std::uint64_t, int>> factors =
prime_factorize(index);
std::uint64_t degree = index;
for (const auto& factor : factors) {
degree = degree / factor.first * (factor.first - 1);
}
assert(degree < std::numeric_limits<std::size_t>::max());
std::vector<T> result(static_cast<std::size_t>(degree) + 1, T(0));
result[0] = T(1);
const std::size_t subset_count = std::size_t(1) << factors.size();
for (std::size_t mask = 0; mask < subset_count; mask++) {
std::uint64_t exponent = index;
bool negative_mobius = false;
for (std::size_t i = 0; i < factors.size(); i++) {
if ((mask >> i) & 1) {
exponent /= factors[i].first;
negative_mobius = !negative_mobius;
}
}
if (exponent > degree) continue;
const std::size_t shift = static_cast<std::size_t>(exponent);
if (negative_mobius) {
// Divide by 1 - x^shift as a truncated formal power series.
for (std::size_t i = shift; i <= degree; i++) {
result[i] += result[i - shift];
}
} else {
// Multiply by 1 - x^shift.
for (std::size_t i = static_cast<std::size_t>(degree);
i >= shift;
i--) {
result[i] -= result[i - shift];
if (i == shift) break;
}
}
}
return result;
}
} // namespace math
} // namespace m1une