Cyclotomic polynomial: Difference between revisions
(Added C#) |
|||
Line 567: | Line 567: | ||
Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient: |
Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient: |
||
CP[1] has coefficient with magnitude = 1 |
|||
CP[105] has coefficient with magnitude = 2 |
|||
CP[385] has coefficient with magnitude = 3 |
|||
CP[1365] has coefficient with magnitude = 4 |
|||
CP[1785] has coefficient with magnitude = 5 |
|||
CP[2805] has coefficient with magnitude = 6 |
|||
CP[3135] has coefficient with magnitude = 7 |
|||
CP[6545] has coefficient with magnitude = 8 |
|||
CP[6545] has coefficient with magnitude = 9 |
|||
CP[10465] has coefficient with magnitude = 10</pre> |
|||
=={{header|C sharp}}== |
|||
{{trans|Java}} |
|||
{{works with|C sharp|8}} |
|||
<lang csharp>using System; |
|||
using System.Collections; |
|||
using System.Collections.Generic; |
|||
using System.Linq; |
|||
using IntMap = System.Collections.Generic.Dictionary<int, int>; |
|||
public static class CyclotomicPolynomial |
|||
{ |
|||
public static void Main2() { |
|||
Console.WriteLine("Task 1: Cyclotomic polynomials for n <= 30:"); |
|||
for (int i = 1; i <= 30; i++) { |
|||
var p = GetCyclotomicPolynomial(i); |
|||
Console.WriteLine($"CP[{i}] = {p.ToString()}"); |
|||
} |
|||
Console.WriteLine(); |
|||
Console.WriteLine("Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:"); |
|||
for (int i = 1, n = 0; i <= 10; i++) { |
|||
while (true) { |
|||
n++; |
|||
var p = GetCyclotomicPolynomial(n); |
|||
if (p.Any(t => Math.Abs(t.Coefficient) == i)) { |
|||
Console.WriteLine($"CP[{n}] has coefficient with magnitude = {i}"); |
|||
n--; |
|||
break; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
private const int MaxFactors = 100_000; |
|||
private const int Algorithm = 2; |
|||
private static readonly Term x = new Term(1, 1); |
|||
private static readonly Dictionary<int, Polynomial> polyCache = |
|||
new Dictionary<int, Polynomial> { [1] = x - 1 }; |
|||
private static readonly Dictionary<int, IntMap> factorCache = |
|||
new Dictionary<int, IntMap> { [2] = new IntMap { [2] = 1 } }; |
|||
private static Polynomial GetCyclotomicPolynomial(in int n) { |
|||
if (polyCache.TryGetValue(n, out var result)) return result; |
|||
var factors = GetFactors(n); |
|||
if (factors.ContainsKey(n)) { //n is prime |
|||
result = new Polynomial(from exp in ..n select x[exp]); |
|||
} else if (factors.Count == 2 && factors.Contains(2, 1) && factors.Contains(n/2, 1)) { //n = 2p |
|||
result = new Polynomial(from i in ..(n/2) select (IsOdd(i) ? -x : x)[i]); |
|||
} else if (factors.Count == 1 && factors.TryGetValue(2, out int h)) { //n = 2^h |
|||
result = x[1<<(h-1)] + 1; |
|||
} else if (factors.Count == 1 && !factors.ContainsKey(n)) { // n = p^k |
|||
(int p, int k) = factors.First(); |
|||
result = new Polynomial(from i in ..p select x[i * (int)Math.Pow(p, k-1)]); |
|||
} else if (factors.Count == 2 && factors.ContainsKey(2)) { //n = 2^h * p^k |
|||
(int p, int k) = factors.First(entry => entry.Key != 2); |
|||
int twoExp = 1 << (factors[2] - 1); |
|||
result = new Polynomial(from i in ..p select (IsOdd(i) ? -x : x)[i * twoExp * (int)Math.Pow(p, k-1)]); |
|||
} else if (factors.ContainsKey(2) && IsOdd(n/2) && n/2 > 1) { // CP(2m)[x] = CP[-m][x], n is odd > 1 |
|||
Polynomial cycloDiv2 = GetCyclotomicPolynomial(n/2); |
|||
result = new Polynomial(from term in cycloDiv2 select IsOdd(term.Exponent) ? -term : term); |
|||
#pragma warning disable CS0162 |
|||
} else if (Algorithm == 0) { |
|||
var divisors = GetDivisors(n); |
|||
result = x[n] - 1; |
|||
foreach (int d in divisors) result /= GetCyclotomicPolynomial(d); |
|||
} else if (Algorithm == 1) { |
|||
var divisors = GetDivisors(n).ToList(); |
|||
int maxDivisor = divisors.Max(); |
|||
result = (x[n] - 1) / (x[maxDivisor] - 1); |
|||
foreach (int d in divisors.Where(div => maxDivisor % div == 0)) { |
|||
result /= GetCyclotomicPolynomial(d); |
|||
} |
|||
} else if (Algorithm == 2) { |
|||
int m = 1; |
|||
result = GetCyclotomicPolynomial(m); |
|||
var primes = factors.Keys.ToList(); |
|||
primes.Sort(); |
|||
foreach (int prime in primes) { |
|||
var cycloM = result; |
|||
result = new Polynomial(from term in cycloM select term.Coefficient * x[term.Exponent * prime]); |
|||
result /= cycloM; |
|||
m *= prime; |
|||
} |
|||
int s = n / m; |
|||
result = new Polynomial(from term in result select term.Coefficient * x[term.Exponent * s]); |
|||
#pragma warning restore CS0162 |
|||
} else { |
|||
throw new InvalidOperationException("Invalid algorithm"); |
|||
} |
|||
polyCache[n] = result; |
|||
return result; |
|||
} |
|||
private static bool IsOdd(int i) => (i & 1) != 0; |
|||
private static bool Contains(this IntMap map, int key, int value) => map.TryGetValue(key, out int v) && v == value; |
|||
private static int GetOrZero(this IntMap map, int key) => map.TryGetValue(key, out int v) ? v : 0; |
|||
private static IEnumerable<T> Select<T>(this Range r, Func<int, T> f) => Enumerable.Range(r.Start.Value, r.End.Value - r.Start.Value).Select(f); |
|||
private static IntMap GetFactors(in int n) { |
|||
if (factorCache.TryGetValue(n, out var factors)) return factors; |
|||
factors = new IntMap(); |
|||
if (!IsOdd(n)) { |
|||
foreach (var entry in GetFactors(n/2)) factors.Add(entry.Key, entry.Value); |
|||
factors[2] = factors.GetOrZero(2) + 1; |
|||
return Cache(n, factors); |
|||
} |
|||
for (int i = 3; i * i <= n; i+=2) { |
|||
if (n % i == 0) { |
|||
foreach (var entry in GetFactors(n/i)) factors.Add(entry.Key, entry.Value); |
|||
factors[i] = factors.GetOrZero(i) + 1; |
|||
return Cache(n, factors); |
|||
} |
|||
} |
|||
factors[n] = 1; |
|||
return Cache(n, factors); |
|||
} |
|||
private static IntMap Cache(int n, IntMap factors) { |
|||
if (n < MaxFactors) factorCache[n] = factors; |
|||
return factors; |
|||
} |
|||
private static IEnumerable<int> GetDivisors(int n) { |
|||
for (int i = 1; i * i <= n; i++) { |
|||
if (n % i == 0) { |
|||
yield return i; |
|||
int div = n / i; |
|||
if (div != i && div != n) yield return div; |
|||
} |
|||
} |
|||
} |
|||
public sealed class Polynomial : IEnumerable<Term> |
|||
{ |
|||
public Polynomial() { } |
|||
public Polynomial(params Term[] terms) : this(terms.AsEnumerable()) { } |
|||
public Polynomial(IEnumerable<Term> terms) { |
|||
Terms.AddRange(terms); |
|||
Simplify(); |
|||
} |
|||
private List<Term>? terms; |
|||
private List<Term> Terms => terms ??= new List<Term>(); |
|||
public int Count => terms?.Count ?? 0; |
|||
public int Degree => Count == 0 ? -1 : Terms[0].Exponent; |
|||
public int LeadingCoefficient => Count == 0 ? 0 : Terms[0].Coefficient; |
|||
public IEnumerator<Term> GetEnumerator() => Terms.GetEnumerator(); |
|||
IEnumerator IEnumerable.GetEnumerator() => GetEnumerator(); |
|||
public override string ToString() => Count == 0 ? "0" : string.Join(" + ", Terms).Replace("+ -", "- "); |
|||
public static Polynomial operator *(Polynomial p, Term t) => new Polynomial(from s in p select s * t); |
|||
public static Polynomial operator +(Polynomial p, Polynomial q) => new Polynomial(p.Terms.Concat(q.Terms)); |
|||
public static Polynomial operator -(Polynomial p, Polynomial q) => new Polynomial(p.Terms.Concat(q.Terms.Select(t => -t))); |
|||
public static Polynomial operator *(Polynomial p, Polynomial q) => new Polynomial(from s in p from t in q select s * t); |
|||
public static Polynomial operator /(Polynomial p, Polynomial q) => p.Divide(q).quotient; |
|||
public (Polynomial quotient, Polynomial remainder) Divide(Polynomial divisor) { |
|||
if (Degree < 0) return (new Polynomial(), this); |
|||
Polynomial quotient = new Polynomial(); |
|||
Polynomial remainder = this; |
|||
int lcv = divisor.LeadingCoefficient; |
|||
int dv = divisor.Degree; |
|||
while (remainder.Degree >= divisor.Degree) { |
|||
int lcr = remainder.LeadingCoefficient; |
|||
Term div = new Term(lcr / lcv, remainder.Degree - dv); |
|||
quotient.Terms.Add(div); |
|||
remainder += divisor * -div; |
|||
} |
|||
quotient.Simplify(); |
|||
remainder.Simplify(); |
|||
return (quotient, remainder); |
|||
} |
|||
private void Simplify() { |
|||
if (Count < 2) return; |
|||
Terms.Sort((a, b) => -a.CompareTo(b)); |
|||
for (int i = Terms.Count - 1; i > 0; i--) { |
|||
Term s = Terms[i-1]; |
|||
Term t = Terms[i]; |
|||
if (t.Exponent == s.Exponent) { |
|||
Terms[i-1] = new Term(s.Coefficient + t.Coefficient, s.Exponent); |
|||
Terms.RemoveAt(i); |
|||
} |
|||
} |
|||
Terms.RemoveAll(t => t.IsZero); |
|||
} |
|||
} |
|||
public readonly struct Term : IEquatable<Term>, IComparable<Term> |
|||
{ |
|||
public Term(int coefficient, int exponent = 0) => (Coefficient, Exponent) = (coefficient, exponent); |
|||
public Term this[int exponent] => new Term(Coefficient, exponent); //Using x[exp] because x^exp has low precedence |
|||
public int Coefficient { get; } |
|||
public int Exponent { get; } |
|||
public bool IsZero => Coefficient == 0; |
|||
public static Polynomial operator +(Term left, Term right) => new Polynomial(left, right); |
|||
public static Polynomial operator -(Term left, Term right) => new Polynomial(left, -right); |
|||
public static implicit operator Term(int coefficient) => new Term(coefficient); |
|||
public static Term operator -(Term t) => new Term(-t.Coefficient, t.Exponent); |
|||
public static Term operator *(Term left, Term right) => new Term(left.Coefficient * right.Coefficient, left.Exponent + right.Exponent); |
|||
public static bool operator ==(Term left, Term right) => left.Equals(right); |
|||
public static bool operator !=(Term left, Term right) => !left.Equals(right); |
|||
public static bool operator <(Term left, Term right) => left.CompareTo(right) < 0; |
|||
public static bool operator >(Term left, Term right) => left.CompareTo(right) > 0; |
|||
public static bool operator <=(Term left, Term right) => left.CompareTo(right) <= 0; |
|||
public static bool operator >=(Term left, Term right) => left.CompareTo(right) >= 0; |
|||
public bool Equals(Term other) => Exponent == other.Exponent && Coefficient == other.Coefficient; |
|||
public override bool Equals(object? obj) => obj is Term t && Equals(t); |
|||
public override int GetHashCode() => Coefficient.GetHashCode() * 31 + Exponent.GetHashCode(); |
|||
public int CompareTo(Term other) { |
|||
int c = Exponent.CompareTo(other.Exponent); |
|||
if (c != 0) return c; |
|||
return Coefficient.CompareTo(other.Coefficient); |
|||
} |
|||
public override string ToString() => (Coefficient, Exponent) switch { |
|||
(0, _) => "0", |
|||
(_, 0) => $"{Coefficient}", |
|||
(1, 1) => "x", |
|||
(-1, 1) => "-x", |
|||
(_, 1) => $"{Coefficient}x", |
|||
(1, _) => $"x^{Exponent}", |
|||
(-1, _) => $"-x^{Exponent}", |
|||
_ => $"{Coefficient}x^{Exponent}" |
|||
}; |
|||
} |
|||
}</lang> |
|||
{{out}} |
|||
<pre> |
|||
Task 1: Cyclotomic polynomials for n <= 30: |
|||
CP[1] = x - 1 |
|||
CP[2] = x + 1 |
|||
CP[3] = x^2 + x + 1 |
|||
CP[4] = x^2 + 1 |
|||
CP[5] = x^4 + x^3 + x^2 + x + 1 |
|||
CP[6] = x^2 - x + 1 |
|||
CP[7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 |
|||
CP[8] = x^4 + 1 |
|||
CP[9] = x^6 + x^3 + 1 |
|||
CP[10] = x^4 - x^3 + x^2 - x + 1 |
|||
CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 |
|||
CP[12] = x^4 - x^2 + 1 |
|||
CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 |
|||
CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 |
|||
CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 |
|||
CP[16] = x^8 + 1 |
|||
CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 |
|||
CP[18] = x^6 - x^3 + 1 |
|||
CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 |
|||
CP[20] = x^8 - x^6 + x^4 - x^2 + 1 |
|||
CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 |
|||
CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 |
|||
CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 |
|||
CP[24] = x^8 - x^4 + 1 |
|||
CP[25] = x^20 + x^15 + x^10 + x^5 + 1 |
|||
CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 |
|||
CP[27] = x^18 + x^9 + 1 |
|||
CP[28] = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 |
|||
CP[29] = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 |
|||
CP[30] = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 |
|||
Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient: |
|||
CP[1] has coefficient with magnitude = 1 |
CP[1] has coefficient with magnitude = 1 |
||
CP[105] has coefficient with magnitude = 2 |
CP[105] has coefficient with magnitude = 2 |
Revision as of 13:20, 19 June 2020
You are encouraged to solve this task according to the task description, using any language you may know.
The nth Cyclotomic polynomial, for any positive integer n, is the unique irreducible polynomial of largest degree with integer coefficients that is a divisor of x^n − 1, and is not a divisor of x^k − 1 for any k < n.
- Task
- Find and print the first 30 cyclotomic polynomials.
- Find and print the order of the first 10 cyclotomic polynomials that have n or -n as a coefficient.
- See also
- Wikipedia article, Cyclotomic polynomial, showing ways to calculate them.
- The sequence A013594 with the smallest order of cyclotomic polynomial containing n or -n as a coefficient.
C++
<lang cpp>#include <algorithm>
- include <iostream>
- include <initializer_list>
- include <map>
- include <vector>
const int MAX_ALL_FACTORS = 100'000; const int algorithm = 2; int divisions = 0;
//Note: Cyclotomic Polynomials have small coefficients. Not appropriate for general polynomial usage. class Term { private:
long m_coefficient; long m_exponent;
public:
Term(long c, long e) : m_coefficient(c), m_exponent(e) { // empty }
Term(const Term &t) : m_coefficient(t.m_coefficient), m_exponent(t.m_exponent) { // empty }
long coefficient() const { return m_coefficient; }
long degree() const { return m_exponent; }
Term operator -() const { return { -m_coefficient, m_exponent }; }
Term operator *(const Term &rhs) const { return { m_coefficient * rhs.m_coefficient, m_exponent + rhs.m_exponent }; }
Term operator +(const Term &rhs) const { if (m_exponent != rhs.m_exponent) { throw std::runtime_error("Exponents not equal"); } return { m_coefficient + rhs.m_coefficient, m_exponent }; }
friend std::ostream &operator<<(std::ostream &, const Term &);
};
std::ostream &operator<<(std::ostream &os, const Term &t) {
if (t.m_coefficient == 0) { return os << '0'; } if (t.m_exponent == 0) { return os << t.m_coefficient; } if (t.m_coefficient == 1) { if (t.m_exponent == 1) { return os << 'x'; } return os << "x^" << t.m_exponent; } if (t.m_coefficient == -1) { if (t.m_exponent == 1) { return os << "-x"; } return os << "-x^" << t.m_exponent; } if (t.m_exponent == 1) { return os << t.m_coefficient << 'x'; } return os << t.m_coefficient << "x^" << t.m_exponent;
}
class Polynomial { public:
std::vector<Term> polynomialTerms;
Polynomial() { polynomialTerms.push_back({ 0, 0 }); }
Polynomial(std::initializer_list<int> values) { if (values.size() % 2 != 0) { throw std::runtime_error("Length must be even."); }
bool ready = false; long t; for (auto v : values) { if (ready) { polynomialTerms.push_back({ t, v }); } else { t = v; } ready = !ready; } std::sort( polynomialTerms.begin(), polynomialTerms.end(), [](const Term &t, const Term &u) { return u.degree() < t.degree(); } ); }
Polynomial(const std::vector<Term> &termList) { if (termList.size() == 0) { polynomialTerms.push_back({ 0, 0 }); } else { for (auto t : termList) { if (t.coefficient() != 0) { polynomialTerms.push_back(t); } } if (polynomialTerms.size() == 0) { polynomialTerms.push_back({ 0, 0 }); } std::sort( polynomialTerms.begin(), polynomialTerms.end(), [](const Term &t, const Term &u) { return u.degree() < t.degree(); } ); } }
Polynomial(const Polynomial &p) : Polynomial(p.polynomialTerms) { // empty }
long leadingCoefficient() const { return polynomialTerms[0].coefficient(); }
long degree() const { return polynomialTerms[0].degree(); }
bool hasCoefficientAbs(int coeff) { for (auto term : polynomialTerms) { if (abs(term.coefficient()) == coeff) { return true; } } return false; }
Polynomial operator+(const Term &term) const { std::vector<Term> termList; bool added = false; for (size_t index = 0; index < polynomialTerms.size(); index++) { auto currentTerm = polynomialTerms[index]; if (currentTerm.degree() == term.degree()) { added = true; if (currentTerm.coefficient() + term.coefficient() != 0) { termList.push_back(currentTerm + term); } } else { termList.push_back(currentTerm); } } if (!added) { termList.push_back(term); } return Polynomial(termList); }
Polynomial operator*(const Term &term) const { std::vector<Term> termList; for (size_t index = 0; index < polynomialTerms.size(); index++) { auto currentTerm = polynomialTerms[index]; termList.push_back(currentTerm * term); } return Polynomial(termList); }
Polynomial operator+(const Polynomial &p) const { std::vector<Term> termList; int thisCount = polynomialTerms.size(); int polyCount = p.polynomialTerms.size(); while (thisCount > 0 || polyCount > 0) { if (thisCount == 0) { auto polyTerm = p.polynomialTerms[polyCount - 1]; termList.push_back(polyTerm); polyCount--; } else if (polyCount == 0) { auto thisTerm = polynomialTerms[thisCount - 1]; termList.push_back(thisTerm); thisCount--; } else { auto polyTerm = p.polynomialTerms[polyCount - 1]; auto thisTerm = polynomialTerms[thisCount - 1]; if (thisTerm.degree() == polyTerm.degree()) { auto t = thisTerm + polyTerm; if (t.coefficient() != 0) { termList.push_back(t); } thisCount--; polyCount--; } else if (thisTerm.degree() < polyTerm.degree()) { termList.push_back(thisTerm); thisCount--; } else { termList.push_back(polyTerm); polyCount--; } } } return Polynomial(termList); }
Polynomial operator/(const Polynomial &v) { divisions++;
Polynomial q; Polynomial r(*this); long lcv = v.leadingCoefficient(); long dv = v.degree(); while (r.degree() >= v.degree()) { long lcr = r.leadingCoefficient(); long s = lcr / lcv; Term term(s, r.degree() - dv); q = q + term; r = r + v * -term; }
return q; }
friend std::ostream &operator<<(std::ostream &, const Polynomial &);
};
std::ostream &operator<<(std::ostream &os, const Polynomial &p) {
auto it = p.polynomialTerms.cbegin(); auto end = p.polynomialTerms.cend(); if (it != end) { os << *it; it = std::next(it); } while (it != end) { if (it->coefficient() > 0) { os << " + " << *it; } else { os << " - " << -*it; } it = std::next(it); } return os;
}
std::vector<int> getDivisors(int number) {
std::vector<int> divisiors; long root = (long)sqrt(number); for (int i = 1; i <= root; i++) { if (number % i == 0) { divisiors.push_back(i); int div = number / i; if (div != i && div != number) { divisiors.push_back(div); } } } return divisiors;
}
std::map<int, std::map<int, int>> allFactors;
std::map<int, int> getFactors(int number) {
if (allFactors.find(number) != allFactors.end()) { return allFactors[number]; }
std::map<int, int> factors; if (number % 2 == 0) { auto factorsDivTwo = getFactors(number / 2); factors.insert(factorsDivTwo.begin(), factorsDivTwo.end()); if (factors.find(2) != factors.end()) { factors[2]++; } else { factors.insert(std::make_pair(2, 1)); } if (number < MAX_ALL_FACTORS) { allFactors.insert(std::make_pair(number, factors)); } return factors; } long root = (long)sqrt(number); long i = 3; while (i <= root) { if (number % i == 0) { auto factorsDivI = getFactors(number / i); factors.insert(factorsDivI.begin(), factorsDivI.end()); if (factors.find(i) != factors.end()) { factors[i]++; } else { factors.insert(std::make_pair(i, 1)); } if (number < MAX_ALL_FACTORS) { allFactors.insert(std::make_pair(number, factors)); } return factors; } i += 2; } factors.insert(std::make_pair(number, 1)); if (number < MAX_ALL_FACTORS) { allFactors.insert(std::make_pair(number, factors)); } return factors;
}
std::map<int, Polynomial> COMPUTED; Polynomial cyclotomicPolynomial(int n) {
if (COMPUTED.find(n) != COMPUTED.end()) { return COMPUTED[n]; }
if (n == 1) { // Polynomial: x - 1 Polynomial p({ 1, 1, -1, 0 }); COMPUTED.insert(std::make_pair(1, p)); return p; }
auto factors = getFactors(n); if (factors.find(n) != factors.end()) { // n prime std::vector<Term> termList; for (int index = 0; index < n; index++) { termList.push_back({ 1, index }); }
Polynomial cyclo(termList); COMPUTED.insert(std::make_pair(n, cyclo)); return cyclo; } else if (factors.size() == 2 && factors.find(2) != factors.end() && factors[2] == 1 && factors.find(n / 2) != factors.end() && factors[n / 2] == 1) { // n = 2p int prime = n / 2; std::vector<Term> termList; int coeff = -1;
for (int index = 0; index < prime; index++) { coeff *= -1; termList.push_back({ coeff, index }); }
Polynomial cyclo(termList); COMPUTED.insert(std::make_pair(n, cyclo)); return cyclo; } else if (factors.size() == 1 && factors.find(2) != factors.end()) { // n = 2^h int h = factors[2]; std::vector<Term> termList; termList.push_back({ 1, (int)pow(2, h - 1) }); termList.push_back({ 1, 0 });
Polynomial cyclo(termList); COMPUTED.insert(std::make_pair(n, cyclo)); return cyclo; } else if (factors.size() == 1 && factors.find(n) != factors.end()) { // n = p^k int p = 0; int k = 0; for (auto iter = factors.begin(); iter != factors.end(); ++iter) { p = iter->first; k = iter->second; } std::vector<Term> termList; for (int index = 0; index < p; index++) { termList.push_back({ 1, index * (int)pow(p, k - 1) }); }
Polynomial cyclo(termList); COMPUTED.insert(std::make_pair(n, cyclo)); return cyclo; } else if (factors.size() == 2 && factors.find(2) != factors.end()) { // n = 2^h * p^k int p = 0; for (auto iter = factors.begin(); iter != factors.end(); ++iter) { if (iter->first != 2) { p = iter->first; } }
std::vector<Term> termList; int coeff = -1; int twoExp = (int)pow(2, factors[2] - 1); int k = factors[p]; for (int index = 0; index < p; index++) { coeff *= -1; termList.push_back({ coeff, index * twoExp * (int)pow(p, k - 1) }); }
Polynomial cyclo(termList); COMPUTED.insert(std::make_pair(n, cyclo)); return cyclo; } else if (factors.find(2) != factors.end() && ((n / 2) % 2 == 1) && (n / 2) > 1) { // CP(2m)[x] = CP(-m)[x], n odd integer > 1 auto cycloDiv2 = cyclotomicPolynomial(n / 2); std::vector<Term> termList; for (auto term : cycloDiv2.polynomialTerms) { if (term.degree() % 2 == 0) { termList.push_back(term); } else { termList.push_back(-term); } }
Polynomial cyclo(termList); COMPUTED.insert(std::make_pair(n, cyclo)); return cyclo; }
// General Case
if (algorithm == 0) { // slow - uses basic definition auto divisors = getDivisors(n); // Polynomial: (x^n - 1) Polynomial cyclo({ 1, n, -1, 0 }); for (auto i : divisors) { auto p = cyclotomicPolynomial(i); cyclo = cyclo / p; }
COMPUTED.insert(std::make_pair(n, cyclo)); return cyclo; } else if (algorithm == 1) { // Faster. Remove Max divisor (and all divisors of max divisor) - only one divide for all divisors of Max Divisor auto divisors = getDivisors(n); int maxDivisor = INT_MIN; for (auto div : divisors) { maxDivisor = std::max(maxDivisor, div); } std::vector<int> divisorExceptMax; for (auto div : divisors) { if (maxDivisor % div != 0) { divisorExceptMax.push_back(div); } }
// Polynomial: ( x^n - 1 ) / ( x^m - 1 ), where m is the max divisor auto cyclo = Polynomial({ 1, n, -1, 0 }) / Polynomial({ 1, maxDivisor, -1, 0 }); for (int i : divisorExceptMax) { auto p = cyclotomicPolynomial(i); cyclo = cyclo / p; }
COMPUTED.insert(std::make_pair(n, cyclo)); return cyclo; } else if (algorithm == 2) { // Fastest // Let p ; q be primes such that p does not divide n, and q q divides n. // Then CP(np)[x] = CP(n)[x^p] / CP(n)[x] int m = 1; auto cyclo = cyclotomicPolynomial(m); std::vector<int> primes; for (auto iter = factors.begin(); iter != factors.end(); ++iter) { primes.push_back(iter->first); } std::sort(primes.begin(), primes.end()); for (auto prime : primes) { // CP(m)[x] auto cycloM = cyclo; // Compute CP(m)[x^p]. std::vector<Term> termList; for (auto t : cycloM.polynomialTerms) { termList.push_back({ t.coefficient(), t.degree() * prime }); } cyclo = Polynomial(termList) / cycloM; m = m * prime; } // Now, m is the largest square free divisor of n int s = n / m; // Compute CP(n)[x] = CP(m)[x^s] std::vector<Term> termList; for (auto t : cyclo.polynomialTerms) { termList.push_back({ t.coefficient(), t.degree() * s }); }
cyclo = Polynomial(termList); COMPUTED.insert(std::make_pair(n, cyclo)); return cyclo; } else { throw std::runtime_error("Invalid algorithm"); }
}
int main() {
// initialization std::map<int, int> factors; factors.insert(std::make_pair(2, 1)); allFactors.insert(std::make_pair(2, factors));
// rest of main std::cout << "Task 1: cyclotomic polynomials for n <= 30:\n"; for (int i = 1; i <= 30; i++) { auto p = cyclotomicPolynomial(i); std::cout << "CP[" << i << "] = " << p << '\n'; }
std::cout << "Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:\n"; int n = 0; for (int i = 1; i <= 10; i++) { while (true) { n++; auto cyclo = cyclotomicPolynomial(n); if (cyclo.hasCoefficientAbs(i)) { std::cout << "CP[" << n << "] has coefficient with magnitude = " << i << '\n'; n--; break; } } }
return 0;
}</lang>
- Output:
Task 1: cyclotomic polynomials for n <= 30: CP[1] = x - 1 CP[2] = x + 1 CP[3] = x^2 + x + 1 CP[4] = x^2 + 1 CP[5] = x^4 + x^3 + x^2 + x + 1 CP[6] = x^2 - x + 1 CP[7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[8] = x^4 + 1 CP[9] = x^6 + x^3 + 1 CP[10] = x^4 - x^3 + x^2 - x + 1 CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[12] = x^4 - x^2 + 1 CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 CP[16] = x^8 + 1 CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[18] = x^6 - x^3 + 1 CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[20] = x^8 - x^6 + x^4 - x^2 + 1 CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[24] = x^8 - x^4 + 1 CP[25] = x^20 + x^15 + x^10 + x^5 + 1 CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[27] = x^18 + x^9 + 1 CP[28] = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 CP[29] = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[30] = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient: CP[1] has coefficient with magnitude = 1 CP[105] has coefficient with magnitude = 2 CP[385] has coefficient with magnitude = 3 CP[1365] has coefficient with magnitude = 4 CP[1785] has coefficient with magnitude = 5 CP[2805] has coefficient with magnitude = 6 CP[3135] has coefficient with magnitude = 7 CP[6545] has coefficient with magnitude = 8 CP[6545] has coefficient with magnitude = 9 CP[10465] has coefficient with magnitude = 10
C#
<lang csharp>using System; using System.Collections; using System.Collections.Generic; using System.Linq; using IntMap = System.Collections.Generic.Dictionary<int, int>;
public static class CyclotomicPolynomial {
public static void Main2() { Console.WriteLine("Task 1: Cyclotomic polynomials for n <= 30:"); for (int i = 1; i <= 30; i++) { var p = GetCyclotomicPolynomial(i); Console.WriteLine($"CP[{i}] = {p.ToString()}"); } Console.WriteLine();
Console.WriteLine("Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:"); for (int i = 1, n = 0; i <= 10; i++) { while (true) { n++; var p = GetCyclotomicPolynomial(n); if (p.Any(t => Math.Abs(t.Coefficient) == i)) { Console.WriteLine($"CP[{n}] has coefficient with magnitude = {i}"); n--; break; } } } }
private const int MaxFactors = 100_000; private const int Algorithm = 2; private static readonly Term x = new Term(1, 1); private static readonly Dictionary<int, Polynomial> polyCache = new Dictionary<int, Polynomial> { [1] = x - 1 }; private static readonly Dictionary<int, IntMap> factorCache = new Dictionary<int, IntMap> { [2] = new IntMap { [2] = 1 } };
private static Polynomial GetCyclotomicPolynomial(in int n) { if (polyCache.TryGetValue(n, out var result)) return result;
var factors = GetFactors(n); if (factors.ContainsKey(n)) { //n is prime result = new Polynomial(from exp in ..n select x[exp]); } else if (factors.Count == 2 && factors.Contains(2, 1) && factors.Contains(n/2, 1)) { //n = 2p result = new Polynomial(from i in ..(n/2) select (IsOdd(i) ? -x : x)[i]); } else if (factors.Count == 1 && factors.TryGetValue(2, out int h)) { //n = 2^h result = x[1<<(h-1)] + 1; } else if (factors.Count == 1 && !factors.ContainsKey(n)) { // n = p^k (int p, int k) = factors.First(); result = new Polynomial(from i in ..p select x[i * (int)Math.Pow(p, k-1)]); } else if (factors.Count == 2 && factors.ContainsKey(2)) { //n = 2^h * p^k (int p, int k) = factors.First(entry => entry.Key != 2); int twoExp = 1 << (factors[2] - 1); result = new Polynomial(from i in ..p select (IsOdd(i) ? -x : x)[i * twoExp * (int)Math.Pow(p, k-1)]); } else if (factors.ContainsKey(2) && IsOdd(n/2) && n/2 > 1) { // CP(2m)[x] = CP[-m][x], n is odd > 1 Polynomial cycloDiv2 = GetCyclotomicPolynomial(n/2); result = new Polynomial(from term in cycloDiv2 select IsOdd(term.Exponent) ? -term : term); #pragma warning disable CS0162 } else if (Algorithm == 0) { var divisors = GetDivisors(n); result = x[n] - 1; foreach (int d in divisors) result /= GetCyclotomicPolynomial(d); } else if (Algorithm == 1) { var divisors = GetDivisors(n).ToList(); int maxDivisor = divisors.Max(); result = (x[n] - 1) / (x[maxDivisor] - 1); foreach (int d in divisors.Where(div => maxDivisor % div == 0)) { result /= GetCyclotomicPolynomial(d); } } else if (Algorithm == 2) { int m = 1; result = GetCyclotomicPolynomial(m); var primes = factors.Keys.ToList(); primes.Sort(); foreach (int prime in primes) { var cycloM = result; result = new Polynomial(from term in cycloM select term.Coefficient * x[term.Exponent * prime]); result /= cycloM; m *= prime; } int s = n / m; result = new Polynomial(from term in result select term.Coefficient * x[term.Exponent * s]); #pragma warning restore CS0162 } else { throw new InvalidOperationException("Invalid algorithm"); } polyCache[n] = result; return result; }
private static bool IsOdd(int i) => (i & 1) != 0; private static bool Contains(this IntMap map, int key, int value) => map.TryGetValue(key, out int v) && v == value; private static int GetOrZero(this IntMap map, int key) => map.TryGetValue(key, out int v) ? v : 0; private static IEnumerable<T> Select<T>(this Range r, Func<int, T> f) => Enumerable.Range(r.Start.Value, r.End.Value - r.Start.Value).Select(f);
private static IntMap GetFactors(in int n) { if (factorCache.TryGetValue(n, out var factors)) return factors;
factors = new IntMap(); if (!IsOdd(n)) { foreach (var entry in GetFactors(n/2)) factors.Add(entry.Key, entry.Value); factors[2] = factors.GetOrZero(2) + 1; return Cache(n, factors); } for (int i = 3; i * i <= n; i+=2) { if (n % i == 0) { foreach (var entry in GetFactors(n/i)) factors.Add(entry.Key, entry.Value); factors[i] = factors.GetOrZero(i) + 1; return Cache(n, factors); } } factors[n] = 1; return Cache(n, factors); }
private static IntMap Cache(int n, IntMap factors) { if (n < MaxFactors) factorCache[n] = factors; return factors; }
private static IEnumerable<int> GetDivisors(int n) { for (int i = 1; i * i <= n; i++) { if (n % i == 0) { yield return i; int div = n / i; if (div != i && div != n) yield return div; } } }
public sealed class Polynomial : IEnumerable<Term> { public Polynomial() { } public Polynomial(params Term[] terms) : this(terms.AsEnumerable()) { }
public Polynomial(IEnumerable<Term> terms) { Terms.AddRange(terms); Simplify(); }
private List<Term>? terms; private List<Term> Terms => terms ??= new List<Term>();
public int Count => terms?.Count ?? 0; public int Degree => Count == 0 ? -1 : Terms[0].Exponent; public int LeadingCoefficient => Count == 0 ? 0 : Terms[0].Coefficient;
public IEnumerator<Term> GetEnumerator() => Terms.GetEnumerator(); IEnumerator IEnumerable.GetEnumerator() => GetEnumerator();
public override string ToString() => Count == 0 ? "0" : string.Join(" + ", Terms).Replace("+ -", "- ");
public static Polynomial operator *(Polynomial p, Term t) => new Polynomial(from s in p select s * t); public static Polynomial operator +(Polynomial p, Polynomial q) => new Polynomial(p.Terms.Concat(q.Terms)); public static Polynomial operator -(Polynomial p, Polynomial q) => new Polynomial(p.Terms.Concat(q.Terms.Select(t => -t))); public static Polynomial operator *(Polynomial p, Polynomial q) => new Polynomial(from s in p from t in q select s * t); public static Polynomial operator /(Polynomial p, Polynomial q) => p.Divide(q).quotient;
public (Polynomial quotient, Polynomial remainder) Divide(Polynomial divisor) { if (Degree < 0) return (new Polynomial(), this); Polynomial quotient = new Polynomial(); Polynomial remainder = this; int lcv = divisor.LeadingCoefficient; int dv = divisor.Degree; while (remainder.Degree >= divisor.Degree) { int lcr = remainder.LeadingCoefficient; Term div = new Term(lcr / lcv, remainder.Degree - dv); quotient.Terms.Add(div); remainder += divisor * -div; } quotient.Simplify(); remainder.Simplify(); return (quotient, remainder); }
private void Simplify() { if (Count < 2) return; Terms.Sort((a, b) => -a.CompareTo(b)); for (int i = Terms.Count - 1; i > 0; i--) { Term s = Terms[i-1]; Term t = Terms[i]; if (t.Exponent == s.Exponent) { Terms[i-1] = new Term(s.Coefficient + t.Coefficient, s.Exponent); Terms.RemoveAt(i); } } Terms.RemoveAll(t => t.IsZero); }
} public readonly struct Term : IEquatable<Term>, IComparable<Term> { public Term(int coefficient, int exponent = 0) => (Coefficient, Exponent) = (coefficient, exponent);
public Term this[int exponent] => new Term(Coefficient, exponent); //Using x[exp] because x^exp has low precedence public int Coefficient { get; } public int Exponent { get; } public bool IsZero => Coefficient == 0;
public static Polynomial operator +(Term left, Term right) => new Polynomial(left, right); public static Polynomial operator -(Term left, Term right) => new Polynomial(left, -right); public static implicit operator Term(int coefficient) => new Term(coefficient); public static Term operator -(Term t) => new Term(-t.Coefficient, t.Exponent); public static Term operator *(Term left, Term right) => new Term(left.Coefficient * right.Coefficient, left.Exponent + right.Exponent);
public static bool operator ==(Term left, Term right) => left.Equals(right); public static bool operator !=(Term left, Term right) => !left.Equals(right); public static bool operator <(Term left, Term right) => left.CompareTo(right) < 0; public static bool operator >(Term left, Term right) => left.CompareTo(right) > 0; public static bool operator <=(Term left, Term right) => left.CompareTo(right) <= 0; public static bool operator >=(Term left, Term right) => left.CompareTo(right) >= 0;
public bool Equals(Term other) => Exponent == other.Exponent && Coefficient == other.Coefficient; public override bool Equals(object? obj) => obj is Term t && Equals(t); public override int GetHashCode() => Coefficient.GetHashCode() * 31 + Exponent.GetHashCode();
public int CompareTo(Term other) { int c = Exponent.CompareTo(other.Exponent); if (c != 0) return c; return Coefficient.CompareTo(other.Coefficient); }
public override string ToString() => (Coefficient, Exponent) switch { (0, _) => "0", (_, 0) => $"{Coefficient}", (1, 1) => "x", (-1, 1) => "-x", (_, 1) => $"{Coefficient}x", (1, _) => $"x^{Exponent}", (-1, _) => $"-x^{Exponent}", _ => $"{Coefficient}x^{Exponent}" }; }
}</lang>
- Output:
Task 1: Cyclotomic polynomials for n <= 30: CP[1] = x - 1 CP[2] = x + 1 CP[3] = x^2 + x + 1 CP[4] = x^2 + 1 CP[5] = x^4 + x^3 + x^2 + x + 1 CP[6] = x^2 - x + 1 CP[7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[8] = x^4 + 1 CP[9] = x^6 + x^3 + 1 CP[10] = x^4 - x^3 + x^2 - x + 1 CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[12] = x^4 - x^2 + 1 CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 CP[16] = x^8 + 1 CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[18] = x^6 - x^3 + 1 CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[20] = x^8 - x^6 + x^4 - x^2 + 1 CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[24] = x^8 - x^4 + 1 CP[25] = x^20 + x^15 + x^10 + x^5 + 1 CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[27] = x^18 + x^9 + 1 CP[28] = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 CP[29] = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[30] = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient: CP[1] has coefficient with magnitude = 1 CP[105] has coefficient with magnitude = 2 CP[385] has coefficient with magnitude = 3 CP[1365] has coefficient with magnitude = 4 CP[1785] has coefficient with magnitude = 5 CP[2805] has coefficient with magnitude = 6 CP[3135] has coefficient with magnitude = 7 CP[6545] has coefficient with magnitude = 8 CP[6545] has coefficient with magnitude = 9 CP[10465] has coefficient with magnitude = 10
Go
<lang go>package main
import (
"fmt" "log" "math" "sort" "strings"
)
const (
algo = 2 maxAllFactors = 100000
)
func iabs(i int) int {
if i < 0 { return -i } return i
}
type term struct{ coef, exp int }
func (t term) mul(t2 term) term {
return term{t.coef * t2.coef, t.exp + t2.exp}
}
func (t term) add(t2 term) term {
if t.exp != t2.exp { log.Fatal("exponents unequal in term.add method") } return term{t.coef + t2.coef, t.exp}
}
func (t term) negate() term { return term{-t.coef, t.exp} }
func (t term) String() string {
switch { case t.coef == 0: return "0" case t.exp == 0: return fmt.Sprintf("%d", t.coef) case t.coef == 1: if t.exp == 1 { return "x" } else { return fmt.Sprintf("x^%d", t.exp) } case t.exp == 1: return fmt.Sprintf("%dx", t.coef) } return fmt.Sprintf("%dx^%d", t.coef, t.exp)
}
type poly struct{ terms []term }
// pass coef, exp in pairs as parameters func newPoly(values ...int) poly {
le := len(values) if le == 0 { return poly{[]term{term{0, 0}}} } if le%2 != 0 { log.Fatalf("odd number of parameters (%d) passed to newPoly function", le) } var terms []term for i := 0; i < le; i += 2 { terms = append(terms, term{values[i], values[i+1]}) } p := poly{terms}.tidy() return p
}
func (p poly) hasCoefAbs(coef int) bool {
for _, t := range p.terms { if iabs(t.coef) == coef { return true } } return false
}
func (p poly) add(p2 poly) poly {
p3 := newPoly() le, le2 := len(p.terms), len(p2.terms) for le > 0 || le2 > 0 { if le == 0 { p3.terms = append(p3.terms, p2.terms[le2-1]) le2-- } else if le2 == 0 { p3.terms = append(p3.terms, p.terms[le-1]) le-- } else { t := p.terms[le-1] t2 := p2.terms[le2-1] if t.exp == t2.exp { t3 := t.add(t2) if t3.coef != 0 { p3.terms = append(p3.terms, t3) } le-- le2-- } else if t.exp < t2.exp { p3.terms = append(p3.terms, t) le-- } else { p3.terms = append(p3.terms, t2) le2-- } } } return p3.tidy()
}
func (p poly) addTerm(t term) poly {
q := newPoly() added := false for i := 0; i < len(p.terms); i++ { ct := p.terms[i] if ct.exp == t.exp { added = true if ct.coef+t.coef != 0 { q.terms = append(q.terms, ct.add(t)) } } else { q.terms = append(q.terms, ct) } } if !added { q.terms = append(q.terms, t) } return q.tidy()
}
func (p poly) mulTerm(t term) poly {
q := newPoly() for i := 0; i < len(p.terms); i++ { ct := p.terms[i] q.terms = append(q.terms, ct.mul(t)) } return q.tidy()
}
func (p poly) div(v poly) poly {
q := newPoly() lcv := v.leadingCoef() dv := v.degree() for p.degree() >= v.degree() { lcp := p.leadingCoef() s := lcp / lcv t := term{s, p.degree() - dv} q = q.addTerm(t) p = p.add(v.mulTerm(t.negate())) } return q.tidy()
}
func (p poly) leadingCoef() int {
return p.terms[0].coef
}
func (p poly) degree() int {
return p.terms[0].exp
}
func (p poly) String() string {
var sb strings.Builder first := true for _, t := range p.terms { if first { sb.WriteString(t.String()) first = false } else { sb.WriteString(" ") if t.coef > 0 { sb.WriteString("+ ") sb.WriteString(t.String()) } else { sb.WriteString("- ") sb.WriteString(t.negate().String()) } } } return sb.String()
}
// in place descending sort by term.exp func (p poly) sortTerms() {
sort.Slice(p.terms, func(i, j int) bool { return p.terms[i].exp > p.terms[j].exp })
}
// sort terms and remove any unnecesary zero terms func (p poly) tidy() poly {
p.sortTerms() if p.degree() == 0 { return p } for i := len(p.terms) - 1; i >= 0; i-- { if p.terms[i].coef == 0 { copy(p.terms[i:], p.terms[i+1:]) p.terms[len(p.terms)-1] = term{0, 0} p.terms = p.terms[:len(p.terms)-1] } } if len(p.terms) == 0 { p.terms = append(p.terms, term{0, 0}) } return p
}
func getDivisors(n int) []int {
var divs []int sqrt := int(math.Sqrt(float64(n))) for i := 1; i <= sqrt; i++ { if n%i == 0 { divs = append(divs, i) d := n / i if d != i && d != n { divs = append(divs, d) } } } return divs
}
var (
computed = make(map[int]poly) allFactors = make(map[int]map[int]int)
)
func init() {
f := map[int]int{2: 1} allFactors[2] = f
}
func getFactors(n int) map[int]int {
if f, ok := allFactors[n]; ok { return f } factors := make(map[int]int) if n%2 == 0 { factorsDivTwo := getFactors(n / 2) for k, v := range factorsDivTwo { factors[k] = v } factors[2]++ if n < maxAllFactors { allFactors[n] = factors } return factors } prime := true sqrt := int(math.Sqrt(float64(n))) for i := 3; i <= sqrt; i += 2 { if n%i == 0 { prime = false for k, v := range getFactors(n / i) { factors[k] = v } factors[i]++ if n < maxAllFactors { allFactors[n] = factors } return factors } } if prime { factors[n] = 1 if n < maxAllFactors { allFactors[n] = factors } } return factors
}
func cycloPoly(n int) poly {
if p, ok := computed[n]; ok { return p } if n == 1 { // polynomial: x - 1 p := newPoly(1, 1, -1, 0) computed[1] = p return p } factors := getFactors(n) cyclo := newPoly() if _, ok := factors[n]; ok { // n is prime for i := 0; i < n; i++ { cyclo.terms = append(cyclo.terms, term{1, i}) } } else if len(factors) == 2 && factors[2] == 1 && factors[n/2] == 1 { // n == 2p prime := n / 2 coef := -1 for i := 0; i < prime; i++ { coef *= -1 cyclo.terms = append(cyclo.terms, term{coef, i}) } } else if len(factors) == 1 { if h, ok := factors[2]; ok { // n == 2^h cyclo.terms = append(cyclo.terms, term{1, 1 << (h - 1)}, term{1, 0}) } else if _, ok := factors[n]; !ok { // n == p ^ k p := 0 for prime := range factors { p = prime } k := factors[p] for i := 0; i < p; i++ { pk := int(math.Pow(float64(p), float64(k-1))) cyclo.terms = append(cyclo.terms, term{1, i * pk}) } } } else if len(factors) == 2 && factors[2] != 0 { // n = 2^h * p^k p := 0 for prime := range factors { if prime != 2 { p = prime } } coef := -1 twoExp := 1 << (factors[2] - 1) k := factors[p] for i := 0; i < p; i++ { coef *= -1 pk := int(math.Pow(float64(p), float64(k-1))) cyclo.terms = append(cyclo.terms, term{coef, i * twoExp * pk}) } } else if factors[2] != 0 && ((n/2)%2 == 1) && (n/2) > 1 { // CP(2m)[x] == CP(-m)[x], n odd integer > 1 cycloDiv2 := cycloPoly(n / 2) for _, t := range cycloDiv2.terms { t2 := t if t.exp%2 != 0 { t2 = t.negate() } cyclo.terms = append(cyclo.terms, t2) } } else if algo == 0 { // slow - uses basic definition divs := getDivisors(n) // polynomial: x^n - 1 cyclo = newPoly(1, n, -1, 0) for _, i := range divs { p := cycloPoly(i) cyclo = cyclo.div(p) } } else if algo == 1 { // faster - remove max divisor (and all divisors of max divisor) // only one divide for all divisors of max divisor divs := getDivisors(n) maxDiv := math.MinInt32 for _, d := range divs { if d > maxDiv { maxDiv = d } } var divsExceptMax []int for _, d := range divs { if maxDiv%d != 0 { divsExceptMax = append(divsExceptMax, d) } } // polynomial: ( x^n - 1 ) / ( x^m - 1 ), where m is the max divisor cyclo = newPoly(1, n, -1, 0) cyclo = cyclo.div(newPoly(1, maxDiv, -1, 0)) for _, i := range divsExceptMax { p := cycloPoly(i) cyclo = cyclo.div(p) } } else if algo == 2 { // fastest // let p, q be primes such that p does not divide n, and q divides n // then CP(np)[x] = CP(n)[x^p] / CP(n)[x] m := 1 cyclo = cycloPoly(m) var primes []int for prime := range factors { primes = append(primes, prime) } sort.Ints(primes) for _, prime := range primes { // CP(m)[x] cycloM := cyclo // compute CP(m)[x^p] var terms []term for _, t := range cycloM.terms { terms = append(terms, term{t.coef, t.exp * prime}) } cyclo = newPoly() cyclo.terms = append(cyclo.terms, terms...) cyclo = cyclo.tidy() cyclo = cyclo.div(cycloM) m *= prime } // now, m is the largest square free divisor of n s := n / m // Compute CP(n)[x] = CP(m)[x^s] var terms []term for _, t := range cyclo.terms { terms = append(terms, term{t.coef, t.exp * s}) } cyclo = newPoly() cyclo.terms = append(cyclo.terms, terms...) } else { log.Fatal("invalid algorithm") } cyclo = cyclo.tidy() computed[n] = cyclo return cyclo
}
func main() {
fmt.Println("Task 1: cyclotomic polynomials for n <= 30:") for i := 1; i <= 30; i++ { p := cycloPoly(i) fmt.Printf("CP[%2d] = %s\n", i, p) }
fmt.Println("\nTask 2: Smallest cyclotomic polynomial with n or -n as a coefficient:") n := 0 for i := 1; i <= 10; i++ { for { n++ cyclo := cycloPoly(n) if cyclo.hasCoefAbs(i) { fmt.Printf("CP[%d] has coefficient with magnitude = %d\n", n, i) n-- break } } }
}</lang>
- Output:
Task 1: cyclotomic polynomials for n <= 30: CP[ 1] = x - 1 CP[ 2] = x + 1 CP[ 3] = x^2 + x + 1 CP[ 4] = x^2 + 1 CP[ 5] = x^4 + x^3 + x^2 + x + 1 CP[ 6] = x^2 - x + 1 CP[ 7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[ 8] = x^4 + 1 CP[ 9] = x^6 + x^3 + 1 CP[10] = x^4 - x^3 + x^2 - x + 1 CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[12] = x^4 - x^2 + 1 CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 CP[16] = x^8 + 1 CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[18] = x^6 - x^3 + 1 CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[20] = x^8 - x^6 + x^4 - x^2 + 1 CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[24] = x^8 - x^4 + 1 CP[25] = x^20 + x^15 + x^10 + x^5 + 1 CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[27] = x^18 + x^9 + 1 CP[28] = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 CP[29] = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[30] = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient: CP[1] has coefficient with magnitude = 1 CP[105] has coefficient with magnitude = 2 CP[385] has coefficient with magnitude = 3 CP[1365] has coefficient with magnitude = 4 CP[1785] has coefficient with magnitude = 5 CP[2805] has coefficient with magnitude = 6 CP[3135] has coefficient with magnitude = 7 CP[6545] has coefficient with magnitude = 8 CP[6545] has coefficient with magnitude = 9 CP[10465] has coefficient with magnitude = 10
Java
<lang java> import java.util.ArrayList; import java.util.Collections; import java.util.Comparator; import java.util.HashMap; import java.util.List; import java.util.Map; import java.util.TreeMap;
public class CyclotomicPolynomial {
@SuppressWarnings("unused") private static int divisions = 0; private static int algorithm = 2; public static void main(String[] args) throws Exception { System.out.println("Task 1: cyclotomic polynomials for n <= 30:"); for ( int i = 1 ; i <= 30 ; i++ ) { Polynomial p = cyclotomicPolynomial(i); System.out.printf("CP[%d] = %s%n", i, p); } System.out.println("Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:"); int n = 0; for ( int i = 1 ; i <= 10 ; i++ ) { while ( true ) { n++; Polynomial cyclo = cyclotomicPolynomial(n); if ( cyclo.hasCoefficientAbs(i) ) { System.out.printf("CP[%d] has coefficient with magnitude = %d%n", n, i); n--; break; } } } }
private static final Map<Integer, Polynomial> COMPUTED = new HashMap<>(); private static Polynomial cyclotomicPolynomial(int n) { if ( COMPUTED.containsKey(n) ) { return COMPUTED.get(n); } //System.out.println("COMPUTE: n = " + n); if ( n == 1 ) { // Polynomial: x - 1 Polynomial p = new Polynomial(1, 1, -1, 0); COMPUTED.put(1, p); return p; }
Map<Integer,Integer> factors = getFactors(n); if ( factors.containsKey(n) ) { // n prime List<Term> termList = new ArrayList<>(); for ( int index = 0 ; index < n ; index++ ) { termList.add(new Term(1, index)); } Polynomial cyclo = new Polynomial(termList); COMPUTED.put(n, cyclo); return cyclo; } else if ( factors.size() == 2 && factors.containsKey(2) && factors.get(2) == 1 && factors.containsKey(n/2) && factors.get(n/2) == 1 ) { // n = 2p int prime = n/2; List<Term> termList = new ArrayList<>(); int coeff = -1; for ( int index = 0 ; index < prime ; index++ ) { coeff *= -1; termList.add(new Term(coeff, index)); }
Polynomial cyclo = new Polynomial(termList); COMPUTED.put(n, cyclo); return cyclo; } else if ( factors.size() == 1 && factors.containsKey(2) ) { // n = 2^h int h = factors.get(2); List<Term> termList = new ArrayList<>(); termList.add(new Term(1, (int) Math.pow(2, h-1))); termList.add(new Term(1, 0)); Polynomial cyclo = new Polynomial(termList); COMPUTED.put(n, cyclo); return cyclo; } else if ( factors.size() == 1 && ! factors.containsKey(n) ) { // n = p^k int p = 0; for ( int prime : factors.keySet() ) { p = prime; } int k = factors.get(p); List<Term> termList = new ArrayList<>(); for ( int index = 0 ; index < p ; index++ ) { termList.add(new Term(1, index * (int) Math.pow(p, k-1))); }
Polynomial cyclo = new Polynomial(termList); COMPUTED.put(n, cyclo); return cyclo; } else if ( factors.size() == 2 && factors.containsKey(2) ) { // n = 2^h * p^k int p = 0; for ( int prime : factors.keySet() ) { if ( prime != 2 ) { p = prime; } } List<Term> termList = new ArrayList<>(); int coeff = -1; int twoExp = (int) Math.pow(2, factors.get(2)-1); int k = factors.get(p); for ( int index = 0 ; index < p ; index++ ) { coeff *= -1; termList.add(new Term(coeff, index * twoExp * (int) Math.pow(p, k-1))); }
Polynomial cyclo = new Polynomial(termList); COMPUTED.put(n, cyclo); return cyclo; } else if ( factors.containsKey(2) && ((n/2) % 2 == 1) && (n/2) > 1 ) { // CP(2m)[x] = CP(-m)[x], n odd integer > 1 Polynomial cycloDiv2 = cyclotomicPolynomial(n/2); List<Term> termList = new ArrayList<>(); for ( Term term : cycloDiv2.polynomialTerms ) { termList.add(term.exponent % 2 == 0 ? term : term.negate()); } Polynomial cyclo = new Polynomial(termList); COMPUTED.put(n, cyclo); return cyclo; } // General Case if ( algorithm == 0 ) { // Slow - uses basic definition. List<Integer> divisors = getDivisors(n); // Polynomial: ( x^n - 1 ) Polynomial cyclo = new Polynomial(1, n, -1, 0); for ( int i : divisors ) { Polynomial p = cyclotomicPolynomial(i); cyclo = cyclo.divide(p); } COMPUTED.put(n, cyclo); return cyclo; } else if ( algorithm == 1 ) { // Faster. Remove Max divisor (and all divisors of max divisor) - only one divide for all divisors of Max Divisor List<Integer> divisors = getDivisors(n); int maxDivisor = Integer.MIN_VALUE; for ( int div : divisors ) { maxDivisor = Math.max(maxDivisor, div); } List<Integer> divisorsExceptMax = new ArrayList<Integer>(); for ( int div : divisors ) { if ( maxDivisor % div != 0 ) { divisorsExceptMax.add(div); } }
// Polynomial: ( x^n - 1 ) / ( x^m - 1 ), where m is the max divisor Polynomial cyclo = new Polynomial(1, n, -1, 0).divide(new Polynomial(1, maxDivisor, -1, 0)); for ( int i : divisorsExceptMax ) { Polynomial p = cyclotomicPolynomial(i); cyclo = cyclo.divide(p); }
COMPUTED.put(n, cyclo);
return cyclo; } else if ( algorithm == 2 ) { // Fastest // Let p ; q be primes such that p does not divide n, and q q divides n. // Then CP(np)[x] = CP(n)[x^p] / CP(n)[x] int m = 1; Polynomial cyclo = cyclotomicPolynomial(m); List<Integer> primes = new ArrayList<>(factors.keySet()); Collections.sort(primes); for ( int prime : primes ) { // CP(m)[x] Polynomial cycloM = cyclo; // Compute CP(m)[x^p]. List<Term> termList = new ArrayList<>(); for ( Term t : cycloM.polynomialTerms ) { termList.add(new Term(t.coefficient, t.exponent * prime)); } cyclo = new Polynomial(termList).divide(cycloM); m = m * prime; } // Now, m is the largest square free divisor of n int s = n / m; // Compute CP(n)[x] = CP(m)[x^s] List<Term> termList = new ArrayList<>(); for ( Term t : cyclo.polynomialTerms ) { termList.add(new Term(t.coefficient, t.exponent * s)); } cyclo = new Polynomial(termList); COMPUTED.put(n, cyclo);
return cyclo; } else { throw new RuntimeException("ERROR 103: Invalid algorithm."); } } private static final List<Integer> getDivisors(int number) { List<Integer> divisors = new ArrayList<Integer>(); long sqrt = (long) Math.sqrt(number); for ( int i = 1 ; i <= sqrt ; i++ ) { if ( number % i == 0 ) { divisors.add(i); int div = number / i; if ( div != i && div != number ) { divisors.add(div); } } } return divisors; }
private static final Map<Integer,Map<Integer,Integer>> allFactors = new TreeMap<Integer,Map<Integer,Integer>>(); static { Map<Integer,Integer> factors = new TreeMap<Integer,Integer>(); factors.put(2, 1); allFactors.put(2, factors); }
public static Integer MAX_ALL_FACTORS = 100000;
public static final Map<Integer,Integer> getFactors(Integer number) { if ( allFactors.containsKey(number) ) { return allFactors.get(number); } Map<Integer,Integer> factors = new TreeMap<Integer,Integer>(); if ( number % 2 == 0 ) { Map<Integer,Integer> factorsdDivTwo = getFactors(number/2); factors.putAll(factorsdDivTwo); factors.merge(2, 1, (v1, v2) -> v1 + v2); if ( number < MAX_ALL_FACTORS ) allFactors.put(number, factors); return factors; } boolean prime = true; long sqrt = (long) Math.sqrt(number); for ( int i = 3 ; i <= sqrt ; i += 2 ) { if ( number % i == 0 ) { prime = false; factors.putAll(getFactors(number/i)); factors.merge(i, 1, (v1, v2) -> v1 + v2); if ( number < MAX_ALL_FACTORS ) allFactors.put(number, factors); return factors; } } if ( prime ) { factors.put(number, 1); if ( number < MAX_ALL_FACTORS ) allFactors.put(number, factors); } return factors; } private static final class Polynomial {
private List<Term> polynomialTerms; // Format - coeff, exp, coeff, exp, (repeating in pairs) . . . public Polynomial(int ... values) { if ( values.length % 2 != 0 ) { throw new IllegalArgumentException("ERROR 102: Polynomial constructor. Length must be even. Length = " + values.length); } polynomialTerms = new ArrayList<>(); for ( int i = 0 ; i < values.length ; i += 2 ) { Term t = new Term(values[i], values[i+1]); polynomialTerms.add(t); } Collections.sort(polynomialTerms, new TermSorter()); } public Polynomial() { // zero polynomialTerms = new ArrayList<>(); polynomialTerms.add(new Term(0,0)); } private boolean hasCoefficientAbs(int coeff) { for ( Term term : polynomialTerms ) { if ( Math.abs(term.coefficient) == coeff ) { return true; } } return false; } private Polynomial(List<Term> termList) { if ( termList.size() == 0 ) { // zero termList.add(new Term(0,0)); } else { // Remove zero terms if needed for ( int i = 0 ; i < termList.size() ; i++ ) { if ( termList.get(i).coefficient == 0 ) { termList.remove(i); } } } if ( termList.size() == 0 ) { // zero termList.add(new Term(0,0)); } polynomialTerms = termList; Collections.sort(polynomialTerms, new TermSorter()); } public Polynomial divide(Polynomial v) { //long start = System.currentTimeMillis(); divisions++; Polynomial q = new Polynomial(); Polynomial r = this; long lcv = v.leadingCoefficient(); long dv = v.degree(); while ( r.degree() >= v.degree() ) { long lcr = r.leadingCoefficient(); long s = lcr / lcv; // Integer division Term term = new Term(s, r.degree() - dv); q = q.add(term); r = r.add(v.multiply(term.negate())); } //long end = System.currentTimeMillis(); //System.out.printf("Divide: Elapsed = %d, Degree 1 = %d, Degree 2 = %d%n", (end-start), this.degree(), v.degree()); return q; }
public Polynomial add(Polynomial polynomial) { List<Term> termList = new ArrayList<>(); int thisCount = polynomialTerms.size(); int polyCount = polynomial.polynomialTerms.size(); while ( thisCount > 0 || polyCount > 0 ) { Term thisTerm = thisCount == 0 ? null : polynomialTerms.get(thisCount-1); Term polyTerm = polyCount == 0 ? null : polynomial.polynomialTerms.get(polyCount-1); if ( thisTerm == null ) { termList.add(polyTerm.clone()); polyCount--; } else if (polyTerm == null ) { termList.add(thisTerm.clone()); thisCount--; } else if ( thisTerm.degree() == polyTerm.degree() ) { Term t = thisTerm.add(polyTerm); if ( t.coefficient != 0 ) { termList.add(t); } thisCount--; polyCount--; } else if ( thisTerm.degree() < polyTerm.degree() ) { termList.add(thisTerm.clone()); thisCount--; } else { termList.add(polyTerm.clone()); polyCount--; } } return new Polynomial(termList); }
public Polynomial add(Term term) { List<Term> termList = new ArrayList<>(); boolean added = false; for ( int index = 0 ; index < polynomialTerms.size() ; index++ ) { Term currentTerm = polynomialTerms.get(index); if ( currentTerm.exponent == term.exponent ) { added = true; if ( currentTerm.coefficient + term.coefficient != 0 ) { termList.add(currentTerm.add(term)); } } else { termList.add(currentTerm.clone()); } } if ( ! added ) { termList.add(term.clone()); } return new Polynomial(termList); }
public Polynomial multiply(Term term) { List<Term> termList = new ArrayList<>(); for ( int index = 0 ; index < polynomialTerms.size() ; index++ ) { Term currentTerm = polynomialTerms.get(index); termList.add(currentTerm.clone().multiply(term)); } return new Polynomial(termList); }
public Polynomial clone() { List<Term> clone = new ArrayList<>(); for ( Term t : polynomialTerms ) { clone.add(new Term(t.coefficient, t.exponent)); } return new Polynomial(clone); }
public long leadingCoefficient() {
// long coefficient = 0; // long degree = Integer.MIN_VALUE; // for ( Term t : polynomialTerms ) { // if ( t.degree() > degree ) { // coefficient = t.coefficient; // degree = t.degree(); // } // }
return polynomialTerms.get(0).coefficient; }
public long degree() {
// long degree = Integer.MIN_VALUE; // for ( Term t : polynomialTerms ) { // if ( t.degree() > degree ) { // degree = t.degree(); // } // }
return polynomialTerms.get(0).exponent; } @Override public String toString() { StringBuilder sb = new StringBuilder(); boolean first = true; for ( Term term : polynomialTerms ) { if ( first ) { sb.append(term); first = false; } else { sb.append(" "); if ( term.coefficient > 0 ) { sb.append("+ "); sb.append(term); } else { sb.append("- "); sb.append(term.negate()); } } } return sb.toString(); } } private static final class TermSorter implements Comparator<Term> { @Override public int compare(Term o1, Term o2) { return (int) (o2.exponent - o1.exponent); } } // Note: Cyclotomic Polynomials have small coefficients. Not appropriate for general polynomial usage. private static final class Term { long coefficient; long exponent; public Term(long c, long e) { coefficient = c; exponent = e; } public Term clone() { return new Term(coefficient, exponent); } public Term multiply(Term term) { return new Term(coefficient * term.coefficient, exponent + term.exponent); } public Term add(Term term) { if ( exponent != term.exponent ) { throw new RuntimeException("ERROR 102: Exponents not equal."); } return new Term(coefficient + term.coefficient, exponent); }
public Term negate() { return new Term(-coefficient, exponent); } public long degree() { return exponent; } @Override public String toString() { if ( coefficient == 0 ) { return "0"; } if ( exponent == 0 ) { return "" + coefficient; } if ( coefficient == 1 ) { if ( exponent == 1 ) { return "x"; } else { return "x^" + exponent; } } if ( exponent == 1 ) { return coefficient + "x"; } return coefficient + "x^" + exponent; } }
} </lang>
- Output:
Task 1: cyclotomic polynomials for n <= 30: CP[1] = x - 1 CP[2] = x + 1 CP[3] = x^2 + x + 1 CP[4] = x^2 + 1 CP[5] = x^4 + x^3 + x^2 + x + 1 CP[6] = x^2 - x + 1 CP[7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[8] = x^4 + 1 CP[9] = x^6 + x^3 + 1 CP[10] = x^4 - x^3 + x^2 - x + 1 CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[12] = x^4 - x^2 + 1 CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 CP[16] = x^8 + 1 CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[18] = x^6 - x^3 + 1 CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[20] = x^8 - x^6 + x^4 - x^2 + 1 CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[24] = x^8 - x^4 + 1 CP[25] = x^20 + x^15 + x^10 + x^5 + 1 CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[27] = x^18 + x^9 + 1 CP[28] = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 CP[29] = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[30] = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient: CP[1] has coefficient with magnitude = 1 CP[105] has coefficient with magnitude = 2 CP[385] has coefficient with magnitude = 3 CP[1365] has coefficient with magnitude = 4 CP[1785] has coefficient with magnitude = 5 CP[2805] has coefficient with magnitude = 6 CP[3135] has coefficient with magnitude = 7 CP[6545] has coefficient with magnitude = 8 CP[6545] has coefficient with magnitude = 9 CP[10465] has coefficient with magnitude = 10
Julia
<lang julia>using Primes, Polynomials
- memoize cache for recursive calls
const cyclotomics = Dict([1 => Poly([big"-1", big"1"]), 2 => Poly([big"1", big"1"])])
- get all integer divisors of an integer, except itself
function divisors(n::Integer)
f = [one(n)] for (p, e) in factor(n) f = reduce(vcat, [f * p^j for j in 1:e], init=f) end return resize!(f, length(f) - 1)
end
"""
cyclotomic(n::Integer)
Calculate the n -th cyclotomic polynomial. See wikipedia article at bottom of section /wiki/Cyclotomic_polynomial#Fundamental_tools The algorithm is reliable but slow for large n > 1000. """ function cyclotomic(n::Integer)
if haskey(cyclotomics, n) c = cyclotomics[n] elseif isprime(n) c = Poly(ones(BigInt, n)) cyclotomics[n] = c else # recursive formula seen in wikipedia article c = Poly([big"-1"; zeros(BigInt, n - 1); big"1"]) for d in divisors(n) c ÷= cyclotomic(d) end cyclotomics[n] = c end return c
end
println("First 30 cyclotomic polynomials:") for i in 1:30
println(rpad("$i: ", 5), cyclotomic(BigInt(i)))
end
const dig = zeros(BigInt, 10) for i in 1:1000000
if all(x -> x != 0, dig) break end for coef in coeffs(cyclotomic(i)) x = abs(coef) if 0 < x < 11 && dig[Int(x)] == 0 dig[Int(x)] = coef < 0 ? -i : i end end
end for (i, n) in enumerate(dig)
println("The cyclotomic polynomial Φ(", abs(n), ") has a coefficient that is ", n < 0 ? -i : i)
end
</lang>
- Output:
First 30 cyclotomic polynomials: 1: Poly(-1 + x) 2: Poly(1 + x) 3: Poly(1 + x + x^2) 4: Poly(1.0 + 1.0*x^2) 5: Poly(1 + x + x^2 + x^3 + x^4) 6: Poly(1.0 - 1.0*x + 1.0*x^2) 7: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6) 8: Poly(1.0 + 1.0*x^4) 9: Poly(1.0 + 1.0*x^3 + 1.0*x^6) 10: Poly(1.0 - 1.0*x + 1.0*x^2 - 1.0*x^3 + 1.0*x^4) 11: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10) 12: Poly(1.0 - 1.0*x^2 + 1.0*x^4) 13: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12) 14: Poly(1.0 - 1.0*x + 1.0*x^2 - 1.0*x^3 + 1.0*x^4 - 1.0*x^5 + 1.0*x^6) 15: Poly(1.0 - 1.0*x + 1.0*x^3 - 1.0*x^4 + 1.0*x^5 - 1.0*x^7 + 1.0*x^8) 16: Poly(1.0 + 1.0*x^8) 17: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12 + x^13 + x^14 + x^15 + x^16) 18: Poly(1.0 - 1.0*x^3 + 1.0*x^6) 19: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12 + x^13 + x^14 + x^15 + x^16 + x^17 + x^18) 20: Poly(1.0 - 1.0*x^2 + 1.0*x^4 - 1.0*x^6 + 1.0*x^8) 21: Poly(1.0 - 1.0*x + 1.0*x^3 - 1.0*x^4 + 1.0*x^6 - 1.0*x^8 + 1.0*x^9 - 1.0*x^11 + 1.0*x^12) 22: Poly(1.0 - 1.0*x + 1.0*x^2 - 1.0*x^3 + 1.0*x^4 - 1.0*x^5 + 1.0*x^6 - 1.0*x^7 + 1.0*x^8 - 1.0*x^9 + 1.0*x^10) 23: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12 + x^13 + x^14 + x^15 + x^16 + x^17 + x^18 + x^19 + x^20 + x^21 + x^22) 24: Poly(1.0 - 1.0*x^4 + 1.0*x^8) 25: Poly(1.0 + 1.0*x^5 + 1.0*x^10 + 1.0*x^15 + 1.0*x^20) 26: Poly(1.0 - 1.0*x + 1.0*x^2 - 1.0*x^3 + 1.0*x^4 - 1.0*x^5 + 1.0*x^6 - 1.0*x^7 + 1.0*x^8 - 1.0*x^9 + 1.0*x^10 - 1.0*x^11 + 1.0*x^12) 27: Poly(1.0 + 1.0*x^9 + 1.0*x^18) 28: Poly(1.0 - 1.0*x^2 + 1.0*x^4 - 1.0*x^6 + 1.0*x^8 - 1.0*x^10 + 1.0*x^12) 29: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12 + x^13 + x^14 + x^15 + x^16 + x^17 + x^18 + x^19 + x^20 + x^21 + x^22 + x^23 + x^24 + x^25 + x^26 + x^27 + x^28) 30: Poly(1.0 + 1.0*x - 1.0*x^3 - 1.0*x^4 - 1.0*x^5 + 1.0*x^7 + 1.0*x^8) The cyclotomic polynomial Φ(1) has a coefficient that is -1 The cyclotomic polynomial Φ(105) has a coefficient that is -2 The cyclotomic polynomial Φ(385) has a coefficient that is -3 The cyclotomic polynomial Φ(1365) has a coefficient that is -4 The cyclotomic polynomial Φ(1785) has a coefficient that is 5 The cyclotomic polynomial Φ(2805) has a coefficient that is -6 The cyclotomic polynomial Φ(3135) has a coefficient that is 7 The cyclotomic polynomial Φ(6545) has a coefficient that is -8 The cyclotomic polynomial Φ(6545) has a coefficient that is 9 The cyclotomic polynomial Φ(10465) has a coefficient that is 10
Kotlin
<lang scala>import java.util.TreeMap import kotlin.math.abs import kotlin.math.pow import kotlin.math.sqrt
private const val algorithm = 2
fun main() {
println("Task 1: cyclotomic polynomials for n <= 30:") for (i in 1..30) { val p = cyclotomicPolynomial(i) println("CP[$i] = $p") } println()
println("Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:") var n = 0 for (i in 1..10) { while (true) { n++ val cyclo = cyclotomicPolynomial(n) if (cyclo!!.hasCoefficientAbs(i)) { println("CP[$n] has coefficient with magnitude = $i") n-- break } } }
}
private val COMPUTED: MutableMap<Int, Polynomial> = HashMap() private fun cyclotomicPolynomial(n: Int): Polynomial? {
if (COMPUTED.containsKey(n)) { return COMPUTED[n] } if (n == 1) { // Polynomial: x - 1 val p = Polynomial(1, 1, -1, 0) COMPUTED[1] = p return p } val factors = getFactors(n) if (factors.containsKey(n)) { // n prime val termList: MutableList<Term> = ArrayList() for (index in 0 until n) { termList.add(Term(1, index.toLong())) } val cyclo = Polynomial(termList) COMPUTED[n] = cyclo return cyclo } else if (factors.size == 2 && factors.containsKey(2) && factors[2] == 1 && factors.containsKey(n / 2) && factors[n / 2] == 1) { // n = 2p val prime = n / 2 val termList: MutableList<Term> = ArrayList() var coeff = -1 for (index in 0 until prime) { coeff *= -1 termList.add(Term(coeff.toLong(), index.toLong())) } val cyclo = Polynomial(termList) COMPUTED[n] = cyclo return cyclo } else if (factors.size == 1 && factors.containsKey(2)) { // n = 2^h val h = factors[2]!! val termList: MutableList<Term> = ArrayList() termList.add(Term(1, 2.0.pow((h - 1).toDouble()).toLong())) termList.add(Term(1, 0)) val cyclo = Polynomial(termList) COMPUTED[n] = cyclo return cyclo } else if (factors.size == 1 && !factors.containsKey(n)) { // n = p^k var p = 0 for (prime in factors.keys) { p = prime } val k = factors[p]!! val termList: MutableList<Term> = ArrayList() for (index in 0 until p) { termList.add(Term(1, (index * p.toDouble().pow(k - 1.toDouble()).toInt()).toLong())) } val cyclo = Polynomial(termList) COMPUTED[n] = cyclo return cyclo } else if (factors.size == 2 && factors.containsKey(2)) { // n = 2^h * p^k var p = 0 for (prime in factors.keys) { if (prime != 2) { p = prime } } val termList: MutableList<Term> = ArrayList() var coeff = -1 val twoExp = 2.0.pow((factors[2]!!) - 1.toDouble()).toInt() val k = factors[p]!! for (index in 0 until p) { coeff *= -1 termList.add(Term(coeff.toLong(), (index * twoExp * p.toDouble().pow(k - 1.toDouble()).toInt()).toLong())) } val cyclo = Polynomial(termList) COMPUTED[n] = cyclo return cyclo } else if (factors.containsKey(2) && n / 2 % 2 == 1 && n / 2 > 1) { // CP(2m)[x] = CP(-m)[x], n odd integer > 1 val cycloDiv2 = cyclotomicPolynomial(n / 2) val termList: MutableList<Term> = ArrayList() for (term in cycloDiv2!!.polynomialTerms) { termList.add(if (term.exponent % 2 == 0L) term else term.negate()) } val cyclo = Polynomial(termList) COMPUTED[n] = cyclo return cyclo }
// General Case return when (algorithm) { 0 -> { // Slow - uses basic definition. val divisors = getDivisors(n) // Polynomial: ( x^n - 1 ) var cyclo = Polynomial(1, n, -1, 0) for (i in divisors) { val p = cyclotomicPolynomial(i) cyclo = cyclo.divide(p) } COMPUTED[n] = cyclo cyclo } 1 -> { // Faster. Remove Max divisor (and all divisors of max divisor) - only one divide for all divisors of Max Divisor val divisors = getDivisors(n) var maxDivisor = Int.MIN_VALUE for (div in divisors) { maxDivisor = maxDivisor.coerceAtLeast(div) } val divisorsExceptMax: MutableList<Int> = ArrayList() for (div in divisors) { if (maxDivisor % div != 0) { divisorsExceptMax.add(div) } }
// Polynomial: ( x^n - 1 ) / ( x^m - 1 ), where m is the max divisor var cyclo = Polynomial(1, n, -1, 0).divide(Polynomial(1, maxDivisor, -1, 0)) for (i in divisorsExceptMax) { val p = cyclotomicPolynomial(i) cyclo = cyclo.divide(p) } COMPUTED[n] = cyclo cyclo } 2 -> { // Fastest // Let p ; q be primes such that p does not divide n, and q q divides n. // Then CP(np)[x] = CP(n)[x^p] / CP(n)[x] var m = 1 var cyclo = cyclotomicPolynomial(m) val primes = factors.keys.toMutableList() primes.sort() for (prime in primes) { // CP(m)[x] val cycloM = cyclo // Compute CP(m)[x^p]. val termList: MutableList<Term> = ArrayList() for (t in cycloM!!.polynomialTerms) { termList.add(Term(t.coefficient, t.exponent * prime)) } cyclo = Polynomial(termList).divide(cycloM) m *= prime } // Now, m is the largest square free divisor of n val s = n / m // Compute CP(n)[x] = CP(m)[x^s] val termList: MutableList<Term> = ArrayList() for (t in cyclo!!.polynomialTerms) { termList.add(Term(t.coefficient, t.exponent * s)) } cyclo = Polynomial(termList) COMPUTED[n] = cyclo cyclo } else -> { throw RuntimeException("ERROR 103: Invalid algorithm.") } }
}
private fun getDivisors(number: Int): List<Int> {
val divisors: MutableList<Int> = ArrayList() val sqrt = sqrt(number.toDouble()).toLong() for (i in 1..sqrt) { if (number % i == 0L) { divisors.add(i.toInt()) val div = (number / i).toInt() if (div.toLong() != i && div != number) { divisors.add(div) } } } return divisors
}
private fun crutch(): MutableMap<Int, Map<Int, Int>> {
val allFactors: MutableMap<Int, Map<Int, Int>> = TreeMap()
val factors: MutableMap<Int, Int> = TreeMap() factors[2] = 1
allFactors[2] = factors return allFactors
}
private val allFactors = crutch()
var MAX_ALL_FACTORS = 100000
fun getFactors(number: Int): Map<Int, Int> {
if (allFactors.containsKey(number)) { return allFactors[number]!! } val factors: MutableMap<Int, Int> = TreeMap() if (number % 2 == 0) { val factorsDivTwo = getFactors(number / 2) factors.putAll(factorsDivTwo) factors.merge(2, 1) { a: Int?, b: Int? -> Integer.sum(a!!, b!!) } if (number < MAX_ALL_FACTORS) allFactors[number] = factors return factors } val sqrt = sqrt(number.toDouble()).toLong() var i = 3 while (i <= sqrt) { if (number % i == 0) { factors.putAll(getFactors(number / i)) factors.merge(i, 1) { a: Int?, b: Int? -> Integer.sum(a!!, b!!) } if (number < MAX_ALL_FACTORS) { allFactors[number] = factors } return factors } i += 2 } factors[number] = 1 if (number < MAX_ALL_FACTORS) { allFactors[number] = factors } return factors
}
private class Polynomial {
val polynomialTerms: MutableList<Term>
// Format - coeff, exp, coeff, exp, (repeating in pairs) . . . constructor(vararg values: Int) { require(values.size % 2 == 0) { "ERROR 102: Polynomial constructor. Length must be even. Length = " + values.size } polynomialTerms = mutableListOf() var i = 0 while (i < values.size) { val t = Term(values[i].toLong(), values[i + 1].toLong()) polynomialTerms.add(t) i += 2 } polynomialTerms.sortWith(TermSorter()) }
constructor() { // zero polynomialTerms = ArrayList() polynomialTerms.add(Term(0, 0)) }
fun hasCoefficientAbs(coeff: Int): Boolean { for (term in polynomialTerms) { if (abs(term.coefficient) == coeff.toLong()) { return true } } return false }
constructor(termList: MutableList<Term>) { if (termList.isEmpty()) { // zero termList.add(Term(0, 0)) } else { // Remove zero terms if needed termList.removeIf { t -> t.coefficient == 0L } } if (termList.size == 0) { // zero termList.add(Term(0, 0)) } polynomialTerms = termList polynomialTerms.sortWith(TermSorter()) }
fun divide(v: Polynomial?): Polynomial { var q = Polynomial() var r = this val lcv = v!!.leadingCoefficient() val dv = v.degree() while (r.degree() >= v.degree()) { val lcr = r.leadingCoefficient() val s = lcr / lcv // Integer division val term = Term(s, r.degree() - dv) q = q.add(term) r = r.add(v.multiply(term.negate())) } return q }
fun add(polynomial: Polynomial): Polynomial { val termList: MutableList<Term> = ArrayList() var thisCount = polynomialTerms.size var polyCount = polynomial.polynomialTerms.size while (thisCount > 0 || polyCount > 0) { val thisTerm = if (thisCount == 0) null else polynomialTerms[thisCount - 1] val polyTerm = if (polyCount == 0) null else polynomial.polynomialTerms[polyCount - 1] when { thisTerm == null -> { termList.add(polyTerm!!.clone()) polyCount-- } polyTerm == null -> { termList.add(thisTerm.clone()) thisCount-- } thisTerm.degree() == polyTerm.degree() -> { val t = thisTerm.add(polyTerm) if (t.coefficient != 0L) { termList.add(t) } thisCount-- polyCount-- } thisTerm.degree() < polyTerm.degree() -> { termList.add(thisTerm.clone()) thisCount-- } else -> { termList.add(polyTerm.clone()) polyCount-- } } } return Polynomial(termList) }
fun add(term: Term): Polynomial { val termList: MutableList<Term> = ArrayList() var added = false for (currentTerm in polynomialTerms) { if (currentTerm.exponent == term.exponent) { added = true if (currentTerm.coefficient + term.coefficient != 0L) { termList.add(currentTerm.add(term)) } } else { termList.add(currentTerm.clone()) } } if (!added) { termList.add(term.clone()) } return Polynomial(termList) }
fun multiply(term: Term): Polynomial { val termList: MutableList<Term> = ArrayList() for (currentTerm in polynomialTerms) { termList.add(currentTerm.clone().multiply(term)) } return Polynomial(termList) }
fun leadingCoefficient(): Long { return polynomialTerms[0].coefficient }
fun degree(): Long { return polynomialTerms[0].exponent }
override fun toString(): String { val sb = StringBuilder() var first = true for (term in polynomialTerms) { if (first) { sb.append(term) first = false } else { sb.append(" ") if (term.coefficient > 0) { sb.append("+ ") sb.append(term) } else { sb.append("- ") sb.append(term.negate()) } } } return sb.toString() }
}
private class TermSorter : Comparator<Term> {
override fun compare(o1: Term, o2: Term): Int { return (o2.exponent - o1.exponent).toInt() }
}
// Note: Cyclotomic Polynomials have small coefficients. Not appropriate for general polynomial usage. private class Term(var coefficient: Long, var exponent: Long) {
fun clone(): Term { return Term(coefficient, exponent) }
fun multiply(term: Term): Term { return Term(coefficient * term.coefficient, exponent + term.exponent) }
fun add(term: Term): Term { if (exponent != term.exponent) { throw RuntimeException("ERROR 102: Exponents not equal.") } return Term(coefficient + term.coefficient, exponent) }
fun negate(): Term { return Term(-coefficient, exponent) }
fun degree(): Long { return exponent }
override fun toString(): String { if (coefficient == 0L) { return "0" } if (exponent == 0L) { return "" + coefficient } if (coefficient == 1L) { return if (exponent == 1L) { "x" } else { "x^$exponent" } } return if (exponent == 1L) { coefficient.toString() + "x" } else coefficient.toString() + "x^" + exponent }
}</lang>
- Output:
Task 1: cyclotomic polynomials for n <= 30: CP[1] = x - 1 CP[2] = x + 1 CP[3] = x^2 + x + 1 CP[4] = x^2 + 1 CP[5] = x^4 + x^3 + x^2 + x + 1 CP[6] = x^2 - x + 1 CP[7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[8] = x^4 + 1 CP[9] = x^6 + x^3 + 1 CP[10] = x^4 - x^3 + x^2 - x + 1 CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[12] = x^4 - x^2 + 1 CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 CP[16] = x^8 + 1 CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[18] = x^6 - x^3 + 1 CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[20] = x^8 - x^6 + x^4 - x^2 + 1 CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[24] = x^8 - x^4 + 1 CP[25] = x^20 + x^15 + x^10 + x^5 + 1 CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[27] = x^18 + x^9 + 1 CP[28] = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 CP[29] = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[30] = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient: CP[1] has coefficient with magnitude = 1 CP[105] has coefficient with magnitude = 2 CP[385] has coefficient with magnitude = 3 CP[1365] has coefficient with magnitude = 4 CP[1785] has coefficient with magnitude = 5 CP[2805] has coefficient with magnitude = 6 CP[3135] has coefficient with magnitude = 7 CP[6545] has coefficient with magnitude = 8 CP[6545] has coefficient with magnitude = 9 CP[10465] has coefficient with magnitude = 10
Perl
Conveniently, the module Math::Polynomial::Cyclotomic
exists to do all the work. An exponent too large
error prevents reaching the 10th step of the 2nd part of the task.
<lang perl>use feature 'say';
use List::Util qw(first);
use Math::Polynomial::Cyclotomic qw(cyclo_poly_iterate);
say 'First 30 cyclotomic polynomials:'; my $it = cyclo_poly_iterate(1); say "$_: " . $it->() for 1 .. 30;
say "\nSmallest cyclotomic polynomial with n or -n as a coefficient:"; $it = cyclo_poly_iterate(1);
for (my ($n, $k) = (1, 1) ; $n <= 10 ; ++$k) {
my $poly = $it->(); while (my $c = first { abs($_) == $n } $poly->coeff) { say "CP $k has coefficient with magnitude = $n"; $n++; }
}</lang>
- Output:
First 30 cyclotomic polynomials: 1: (x - 1) 2: (x + 1) 3: (x^2 + x + 1) 4: (x^2 + 1) 5: (x^4 + x^3 + x^2 + x + 1) 6: (x^2 - x + 1) 7: (x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 8: (x^4 + 1) 9: (x^6 + x^3 + 1) 10: (x^4 - x^3 + x^2 - x + 1) 11: (x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 12: (x^4 - x^2 + 1) 13: (x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 14: (x^6 - x^5 + x^4 - x^3 + x^2 - x + 1) 15: (x^8 - x^7 + x^5 - x^4 + x^3 - x + 1) 16: (x^8 + 1) 17: (x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 18: (x^6 - x^3 + 1) 19: (x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 20: (x^8 - x^6 + x^4 - x^2 + 1) 21: (x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1) 22: (x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1) 23: (x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 24: (x^8 - x^4 + 1) 25: (x^20 + x^15 + x^10 + x^5 + 1) 26: (x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1) 27: (x^18 + x^9 + 1) 28: (x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1) 29: (x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 30: (x^8 + x^7 - x^5 - x^4 - x^3 + x + 1) Smallest cyclotomic polynomial with n or -n as a coefficient: CP 1 has coefficient with magnitude = 1 CP 105 has coefficient with magnitude = 2 CP 385 has coefficient with magnitude = 3 CP 1365 has coefficient with magnitude = 4 CP 1785 has coefficient with magnitude = 5 CP 2805 has coefficient with magnitude = 6 CP 3135 has coefficient with magnitude = 7 CP 6545 has coefficient with magnitude = 8 CP 6545 has coefficient with magnitude = 9
Phix
Uses several routines from Polynomial_long_division#Phix, tweaked slightly to check remainder is zero and trim the quotient. <lang Phix>-- demo\rosetta\Cyclotomic_Polynomial.exw function degree(sequence p)
for i=length(p) to 1 by -1 do if p[i]!=0 then return i end if end for return -1
end function
function poly_div(sequence n, d)
while length(d)<length(n) do d &=0 end while integer dn = degree(n), dd = degree(d) if dd<0 then throw("divide by zero") end if sequence quot = repeat(0,dn) while dn>=dd do integer k = dn-dd integer qk = n[dn]/d[dd] quot[k+1] = qk sequence d2 = d[1..length(d)-k] for i=1 to length(d2) do n[-i] -= d2[-i]*qk end for dn = degree(n) end while
-- return {quot,n} -- (n is now the remainder)
if n!=repeat(0,length(n)) then ?9/0 end if while quot[$]=0 do quot = quot[1..$-1] end while return quot
end function
function poly(sequence si) -- display helper
string r = "" for t=length(si) to 1 by -1 do integer sit = si[t] if sit!=0 then if sit=1 and t>1 then r &= iff(r=""? "":" + ") elsif sit=-1 and t>1 then r &= iff(r=""?"-":" - ") else if r!="" then r &= iff(sit<0?" - ":" + ") sit = abs(sit) end if r &= sprintf("%d",sit) end if r &= iff(t>1?"x"&iff(t>2?sprintf("^%d",t-1):""):"") end if end for if r="" then r="0" end if return r
end function --</Polynomial_long_division.exw>
--# memoize cache for recursive calls constant cyclotomics = new_dict({{1,{-1,1}},{2,{1,1}}})
function cyclotomic(integer n) -- -- Calculate nth cyclotomic polynomial. -- See wikipedia article at bottom of section /wiki/Cyclotomic_polynomial#Fundamental_tools -- The algorithm is reliable but slow for large n > 1000. --
sequence c if getd_index(n,cyclotomics)!=NULL then c = getd(n,cyclotomics) else if is_prime(n) then c = repeat(1,n) else -- recursive formula seen in wikipedia article c = -1&repeat(0,n-1)&1 sequence f = factors(n,-1) for i=1 to length(f) do c = poly_div(c,cyclotomic(f[i])) end for end if setd(n,c,cyclotomics) end if return c
end function
for i=1 to 30 do
sequence z = cyclotomic(i) string s = poly(z) printf(1,"cp(%2d) = %s\n",{i,s}) if i>1 and z!=reverse(z) then ?9/0 end if -- sanity check
end for
integer found = 0, n = 1, cheat = 0 sequence fn = repeat(false,10),
nxt = {105,385,1365,1785,2805,3135,6545,6545,10465,10465}
atom t1 = time()+1 puts(1,"\n") while found<10 do
sequence z = cyclotomic(n) for i=1 to length(z) do atom azi = abs(z[i]) if azi>=1 and azi<=10 and fn[azi]=0 then printf(1,"cp(%d) has a coefficient with magnitude %d\n",{n,azi}) cheat = azi -- (comment this out to prevent cheating!) found += 1 fn[azi] = true t1 = time()+1 end if end for if cheat then {n,cheat} = {nxt[cheat],0} else n += iff(n=1?4:10) end if if time()>t1 then printf(1,"working (%d) ...\r",n) t1 = time()+1 end if
end while</lang>
- Output:
If you disable the cheating, and if in a particularly self harming mood replace it with n+=1, you will get exactly the same output, eventually.
(The distributed version contains simple instrumentation showing cp(1260) executes the line in the heart of poly_div() that subtracts a multiple of qk over 15 million times.)
cp( 1) = x - 1 cp( 2) = x + 1 cp( 3) = x^2 + x + 1 cp( 4) = x^2 + 1 cp( 5) = x^4 + x^3 + x^2 + x + 1 cp( 6) = x^2 - x + 1 cp( 7) = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp( 8) = x^4 + 1 cp( 9) = x^6 + x^3 + 1 cp(10) = x^4 - x^3 + x^2 - x + 1 cp(11) = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(12) = x^4 - x^2 + 1 cp(13) = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(14) = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 cp(15) = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 cp(16) = x^8 + 1 cp(17) = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(18) = x^6 - x^3 + 1 cp(19) = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(20) = x^8 - x^6 + x^4 - x^2 + 1 cp(21) = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 cp(22) = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 cp(23) = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(24) = x^8 - x^4 + 1 cp(25) = x^20 + x^15 + x^10 + x^5 + 1 cp(26) = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 cp(27) = x^18 + x^9 + 1 cp(28) = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 cp(29) = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(30) = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 cp(1) has a coefficient with magnitude 1 cp(105) has a coefficient with magnitude 2 cp(385) has a coefficient with magnitude 3 cp(1365) has a coefficient with magnitude 4 cp(1785) has a coefficient with magnitude 5 cp(2805) has a coefficient with magnitude 6 cp(3135) has a coefficient with magnitude 7 cp(6545) has a coefficient with magnitude 8 cp(6545) has a coefficient with magnitude 9 cp(10465) has a coefficient with magnitude 10
Python
<lang python>from itertools import count, chain from collections import deque
def primes(_cache=[2, 3]):
yield from _cache for n in count(_cache[-1]+2, 2): if isprime(n): _cache.append(n) yield n
def isprime(n):
for p in primes(): if n%p == 0: return False if p*p > n: return True
def factors(n):
for p in primes(): # prime factoring is such a non-issue for small numbers that, for # this example, we might even just say # for p in count(2): if p*p > n: if n > 1: yield(n, 1, 1) break
if n%p == 0: cnt = 0 while True: n, cnt = n//p, cnt+1 if n%p != 0: break yield p, cnt, n
- ^^ not the most sophisticated prime number routines, because no need
- Returns (list1, list2) representing the division between
- two polinomials. A list p of integers means the product
- (x^p[0] - 1) * (x^p[1] - 1) * ...
def cyclotomic(n):
def poly_div(num, den): return (num[0] + den[1], num[1] + den[0])
def elevate(poly, n): # replace poly p(x) with p(x**n) powerup = lambda p, n: [a*n for a in p] return poly if n == 1 else (powerup(poly[0], n), powerup(poly[1], n))
if n == 0: return ([], []) if n == 1: return ([1], [])
p, m, r = next(factors(n)) poly = cyclotomic(r) return elevate(poly_div(elevate(poly, p), poly), p**(m-1))
def to_text(poly):
def getx(c, e): if e == 0: return '1' elif e == 1: return 'x' return 'x' + (.join('⁰¹²³⁴⁵⁶⁷⁸⁹'[i] for i in map(int, str(e))))
parts = [] for (c,e) in (poly): if c < 0: coef = ' - ' if c == -1 else f' - {-c} ' else: coef = (parts and ' + ' or ) if c == 1 else f' + {c}' parts.append(coef + getx(c,e)) return .join(parts)
def terms(poly):
# convert above representation of division to (coef, power) pairs
def merge(a, b): # a, b should be deques. They may change during the course. while a or b: l = a[0] if a else (0, -1) # sentinel value r = b[0] if b else (0, -1) if l[1] > r[1]: a.popleft() elif l[1] < r[1]: b.popleft() l = r else: a.popleft() b.popleft() l = (l[0] + r[0], l[1]) yield l
def mul(poly, p): # p means polynomial x^p - 1 poly = list(poly) return merge(deque((c, e+p) for c,e in poly), deque((-c, e) for c,e in poly))
def div(poly, p): # p means polynomial x^p - 1 q = deque() for c,e in merge(deque(poly), q): if c: q.append((c, e - p)) yield (c, e - p) if e == p: break
p = [(1, 0)] # 1*x^0, i.e. 1
for x in poly[0]: # numerator p = mul(p, x) for x in sorted(poly[1], reverse=True): # denominator p = div(p, x) return p
for n in chain(range(11), [2]):
print(f'{n}: {to_text(terms(cyclotomic(n)))}')
want = 1 for n in count():
c = [c for c,_ in terms(cyclotomic(n))] while want in c or -want in c: print(f'C[{want}]: {n}') want += 1</lang>
- Output:
Only showing first 10 polynomials to avoid clutter.
0: 1 1: x - 1 2: x + 1 3: x² + x + 1 4: x² + 1 5: x⁴ + x³ + x² + x + 1 6: x² - x + 1 7: x⁶ + x⁵ + x⁴ + x³ + x² + x + 1 8: x⁴ + 1 9: x⁶ + x³ + 1 10: x⁴ - x³ + x² - x + 1 105: x⁴⁸ + x⁴⁷ + x⁴⁶ - x⁴³ - x⁴² - 2 x⁴¹ - x⁴⁰ - x³⁹ + x³⁶ + x³⁵ + x³⁴ + x³³ + x³² + x³¹ - x²⁸ - x²⁶ - x²⁴ - x²² - x²⁰ + x¹⁷ + x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² - x⁹ - x⁸ - 2 x⁷ - x⁶ - x⁵ + x² + x + 1 C[1]: 0 C[2]: 105 C[3]: 385 C[4]: 1365 C[5]: 1785 C[6]: 2805 C[7]: 3135 C[8]: 6545 C[9]: 6545 C[10]: 10465 C[11]: 10465 C[12]: 10465 C[13]: 10465 C[14]: 10465 C[15]: 11305 C[16]: 11305 C[17]: 11305 C[18]: 11305 C[19]: 11305 C[20]: 11305 C[21]: 11305 C[22]: 15015 C[23]: 15015
Raku
(formerly Perl 6)
Uses the same library as Perl, so comes with the same caveats. <lang perl6>use Math::Polynomial::Cyclotomic:from<Perl5> <cyclo_poly_iterate cyclo_poly>;
say 'First 30 cyclotomic polynomials:'; my $iterator = cyclo_poly_iterate(1); say "Φ($_) = " ~ super $iterator().Str for 1..30;
say "\nSmallest cyclotomic polynomial with |n| as a coefficient:"; say "Φ(1) has a coefficient magnitude: 1";
my $index = 0; for 2..9 -> $coefficient {
loop { $index += 5; my \Φ = cyclo_poly($index); next unless Φ ~~ / $coefficient\* /; say "Φ($index) has a coefficient magnitude: $coefficient"; $index -= 5; last; }
}
sub super ($str) {
$str.subst( / '^' (\d+) /, { $0.trans([<0123456789>.comb] => [<⁰¹²³⁴⁵⁶⁷⁸⁹>.comb]) }, :g)
}</lang>
First 30 cyclotomic polynomials: Φ(1) = (x - 1) Φ(2) = (x + 1) Φ(3) = (x² + x + 1) Φ(4) = (x² + 1) Φ(5) = (x⁴ + x³ + x² + x + 1) Φ(6) = (x² - x + 1) Φ(7) = (x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(8) = (x⁴ + 1) Φ(9) = (x⁶ + x³ + 1) Φ(10) = (x⁴ - x³ + x² - x + 1) Φ(11) = (x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(12) = (x⁴ - x² + 1) Φ(13) = (x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(14) = (x⁶ - x⁵ + x⁴ - x³ + x² - x + 1) Φ(15) = (x⁸ - x⁷ + x⁵ - x⁴ + x³ - x + 1) Φ(16) = (x⁸ + 1) Φ(17) = (x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(18) = (x⁶ - x³ + 1) Φ(19) = (x¹⁸ + x¹⁷ + x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(20) = (x⁸ - x⁶ + x⁴ - x² + 1) Φ(21) = (x¹² - x¹¹ + x⁹ - x⁸ + x⁶ - x⁴ + x³ - x + 1) Φ(22) = (x¹⁰ - x⁹ + x⁸ - x⁷ + x⁶ - x⁵ + x⁴ - x³ + x² - x + 1) Φ(23) = (x²² + x²¹ + x²⁰ + x¹⁹ + x¹⁸ + x¹⁷ + x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(24) = (x⁸ - x⁴ + 1) Φ(25) = (x²⁰ + x¹⁵ + x¹⁰ + x⁵ + 1) Φ(26) = (x¹² - x¹¹ + x¹⁰ - x⁹ + x⁸ - x⁷ + x⁶ - x⁵ + x⁴ - x³ + x² - x + 1) Φ(27) = (x¹⁸ + x⁹ + 1) Φ(28) = (x¹² - x¹⁰ + x⁸ - x⁶ + x⁴ - x² + 1) Φ(29) = (x²⁸ + x²⁷ + x²⁶ + x²⁵ + x²⁴ + x²³ + x²² + x²¹ + x²⁰ + x¹⁹ + x¹⁸ + x¹⁷ + x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(30) = (x⁸ + x⁷ - x⁵ - x⁴ - x³ + x + 1) Smallest cyclotomic polynomial with |n| as a coefficient: Φ(1) has a coefficient magnitude: 1 Φ(105) has a coefficient magnitude: 2 Φ(385) has a coefficient magnitude: 3 Φ(1365) has a coefficient magnitude: 4 Φ(1785) has a coefficient magnitude: 5 Φ(2805) has a coefficient magnitude: 6 Φ(3135) has a coefficient magnitude: 7 Φ(6545) has a coefficient magnitude: 8 Φ(6545) has a coefficient magnitude: 9
Sidef
Solution based on polynomial interpolation (slow). <lang ruby>var Poly = require('Math::Polynomial') Poly.string_config(Hash(fold_sign => true, prefix => "", suffix => ""))
func poly_interpolation(v) {
v.len.of {|n| v.len.of {|k| n**k } }.msolve(v)
}
say "First 30 cyclotomic polynomials:" for k in (1..30) {
var a = (k+1).of { cyclotomic(k, _) } var Φ = poly_interpolation(a) say ("Φ(#{k}) = ", Poly.new(Φ...))
}
say "\nSmallest cyclotomic polynomial with n or -n as a coefficient:" for n in (1..10) { # very slow
var k = (1..Inf -> first {|k| poly_interpolation((k+1).of { cyclotomic(k, _) }).first { .abs == n } }) say "Φ(#{k}) has coefficient with magnitude #{n}"
}</lang>
Slightly faster solution, using the Math::Polynomial::Cyclotomic Perl module. <lang ruby>var Poly = require('Math::Polynomial')
require('Math::Polynomial::Cyclotomic')
Poly.string_config(Hash(fold_sign => true, prefix => "", suffix => ""))
say "First 30 cyclotomic polynomials:" for k in (1..30) {
say ("Φ(#{k}) = ", Poly.new.cyclotomic(k))
}
say "\nSmallest cyclotomic polynomial with n or -n as a coefficient:" for n in (1..10) {
var p = Poly.new var k = (1..Inf -> first {|k| [p.cyclotomic(k).coeff].first { .abs == n } }) say "Φ(#{k}) has coefficient with magnitude = #{n}"
}</lang>
- Output:
First 30 cyclotomic polynomials: Φ(1) = x - 1 Φ(2) = x + 1 Φ(3) = x^2 + x + 1 Φ(4) = x^2 + 1 Φ(5) = x^4 + x^3 + x^2 + x + 1 Φ(6) = x^2 - x + 1 Φ(7) = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(8) = x^4 + 1 Φ(9) = x^6 + x^3 + 1 Φ(10) = x^4 - x^3 + x^2 - x + 1 Φ(11) = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(12) = x^4 - x^2 + 1 Φ(13) = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(14) = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 Φ(15) = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 Φ(16) = x^8 + 1 Φ(17) = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(18) = x^6 - x^3 + 1 Φ(19) = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(20) = x^8 - x^6 + x^4 - x^2 + 1 Φ(21) = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 Φ(22) = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 Φ(23) = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(24) = x^8 - x^4 + 1 Φ(25) = x^20 + x^15 + x^10 + x^5 + 1 Φ(26) = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 Φ(27) = x^18 + x^9 + 1 Φ(28) = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 Φ(29) = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(30) = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 Smallest cyclotomic polynomial with n or -n as a coefficient: Φ(1) has coefficient with magnitude = 1 Φ(105) has coefficient with magnitude = 2 Φ(385) has coefficient with magnitude = 3 Φ(1365) has coefficient with magnitude = 4 Φ(1785) has coefficient with magnitude = 5 Φ(2805) has coefficient with magnitude = 6 Φ(3135) has coefficient with magnitude = 7 ^C