Cyclotomic polynomial

From Rosetta Code
Task
Cyclotomic polynomial
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++[edit]

Translation of: Java
#include <algorithm>
#include <iostream>
#include <initializer_list>
#include <map>
#include <vector>

const int MAX_ALL_FACTORS = 100000;
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;
}
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#[edit]

Translation of: Java
Works with: C sharp version 8
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}"
        };
    }
}
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

D[edit]

Translation of: Kotlin
import std.algorithm;
import std.exception;
import std.format;
import std.functional;
import std.math;
import std.range;
import std.stdio;

immutable MAX_ALL_FACTORS = 100_000;
immutable ALGORITHM = 2;

//Note: Cyclotomic Polynomials have small coefficients.  Not appropriate for general polynomial usage.

struct Term {
    private long m_coefficient;
    private long m_exponent;

    public this(long c, long e) {
        m_coefficient = c;
        m_exponent = e;
    }

    public long coefficient() const {
        return m_coefficient;
    }

    public long exponent() const {
        return m_exponent;
    }

    public Term opUnary(string op)() const {
        static if (op == "-") {
            return Term(-m_coefficient, m_exponent);
        } else {
            assert(false, "Not implemented");
        }
    }

    public Term opBinary(string op)(Term term) const {
        static if (op == "+") {
            if (exponent() != term.exponent()) {
                assert(false, "Error 102: Exponents not equals.");
            }
            return Term(coefficient() + term.coefficient(), exponent());
        } else if (op == "*") {
            return Term(coefficient() * term.coefficient(), exponent() + term.exponent());
        } else {
            assert(false, "Not implemented: " ~ op);
        }
    }

    public void toString(scope void delegate(const(char)[]) sink) const {
        auto spec = singleSpec("%s");
        if (m_coefficient == 0) {
            sink("0");
        } else if (m_exponent == 0) {
            formatValue(sink, m_coefficient, spec);
        } else if (m_coefficient == 1) {
            if (m_exponent == 1) {
                sink("x");
            } else {
                sink("x^");
                formatValue(sink, m_exponent, spec);
            }
        } else if (m_coefficient == -1) {
            if (m_exponent == 1) {
                sink("-x");
            } else {
                sink("-x^");
                formatValue(sink, m_exponent, spec);
            }
        } else if (m_exponent == 1) {
            formatValue(sink, m_coefficient, spec);
            sink("x");
        } else {
            formatValue(sink, m_coefficient, spec);
            sink("x^");
            formatValue(sink, m_exponent, spec);
        }
    }
}

struct Polynomial {
    private Term[] terms;

    public this(const Term[] ts...) {
        terms = ts.dup;
        terms.sort!"b.exponent < a.exponent";
    }

    bool hasCoefficientAbs(int coeff) const {
        foreach (term; terms) {
            if (abs(term.coefficient) == coeff) {
                return true;
            }
        }
        return false;
    }

    public long leadingCoefficient() const {
        return terms[0].coefficient();
    }

    public long degree() const {
        if (terms.empty) {
            return -1;
        }
        return terms[0].exponent();
    }

    public Polynomial opBinary(string op)(Term term) const {
        static if (op == "+") {
            Term[] newTerms;
            auto added = false;
            foreach (currentTerm; terms) {
                if (currentTerm.exponent() == term.exponent()) {
                    added = true;
                    if (currentTerm.coefficient() + term.coefficient() != 0) {
                        newTerms ~= currentTerm + term;
                    }
                } else {
                    newTerms ~= currentTerm;
                }
            }
            if (!added) {
                newTerms ~= term;
            }
            return Polynomial(newTerms);
        } else if (op == "*") {
            Term[] newTerms;
            foreach (currentTerm; terms) {
                newTerms ~= currentTerm * term;
            }
            return Polynomial(newTerms);
        } else {
            assert(false, "Not implemented: " ~ op);
        }
    }

    public Polynomial opBinary(string op)(Polynomial rhs) const {
        static if (op == "+") {
            Term[] newTerms;
            auto thisCount = terms.length;
            auto polyCount = rhs.terms.length;
            while (thisCount > 0 || polyCount > 0) {
                if (thisCount == 0) {
                    newTerms ~= rhs.terms[polyCount - 1];
                    polyCount--;
                } else if (polyCount == 0) {
                    newTerms ~= terms[thisCount - 1];
                    thisCount--;
                } else {
                    auto thisTerm = terms[thisCount - 1];
                    auto polyTerm = rhs.terms[polyCount - 1];
                    if (thisTerm.exponent() == polyTerm.exponent()) {
                        auto t = thisTerm + polyTerm;
                        if (t.coefficient() != 0) {
                            newTerms ~= t;
                        }
                        thisCount--;
                        polyCount--;
                    } else if (thisTerm.exponent() < polyTerm.exponent()) {
                        newTerms ~= thisTerm;
                        thisCount--;
                    } else {
                        newTerms ~= polyTerm;
                        polyCount--;
                    }
                }
            }
            return Polynomial(newTerms);
        } else if (op == "/") {
            Polynomial q;
            auto r = Polynomial(terms);
            auto lcv = rhs.leadingCoefficient();
            auto dv = rhs.degree();
            while (r.degree() >= rhs.degree()) {
                auto lcr = r.leadingCoefficient();
                auto s = lcr / lcv;
                auto term = Term(s, r.degree() - dv);
                q = q + term;
                r = r + rhs * -term;
            }
            return q;
        } else {
            assert(false, "Not implemented: " ~ op);
        }
    }

    public int opApply(int delegate(Term) dg) const {
        foreach (term; terms) {
            auto rv = dg(term);
            if (rv != 0) {
                return rv;
            }
        }
        return 0;
    }

    public void toString(scope void delegate(const(char)[]) sink) const {
        auto spec = singleSpec("%s");
        if (!terms.empty) {
            formatValue(sink, terms[0], spec);
            foreach (t; terms[1..$]) {
                if (t.coefficient > 0) {
                    sink(" + ");
                    formatValue(sink, t, spec);
                } else {
                    sink(" - ");
                    formatValue(sink, -t, spec);
                }
            }
        }
    }
}

void putAll(K, V)(ref V[K] a, V[K] b) {
    foreach (k, v; b) {
        a[k] = v;
    }
}

void merge(K, V, F)(ref V[K] a, K k, V v, F f) {
    if (k in a) {
        a[k] = f(a[k], v);
    } else {
        a[k] = v;
    }
}

int sum(int a, int b) {
    return a + b;
}

int[int] getFactorsImpl(int number) {
    int[int] factors;
    if (number % 2 == 0) {
        if (number > 2) {
            auto factorsDivTwo = memoize!getFactorsImpl(number / 2);
            factors.putAll(factorsDivTwo);
        }
        factors.merge(2, 1, &sum);
        return factors;
    }
    auto root = sqrt(cast(real) number);
    auto i = 3;
    while (i <= root) {
        if (number % i == 0) {
            factors.putAll(memoize!getFactorsImpl(number / i));
            factors.merge(i, 1, &sum);
            return factors;
        }
        i += 2;
    }
    factors[number] = 1;
    return factors;
}
alias getFactors = memoize!getFactorsImpl;

int[] getDivisors(int number) {
    int[] divisors;
    auto root = cast(int)sqrt(cast(real) number);
    foreach (i; 1..root) {
        if (number % i == 0) {
            divisors ~= i;
        }
        auto div = number / i;
        if (div != i && div != number) {
            divisors ~= div;
        }
    }
    return divisors;
}

Polynomial cyclotomicPolynomialImpl(int n) {
    if (n == 1) {
        //  Polynomial:  x - 1
        return Polynomial(Term(1, 1), Term(-1, 0));
    }
    auto factors = getFactors(n);
    if (n in factors) {
        // n prime
        Term[] terms;
        foreach (i; 0..n) {
            terms ~= Term(1, i);
        }
        return Polynomial(terms);
    } else if (factors.length == 2 && 2 in factors && factors[2] == 1 && (n / 2) in factors && factors[n / 2] == 1) {
        //  n = 2p
        auto prime = n / 2;
        Term[] terms;
        auto coeff = -1;
        foreach (i; 0..prime) {
            coeff *= -1;
            terms ~= Term(coeff, i);
        }
        return Polynomial(terms);
    } else if (factors.length == 1 && 2 in factors) {
        //  n = 2^h
        auto h = factors[2];
        Term[] terms;
        terms ~= Term(1, 2 ^^ (h - 1));
        terms ~= Term(1, 0);
        return Polynomial(terms);
    } else if (factors.length == 1 && n !in factors) {
        // n = p^k
        auto p = 0;
        auto k = 0;
        foreach (prime, v; factors) {
            if (prime > p) {
                p = prime;
                k = v;
            }
        }
        Term[] terms;
        foreach (i; 0..p) {
            terms ~= Term(1, (i * p) ^^ (k - 1));
        }
        return Polynomial(terms);
    } else if (factors.length == 2 && 2 in factors) {
        // n = 2^h * p^k
        auto p = 0;
        auto k = 0;
        foreach (prime, v; factors) {
            if (prime != 2 && prime > p) {
                p = prime;
                k = v;
            }
        }
        Term[] terms;
        auto coeff = -1;
        auto twoExp = 2 ^^ (factors[2] - 1);
        foreach (i; 0..p) {
            coeff *= -1;
            auto exponent = i * twoExp * p ^^ (k - 1);
            terms ~= Term(coeff, exponent);
        }
        return Polynomial(terms);
    } else if (2 in factors && n / 2 % 2 == 1 && n / 2 > 1) {
        //  CP(2m)[x] = CP(-m)[x], n odd integer > 1
        auto cycloDiv2 = memoize!cyclotomicPolynomialImpl(n / 2);
        Term[] terms;
        foreach (term; cycloDiv2) {
            if (term.exponent() % 2 == 0) {
                terms ~= term;
            } else {
                terms ~= -term;
            }
        }
        return Polynomial(terms);
    }

    if (ALGORITHM == 0) {
        //  Slow - uses basic definition.
        auto divisors = getDivisors(n);
        //  Polynomial:  ( x^n - 1 )
        auto cyclo = Polynomial(Term(1, n), Term(-1, 0));
        foreach (i; divisors) {
            auto p = memoize!cyclotomicPolynomialImpl(i);
            cyclo = cyclo / p;
        }
        return cyclo;
    }
    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);
        auto maxDivisor = int.min;
        foreach (div; divisors) {
            maxDivisor = max(maxDivisor, div);
        }
        int[] divisorsExceptMax;
        foreach (div; divisors) {
            if (maxDivisor % div != 0) {
                divisorsExceptMax ~= div;
            }
        }

        //  Polynomial:  ( x^n - 1 ) / ( x^m - 1 ), where m is the max divisor
        auto cyclo = Polynomial(Term(1, n), Term(-1, 0)) / Polynomial(Term(1, maxDivisor), Term(-1, 0));
        foreach (i; divisorsExceptMax) {
            auto p = memoize!cyclotomicPolynomialImpl(i);
            cyclo = cyclo / p;
        }
        return cyclo;
    }
    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]
        auto m = 1;
        auto cyclo = memoize!cyclotomicPolynomialImpl(m);
        auto primes = factors.keys;
        primes.sort;
        foreach (prime; primes) {
            //  CP(m)[x]
            auto cycloM = cyclo;
            //  Compute CP(m)[x^p].
            Term[] terms;
            foreach (term; cycloM) {
                terms ~= Term(term.coefficient(), term.exponent() * prime);
            }
            cyclo = Polynomial(terms) / cycloM;
            m *= prime;
        }
        //  Now, m is the largest square free divisor of n
        auto s = n / m;
        //  Compute CP(n)[x] = CP(m)[x^s]
        Term[] terms;
        foreach (term; cyclo) {
            terms ~= Term(term.coefficient(), term.exponent() * s);
        }
        return Polynomial(terms);
    }
    assert(false, "Error 103: Invalid algorithm");
}
alias cyclotomicPolynomial = memoize!cyclotomicPolynomialImpl;

void main() {
    writeln("Task 1:  cyclotomic polynomials for n <= 30:");
   foreach (i; 1 .. 31) {
        auto p = cyclotomicPolynomial(i);
        writefln("CP[%d] = %s", i, p);
   }
    writeln;

    writeln("Task 2:  Smallest cyclotomic polynomial with n or -n as a coefficient:");
    auto n = 0;
    foreach (i; 1 .. 11) {
         while (true) {
            n++;
            auto cyclo = cyclotomicPolynomial(n);
            if (cyclo.hasCoefficientAbs(i)) {
                writefln("CP[%d] has coefficient with magnitude = %d", n, i);
                n--;
                break;
            }
         }
    }
}
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^36 + 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

Fermat[edit]

This isn't terribly efficient if you have to calculate many cyclotomics- store them in an array rather than using recursion instead if you need to do that- but it showcases Fermat's strength at polynomial expressions.

&(J=x);                                                    {adjoin x as the variable in the polynomials}

Func Cyclotomic(n) =
    if n=1 then x-1 fi;                                    {first cyclotomic polynomial is x^n-1}                           
    r:=x^n-1;                                              {caclulate cyclotomic by division}
    for d = 1 to n-1 do
        if Divides(d,n) then
            r:=r\Cyclotomic(d)
        fi;
    od;
    r.;                                                    {return the polynomial}
   
Func Hascoef(n, k) =
    p:=Cyclotomic(n);
    for d = 0 to Deg(p) do
        if |(Coef(p,d))|=k then Return(1) fi
    od;
    0.;
    
for d = 1 to 30 do
    !!(d,' : ',Cyclotomic(d))
od;

for m = 1 to 10 do
    i:=1;
    while not Hascoef(i, m) do
        i:+
    od;
    !!(m,'   :   ',i);
od;

Go[edit]

Translation of: Java
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
            }
        }
    }
}
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

Haskell[edit]

Uses synthetic polynomial division and simple memoization.

import Data.List
import Data.Numbers.Primes (primeFactors)

negateVar p = zipWith (*) p $ reverse $ take (length p) $ cycle [1,-1]

lift p 1 = p
lift p n = intercalate (replicate (n-1) 0) (pure <$> p)

shortDiv :: [Integer] -> [Integer] -> [Integer]
shortDiv p1 (_:p2) = unfoldr go (length p1 - length p2, p1)
  where
    go (0, _) = Nothing
    go (i, h:t) = Just (h, (i-1, zipWith (+) (map (h *) ker) t))
    ker = negate <$> p2 ++ repeat 0

primePowerFactors = sortOn fst . map (\x-> (head x, length x)) . group . primeFactors
                     
-- simple memoization
cyclotomics :: [[Integer]]
cyclotomics = cyclotomic <$> [0..]

cyclotomic :: Int -> [Integer]
cyclotomic 0 = [0]
cyclotomic 1 = [1, -1]
cyclotomic 2 = [1, 1]
cyclotomic n = case primePowerFactors n of
  -- for n = 2^k
  [(2,h)]       -> 1 : replicate (2 ^ (h-1) - 1) 0 ++ [1]
  -- for prime n
  [(p,1)]       -> replicate n 1
  -- for power of prime n
  [(p,m)]       -> lift (cyclotomics !! p) (p^(m-1))
  -- for n = 2*p and prime p
  [(2,1),(p,1)] -> take (n `div` 2) $ cycle [1,-1]
  -- for n = 2*m and odd m
  (2,1):_       -> negateVar $ cyclotomics !! (n `div` 2)
  -- general case
  (p, m):ps     -> let cm = cyclotomics !! (n `div` (p ^ m))
                   in lift (lift cm p `shortDiv` cm) (p^(m-1))

Simple examples

λ> cyclotomic 7
[1,1,1,1,1,1,1]

λ> cyclotomic 9
[1,0,0,1,0,0,1]

λ> cyclotomic 16
[1,0,0,0,0,0,0,0,1]

The task solution

showPoly [] = "0"
showPoly p = foldl1 (\r -> (r ++) . term) $
             dropWhile null $
             foldMap (\(c, n) -> [show c ++ expt n]) $
             zip (reverse p) [0..]
  where
    expt = \case 0 -> ""
                 1 -> "*x"
                 n -> "*x^" ++ show n

    term = \case [] -> ""
                 '0':'*':t -> ""
                 '-':'1':'*':t -> " - " ++ t
                 '1':'*':t -> " + " ++ t
                 '-':t -> " - " ++ t
                 t -> " + " ++ t     

main = do
  mapM_ (print . showPoly . cyclotomic) [1..30]
  putStrLn $ replicate 40 '-'
  
  mapM_ showLine $ take 4 task2
  where
    showLine (j, i, l) = putStrLn $ concat [ show j
                                            , " appears in CM(", show i
                                            , ") of length ", show l ]

    -- in order to make computations faster we leave only each 5-th polynomial
    task2 = (1,1,2) : tail (search 1 $ zip [0,5..] $ skipBy 5 cyclotomics)
      where
        search i ((k, p):ps) = if i `notElem` (abs <$> p)
                               then search i ps
                               else (i, k, length p) : search (i+1) ((k, p):ps)

skipBy n [] = []
skipBy n lst = let (x:_, b) = splitAt n lst in x:skipBy n b

Result

"-1 + x^1"
"1 + x^1"
"1 + x^1 + x^2"
"1 + x^2"
"1 + x^1 + x^2 + x^3 + x^4"
"1 - x^1 + x^2"
"1 + x^1 + x^2 + x^3 + x^4 + x^5 + x^6"
"1 + x^4"
"1 + x^3 + x^6"
"1 - x^1 + x^2 - x^3 + x^4"
"1 + x^1 + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10"
"1 - x^2 + x^4"
"1 + x^1 + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12"
"1 - x^1 + x^2 - x^3 + x^4 - x^5 + x^6"
"1 - x^1 + x^3 - x^4 + x^5 - x^7 + x^8"
"1 + x^8"
"1 + x^1 + 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"
"1 - x^3 + x^6"
"1 + x^1 + 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"
"1 - x^2 + x^4 - x^6 + x^8"
"1 - x^1 + x^3 - x^4 + x^6 - x^8 + x^9 - x^11 + x^12"
"1 - x^1 + x^2 - x^3 + x^4 - x^5 + x^6 - x^7 + x^8 - x^9 + x^10"
"1 + x^1 + 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"
"1 - x^4 + x^8"
"1 + x^5 + x^10 + x^15 + x^20"
"1 - x^1 + x^2 - x^3 + x^4 - x^5 + x^6 - x^7 + x^8 - x^9 + x^10 - x^11 + x^12"
"1 + x^9 + x^18"
"1 - x^2 + x^4 - x^6 + x^8 - x^10 + x^12"
"1 + x^1 + 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"
"1 + x^1 - x^3 - x^4 - x^5 + x^7 + x^8"
----------------------------------------
1 appears in CM(1) having 2 terms
2 appears in CM(105) having 49 terms
3 appears in CM(385) having 241 terms
4 appears in CM(1365) having 577 terms
5 appears in CM(1785) having 769 terms
6 appears in CM(2805) having 1281 terms
7 appears in CM(3135) having 1441 terms
8 appears in CP(6545) having 3841 terms
9 appears in CP(6545) having 3841 terms
10 appears in CP(10465) having 6337 terms

Computations take a while...

J[edit]

For values up to 70, we can find cyclotomic polynomials by finding a polynomial with roots of unity relatively prime to the order of the polynomial:

cyclo=: {{<.-:1+(++) p. 1;^0j2p1* y%~1+I.1=y+.1+i.y}}

This approach suggests that cyclotomic polynomial zero should be f0(x)= 1

Routine to find the nth cyclotomic polynomial:

{{ if.0>nc<'cache' do.cache=:y end.}} (,1);_1 1

cyclotomic=: {{
   if.y<#cache do.
     if.#c=. y{::cache do.
       c return.
     end.
   end.
   c=. unpad cyclotomic000 y
   if. y>:#cache do. cache=:(100+y){.cache end.
   cache=: (<c) y} cache
   c
}}

cyclotomic000=:  {{ assert.0<y
   'q p'=. __ q: y
   if. 1=#q do.
     ,(y%*/q) {."0 q#1
   elseif.2={.q do.
     ,(y%*/q) {."0 (* 1 _1 $~ #) cyclotomic */}.q
   elseif. 1 e. 1 < p do.
     ,(y%*/q) {."0 cyclotomic */q
   else.
     (_1,(-y){.1) pDiv ;+//.@(*/)each/ cyclotomic each}:*/@>,{1,each q
   end.
}}


NB. discard high order zero coefficients in representation of polynomial
unpad=: {.~ 1+0 i:~0=]

NB. polynomial division, optimized for somewhat sparse polynomials
pDiv=: {{
  q=. $j=. 2 + x -&# y
  'x y'=. x,:y
  while. j=. j-1 do.
    if. 0={.x do. j=. j-<:i=. 0 i.~ 0=x
      q=. q,i#0
      x=. i |.!.0 x
    else.
      q=. q, r=. x %&{. y
      x=. 1 |.!.0 x - y*r
    end.
  end.q
}}

If you take all the divisors of a number. (For example, for 12, the divisors are: 1, 2, 3, 4, 6 and 12) and find the product of their cyclotomic polynomials (for example, for 12, x-1, x+1, x2+x+1, x2+1, x2-x+1, and x4-x2+1) you get xn-1 (for 12, that would of course be x12-1).

Notes:

  • the coefficients of cyclotomic polynomials after 1 form a palindrome (that's the q#1 phrase in the implementation).
  • the cyclotomic polynomial for a prime number has as many terms as that number, and the coefficients are all 1 (with no intervening zeros -- the highest power is one less than that prime).
  • powers of primes add zero coefficients to the polynomial (that's the ,(y%*/q) {."0 ... phrase in the implementation). This means that we can mostly ignore powers of prime numbers -- they're just going to correspond to zeros we add to the base polynomial.
  • an even base cyclotomic polynomial is the same as the corresponding odd base cyclotomic polynomial except with x replaced by negative x. (that's the (* 1 _1 $~ #) phrase in the implementation.
  • To deal with the general case, we use polynomial division, xn-1 divided by the polynomial product the cyclotomic polynomials of the proper divisors of number we're looking for.
  • +//.@(*/) is polynomial product in J.

Task examples:

taskfmt=: {{
  c=. ":each j=.cyclotomic y
  raw=. rplc&'_-' ;:inv}.,'+';"0|.(*|j)#c('(',[,],')'"_)each '*x^',&":L:0 <"0 i.#c
  txt=. raw rplc'(1*x^0)';'1';'(-1*x^0)';'(-1)';'*x^1)';'*x)'
  LF,~'CP[',y,&":']= ',rplc&('(x)';'x';'+ (-1)';'- 1')txt rplc'(1*';'(';'(-1*';'(-'
}}

taskorder=: {{
  r=.$k=.0
  while.y>#r do.k=.k+1
    if.(1+#r) e.|cyclotomic k do.
      r=. r,k
      k=. k-1
    end.
  end.r
}}

   ;taskfmt each 1+i.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

   (,.~#\) taskorder 10
 1     1
 2   105
 3   385
 4  1365
 5  1785
 6  2805
 7  3135
 8  6545
 9  6545
10 10465

Another approach[edit]

As noted in the J programming forum, we can improve the big-O character of this algorithm by using the fast fourier transform for polynomial multiplication and division.

NB. install'math/fftw'
require'math/fftw'

cyclotomic000=:  {{ assert.0<y
  if. y = 1 do. _1 1 return. end.
  'q p'=. __ q: y
  if. 1=#q do.
    ,(y%*/q) {."0 q#1
  elseif.2={.q do.
    ,(y%*/q) {."0 (* 1 _1 $~ #) cyclotomic */}.q
  elseif. 1 e. 1 < p do.
    ,(y%*/q) {."0 cyclotomic */q
  else.
    lgl=. {:$ ctlist=. cyclotomic "0 }:*/@>,{1,each q    NB. ctlist is 2-d table of polynomial divisors
    lgd=. # dividend=. _1,(-y){.1                        NB. (x^n) - 1, and its size
    lg=.  >.&.(2&^.)  lgl >. lgd                         NB. required lengths of all polynomials for fft transforms
                      NB. really, "divisor" is the fft of the divisor!
            divisor=. */ fftw"1 lg{."1 ctlist            NB. FFT article doesn't deal with lists of multiplicands
    unpad roundreal ifftw"1 divisor %~ fftw lg{.dividend NB. similar to article's multiplication
  end.
}}

roundreal =: [: <. 0.5 + 9&o.

This variation for polynomial division is only valid when there's no remainder to be concerned with (which is the case, here). The article mentioned in the comments is an essay on using fft for polynomial multiplication

This approach gave slightly over a 16x speedup for taskorder 10, from a 2 element cache, with an approximately 50% increased memory footprint. (Remember, of course, that benchmarks and benchmark ratios have dependencies on computer architecture and language implementation, and the host environment.)

Java[edit]

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;
        }
    }

}
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[edit]

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
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[edit]

Translation of: Java
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
    }
}
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

Maple[edit]

with(NumberTheory):
for n to 30 do lprint(Phi(n,x)) od:

x-1
x+1
x^2+x+1
x^2+1
x^4+x^3+x^2+x+1
x^2-x+1
x^6+x^5+x^4+x^3+x^2+x+1
x^4+1
x^6+x^3+1
x^4-x^3+x^2-x+1
x^10+x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2+x+1
x^4-x^2+1
x^12+x^11+x^10+x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2+x+1
x^6-x^5+x^4-x^3+x^2-x+1
x^8-x^7+x^5-x^4+x^3-x+1
x^8+1
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
x^6-x^3+1
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
x^8-x^6+x^4-x^2+1
x^12-x^11+x^9-x^8+x^6-x^4+x^3-x+1
x^10-x^9+x^8-x^7+x^6-x^5+x^4-x^3+x^2-x+1
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
x^8-x^4+1
x^20+x^15+x^10+x^5+1
x^12-x^11+x^10-x^9+x^8-x^7+x^6-x^5+x^4-x^3+x^2-x+1
x^18+x^9+1
x^12-x^10+x^8-x^6+x^4-x^2+1
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
x^8+x^7-x^5-x^4-x^3+x+1

PhiSet:=[seq(map(abs,{coeffs(Phi(k,x),x)}),k=1..15000)]:
[seq(ListTools:-SelectFirst(s->member(n,s),PhiSet,output=indices),n=1..20)];
#[1, 105, 385, 1365, 1785, 2805, 3135, 6545, 6545, 10465, 10465, 
#  10465, 10465, 10465, 11305, 11305, 11305, 11305, 11305, 11305]

Mathematica / Wolfram Language[edit]

Cyclotomic[#, x] & /@ Range[30] // Column
i = 1;
n = 10;
PrintTemporary[Dynamic[{magnitudes, i}]];
magnitudes = ConstantArray[True, n];
While[Or @@ magnitudes,
 coeff = Abs[CoefficientList[Cyclotomic[i, x], x]];
 coeff = Select[coeff, Between[{1, n}]];
 coeff = DeleteDuplicates[coeff];
 If[Or @@ magnitudes[[coeff]],
  Do[
   If[magnitudes[[c]] == True,
    Print["CyclotomicPolynomial(", i, 
     ") has coefficient with magnitude ", c]
    ]
   ,
   {c, coeff}
   ];
  magnitudes[[coeff]] = False;
  ];
 i++;
 ]
Output:
-1+x
1+x
1+x+x^2
1+x^2
1+x+x^2+x^3+x^4
1-x+x^2
1+x+x^2+x^3+x^4+x^5+x^6
1+x^4
1+x^3+x^6
1-x+x^2-x^3+x^4
1+x+x^2+x^3+x^4+x^5+x^6+x^7+x^8+x^9+x^10
1-x^2+x^4
1+x+x^2+x^3+x^4+x^5+x^6+x^7+x^8+x^9+x^10+x^11+x^12
1-x+x^2-x^3+x^4-x^5+x^6
1-x+x^3-x^4+x^5-x^7+x^8
1+x^8
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
1-x^3+x^6
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
1-x^2+x^4-x^6+x^8
1-x+x^3-x^4+x^6-x^8+x^9-x^11+x^12
1-x+x^2-x^3+x^4-x^5+x^6-x^7+x^8-x^9+x^10
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
1-x^4+x^8
1+x^5+x^10+x^15+x^20
1-x+x^2-x^3+x^4-x^5+x^6-x^7+x^8-x^9+x^10-x^11+x^12
1+x^9+x^18
1-x^2+x^4-x^6+x^8-x^10+x^12
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
1+x-x^3-x^4-x^5+x^7+x^8

CyclotomicPolynomial(1) has coefficient with magnitude 1
CyclotomicPolynomial(105) has coefficient with magnitude 2
CyclotomicPolynomial(385) has coefficient with magnitude 3
CyclotomicPolynomial(1365) has coefficient with magnitude 4
CyclotomicPolynomial(1785) has coefficient with magnitude 5
CyclotomicPolynomial(2805) has coefficient with magnitude 6
CyclotomicPolynomial(3135) has coefficient with magnitude 7
CyclotomicPolynomial(6545) has coefficient with magnitude 8
CyclotomicPolynomial(6545) has coefficient with magnitude 9
CyclotomicPolynomial(10465) has coefficient with magnitude 10

Nim[edit]

Translation of: Java

We use Java algorithm with ideas from C#, D, Go and Kotlin. We have kept only algorithm number 2 as other algorithms are much less efficient. We have also done some Nim specific improvements in order to get better performances.

import algorithm, math, sequtils, strformat, tables

type

  Term = tuple[coeff: int; exp: Natural]
  Polynomial = seq[Term]

  # Table used to represent the list of factors of a number.
  # If, for a number "n", "k" is present in the table "f" of its factors,
  # "f[k]" contains the exponent of "k" in the prime factor decomposition.
  Factors = Table[int, int]


####################################################################################################
# Miscellaneous.

## Parity tests.
template isOdd(n: int): bool = (n and 1) != 0
template isEven(n: int): bool = (n and 1) == 0

#---------------------------------------------------------------------------------------------------

proc sort(poly: var Polynomial) {.inline.} =
  ## Sort procedure for the terms of a polynomial (high degree first).
  algorithm.sort(poly, proc(x, y: Term): int = cmp(x.exp, y.exp), Descending)


####################################################################################################
# Superscripts.

const Superscripts: array['0'..'9', string] = ["⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹"]

func superscript(n: Natural): string =
  ## Return the Unicode string to use to represent an exponent.
  if n == 1:
    return ""
  for d in $n:
    result.add(Superscripts[d])


####################################################################################################
# Term operations.

func term(coeff, exp: int): Term =
  ## Create a term.
  if exp < 0:
    raise newException(ValueError, "term exponent cannot be negative")
  (coeff, Natural exp)

#---------------------------------------------------------------------------------------------------

func `*`(a, b: Term): Term =
  ## Multiply two terms.
  (a.coeff * b.coeff, Natural a.exp + b.exp)

#---------------------------------------------------------------------------------------------------

func `+`(a, b: Term): Term =
  ## Add two terms.

  if a.exp != b.exp:
    raise newException(ValueError, "addition of terms with unequal exponents")
  (a.coeff + b.coeff, a.exp)

#---------------------------------------------------------------------------------------------------

func `-`(a: Term): Term =
  ## Return the opposite of a term.
  (-a.coeff, a.exp)

#---------------------------------------------------------------------------------------------------

func `$`(a: Term): string =
  ## Return the string representation of a term.
  if a.coeff == 0: "0"
  elif a.exp == 0: $a.coeff
  elif a.coeff == 1: 'x' & superscript(a.exp)
  elif a.coeff == -1: "-x" & superscript(a.exp)
  else: $a.coeff & 'x' & superscript(a.exp)


####################################################################################################
# Polynomial.

func polynomial(terms: varargs[Term]): Polynomial =
  ## Create a polynomial described by its terms.
  for t in terms:
    if t.coeff != 0:
      result.add(t)
  if result.len == 0:
    return @[term(0, 0)]
  sort(result)

#---------------------------------------------------------------------------------------------------

func hasCoeffAbs(poly: Polynomial; coeff: int): bool =
  ## Return true if the polynomial contains a given coefficient.
  for t in poly:
    if abs(t.coeff) == coeff:
      return true

#---------------------------------------------------------------------------------------------------

func leadingCoeff(poly: Polynomial): int {.inline.} =
  ## Return the coefficient of the term with the highest degree.
  poly[0].coeff

#---------------------------------------------------------------------------------------------------

func degree(poly: Polynomial): int {.inline.} =
  ## Return the degree of the polynomial.
  if poly.len == 0: -1
  else: poly[0].exp

#---------------------------------------------------------------------------------------------------

func `+`(poly: Polynomial; someTerm: Term): Polynomial =
  ## Add a term to a polynomial.

  var added = false
  for currTerm in poly:
    if currterm.exp == someTerm.exp:
      added = true
      if currTerm.coeff + someTerm.coeff != 0:
        result.add(currTerm + someTerm)
    else:
      result.add(currTerm)

  if not added:
    result.add(someTerm)

#---------------------------------------------------------------------------------------------------

func `+`(a, b: Polynomial): Polynomial =
  ## Add two polynomials.

  var aIndex = a.high
  var bIndex = b.high

  while aIndex >= 0 or bIndex >= 0:
    if aIndex < 0:
      result &= b[bIndex]
      dec bIndex
    elif bIndex < 0:
      result &= a[aIndex]
      dec aIndex
    else:
      let t1 = a[aIndex]
      let t2 = b[bIndex]
      if t1.exp == t2.exp:
        let t3 = t1 + t2
        if t3.coeff != 0:
          result.add(t3)
        dec aIndex
        dec bIndex
      elif t1.exp < t2.exp:
        result.add(t1)
        dec aIndex
      else:
        result.add(t2)
        dec bIndex

  sort(result)

#---------------------------------------------------------------------------------------------------

func `*`(poly: Polynomial; someTerm: Term): Polynomial =
  ## Multiply a polynomial by a term.
  for currTerm in poly:
    result.add(currTerm * someTerm)

#---------------------------------------------------------------------------------------------------

func `/`(a, b: Polynomial): Polynomial =
  ## Divide a polynomial by another polynomial.

  var a = a
  let lcb = b.leadingCoeff
  let db = b.degree
  while a.degree >= b.degree:
    let lca = a.leadingCoeff
    let s = lca div lcb
    let t = term(s, a.degree - db)
    result = result + t
    a = a + b * -t

#---------------------------------------------------------------------------------------------------

func `$`(poly: Polynomial): string =
  ## Return the string representation of a polynomial.

  for t in poly:
    if result.len == 0:
      result.add($t)
    else:
      if t.coeff > 0:
        result.add('+')
        result.add($t)
      else:
        result.add('-')
        result.add($(-t))


####################################################################################################
# Cyclotomic polynomial.

var

  # Cache of list of factors.
  factorCache: Table[int, Factors] = {2: {2: 1}.toTable}.toTable

  # Cache of cyclotomic polynomials. Initialized with 1 -> x - 1.
  polyCache: Table[int, Polynomial] = {1: polynomial(term(1, 1), term(-1, 0))}.toTable

#---------------------------------------------------------------------------------------------------

proc getFactors(n: int): Factors =
  ## Return the list of factors of a number.

  if n in factorCache:
    return factorCache[n]

  if n.isEven:
    result = getFactors(n div 2)
    result[2] = result.getOrDefault(2) + 1
    factorCache[n] = result
    return

  var i = 3
  while i * i <= n:
    if n mod i == 0:
      result = getFactors( n div i)
      result[i] = result.getOrDefault(i) + 1
      factorCache[n] = result
      return
    inc i, 2

  result[n] = 1
  factorCache[n] = result

#---------------------------------------------------------------------------------------------------

proc cycloPoly(n: int): Polynomial =
  ## Find the nth cyclotomic polynomial.

  if n in polyCache:
    return polyCache[n]

  let factors = getFactors(n)

  if n in factors:
    # n is prime.
    for i in countdown(n - 1, 0):       # Add the terms by decreasing degrees.
      result.add(term(1, i))

  elif factors.len == 2 and factors.getOrDefault(2) == 1 and factors.getOrDefault(n div 2) == 1:
    # n = 2 x prime.
    let prime = n div 2
    var coeff = -1
    for i in countdown(prime - 1, 0):   # Add the terms by decreasing degrees.
      coeff *= -1
      result.add(term(coeff, i))

  elif factors.len == 1 and 2 in factors:
    # n = 2 ^ h.
    let h = factors[2]
    result.add([term(1, 1 shl (h - 1)), term(1, 0)])

  elif factors.len == 1 and n notin factors:
    # n = prime ^ k.
    var p, k = 0
    for prime, v in factors.pairs:
      if prime > p:
        p = prime
        k = v
    for i in countdown(p - 1, 0):       # Add the terms by decreasing degrees.
      result.add(term(1, i * p^(k-1)))

  elif factors.len == 2 and 2 in factors:
    # n = 2 ^ h x prime ^ k.
    var p, k = 0
    for prime, v in factors.pairs:
      if prime != 2 and prime > p:
        p = prime
        k = v
    var coeff = -1
    let twoExp = 1 shl (factors[2] - 1)
    for i in countdown(p - 1, 0):       # Add the terms by decreasing degrees.
      coeff *= -1
      result.add(term(coeff, i * twoExp * p^(k-1)))

  elif 2 in factors and isOdd(n div 2) and n div 2 > 1:
    # CP(2m)[x] = CP(-m)[x], n odd integer > 1.
    let cycloDiv2 = cycloPoly(n div 2)
    for t in cycloDiv2:
      result.add(if t.exp.isEven: t else: -t)

  else:
    # 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].
    var m = 1
    var cyclo = cycloPoly(m)
    let primes = sorted(toSeq(factors.keys))
    for prime in primes:
      # Compute CP(m)[x^p].
      var terms: Polynomial
      for t in cyclo:
        terms.add(term(t.coeff, t.exp * prime))
      cyclo = terms / cyclo
      m *= prime
    # Now, m is the largest square free divisor of n.
    let s = n div m
    # Compute CP(n)[x] = CP(m)[x^s].
    for t in cyclo:
      result.add(term(t.coeff, t.exp * s))

  polyCache[n] = result


#———————————————————————————————————————————————————————————————————————————————————————————————————

echo "Cyclotomic polynomials for n ⩽ 30:"
for i in 1..30:
  echo &"Φ{'(' & $i & ')':4} = {cycloPoly(i)}"

echo ""
echo "Smallest cyclotomic polynomial with n or -n as a coefficient:"
var n = 0
for i in 1..10:
  while true:
    inc n
    if cycloPoly(n).hasCoeffAbs(i):
      echo &"Φ{'(' & $n & ')':7} has coefficient with magnitude = {i}"
      dec n
      break
Output:

The program runs in 41 seconds on our reasonably performing laptop.

Cyclotomic polynomials for n ⩽ 30:
Φ(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 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
Φ(6545)  has coefficient with magnitude = 8
Φ(6545)  has coefficient with magnitude = 9
Φ(10465) has coefficient with magnitude = 10

PARI/GP[edit]

Cyclotomic polynomials are a built-in function.

for(n=1,30,print(n," : ",polcyclo(n)))

contains_coeff(n, d) = p=polcyclo(n);for(k=0,poldegree(p),if(abs(polcoef(p,k))==d,return(1)));return(0)

for(d=1,10,i=1; while(contains_coeff(i,d)==0,i=i+1);print(d," : ",i))
Output:

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 1 : 1 2 : 105 3 : 385 4 : 1365 5 : 1785 6 : 2805 7 : 3135 8 : 6545 9 : 6545 10 : 10465

Perl[edit]

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.

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++;
    }
}
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[edit]

Translation of: Julia

Uses several routines from Polynomial_long_division#Phix, tweaked slightly to check remainder is zero and trim the quotient.

-- demo\rosetta\Cyclotomic_Polynomial.exw
with javascript_semantics
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
            integer mi = -i
            n[mi] -= d2[mi]*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 the 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,deep_copy(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<iff(platform()=JS?5: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 and platform()!=JS then
        printf(1,"working (%d) ...\r",n)
        t1 = time()+1
    end if
end while
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[edit]

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
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[edit]

(formerly Perl 6)

Works with: Rakudo version 2020.01
Translation of: Perl

Uses the same library as Perl, so comes with the same caveats.

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)
}
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[edit]

Built-in:

say "First 30 cyclotomic polynomials:"
for k in (1..30) {
    say ("Φ(#{k}) = ", cyclotomic(k))
}

say "\nSmallest cyclotomic polynomial with n or -n as a coefficient:"
for n in (1..10) {  # slow
    var k = (1..Inf -> first {|k|
        cyclotomic(k).coeffs.any { .tail.abs == n }
    })
    say "Φ(#{k}) has coefficient with magnitude #{n}"
}
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

Visual Basic .NET[edit]

Translation of: C++
Imports System.Text

Module Module1
    Private ReadOnly MAX_ALL_FACTORS As Integer = 100_000
#Const ALGORITHM = 2

    Class Term
        Implements IComparable(Of Term)

        Public ReadOnly Property Coefficient As Long
        Public ReadOnly Property Exponent As Long

        Public Sub New(c As Long, Optional e As Long = 0)
            Coefficient = c
            Exponent = e
        End Sub

        Public Shared Operator -(t As Term) As Term
            Return New Term(-t.Coefficient, t.Exponent)
        End Operator

        Public Shared Operator +(lhs As Term, rhs As Term) As Term
            If lhs.Exponent <> rhs.Exponent Then
                Throw New ArgumentException("Exponents not equal")
            End If
            Return New Term(lhs.Coefficient + rhs.Coefficient, lhs.Exponent)
        End Operator

        Public Shared Operator *(lhs As Term, rhs As Term) As Term
            Return New Term(lhs.Coefficient * rhs.Coefficient, lhs.Exponent + rhs.Exponent)
        End Operator

        Public Function CompareTo(other As Term) As Integer Implements IComparable(Of Term).CompareTo
            Return -Exponent.CompareTo(other.Exponent)
        End Function

        Public Overrides Function ToString() As String
            If Coefficient = 0 Then
                Return "0"
            End If
            If Exponent = 0 Then
                Return Coefficient.ToString
            End If
            If Coefficient = 1 Then
                If Exponent = 1 Then
                    Return "x"
                End If
                Return String.Format("x^{0}", Exponent)
            End If
            If Coefficient = -1 Then
                If Exponent = 1 Then
                    Return "-x"
                End If
                Return String.Format("-x^{0}", Exponent)
            End If
            If Exponent = 1 Then
                Return String.Format("{0}x", Coefficient)
            End If
            Return String.Format("{0}x^{1}", Coefficient, Exponent)
        End Function
    End Class

    Class Polynomial
        Implements IEnumerable(Of Term)

        Private ReadOnly polyTerms As New List(Of Term)

        Public Sub New()
            polyTerms.Add(New Term(0))
        End Sub

        Public Sub New(ParamArray values() As Term)
            If values.Length = 0 Then
                polyTerms.Add(New Term(0))
            Else
                polyTerms.AddRange(values)
            End If
            Normalize()
        End Sub

        Public Sub New(values As IEnumerable(Of Term))
            polyTerms.AddRange(values)
            If polyTerms.Count = 0 Then
                polyTerms.Add(New Term(0))
            End If
            Normalize()
        End Sub

        Public Function LeadingCoeficient() As Long
            Return polyTerms(0).Coefficient
        End Function

        Public Function Degree() As Long
            Return polyTerms(0).Exponent
        End Function

        Public Function HasCoefficentAbs(coeff As Long) As Boolean
            For Each t In polyTerms
                If Math.Abs(t.Coefficient) = coeff Then
                    Return True
                End If
            Next
            Return False
        End Function

        Public Function GetEnumerator() As IEnumerator(Of Term) Implements IEnumerable(Of Term).GetEnumerator
            Return polyTerms.GetEnumerator
        End Function

        Private Function IEnumerable_GetEnumerator() As IEnumerator Implements IEnumerable.GetEnumerator
            Return polyTerms.GetEnumerator
        End Function

        Private Sub Normalize()
            polyTerms.Sort(Function(a As Term, b As Term) a.CompareTo(b))
        End Sub

        Public Shared Operator +(lhs As Polynomial, rhs As Term) As Polynomial
            Dim terms As New List(Of Term)
            Dim added = False
            For Each ct In lhs
                If ct.Exponent = rhs.Exponent Then
                    added = True
                    If ct.Coefficient + rhs.Coefficient <> 0 Then
                        terms.Add(ct + rhs)
                    End If
                Else
                    terms.Add(ct)
                End If
            Next
            If Not added Then
                terms.Add(rhs)
            End If
            Return New Polynomial(terms)
        End Operator

        Public Shared Operator *(lhs As Polynomial, rhs As Term) As Polynomial
            Dim terms As New List(Of Term)
            For Each ct In lhs
                terms.Add(ct * rhs)
            Next
            Return New Polynomial(terms)
        End Operator

        Public Shared Operator +(lhs As Polynomial, rhs As Polynomial) As Polynomial
            Dim terms As New List(Of Term)
            Dim thisCount = lhs.polyTerms.Count
            Dim polyCount = rhs.polyTerms.Count
            While thisCount > 0 OrElse polyCount > 0
                If thisCount = 0 Then
                    Dim polyTerm = rhs.polyTerms(polyCount - 1)
                    terms.Add(polyTerm)
                    polyCount -= 1
                ElseIf polyCount = 0 Then
                    Dim thisTerm = lhs.polyTerms(thisCount - 1)
                    terms.Add(thisTerm)
                    thisCount -= 1
                Else
                    Dim polyTerm = rhs.polyTerms(polyCount - 1)
                    Dim thisTerm = lhs.polyTerms(thisCount - 1)
                    If thisTerm.Exponent = polyTerm.Exponent Then
                        Dim t = thisTerm + polyTerm
                        If t.Coefficient <> 0 Then
                            terms.Add(t)
                        End If
                        thisCount -= 1
                        polyCount -= 1
                    ElseIf thisTerm.Exponent < polyTerm.Exponent Then
                        terms.Add(thisTerm)
                        thisCount -= 1
                    Else
                        terms.Add(polyTerm)
                        polyCount -= 1
                    End If
                End If
            End While
            Return New Polynomial(terms)
        End Operator

        Public Shared Operator *(lhs As Polynomial, rhs As Polynomial) As Polynomial
            Throw New Exception("Not implemented")
        End Operator

        Public Shared Operator /(lhs As Polynomial, rhs As Polynomial) As Polynomial
            Dim q As New Polynomial
            Dim r = lhs
            Dim lcv = rhs.LeadingCoeficient
            Dim dv = rhs.Degree
            While r.Degree >= rhs.Degree
                Dim lcr = r.LeadingCoeficient
                Dim s = lcr \ lcv
                Dim t As New Term(s, r.Degree() - dv)
                q += t
                r += rhs * -t
            End While
            Return q
        End Operator

        Public Overrides Function ToString() As String
            Dim builder As New StringBuilder
            Dim it = polyTerms.GetEnumerator()
            If it.MoveNext Then
                builder.Append(it.Current)
            End If
            While it.MoveNext
                If it.Current.Coefficient < 0 Then
                    builder.Append(" - ")
                    builder.Append(-it.Current)
                Else
                    builder.Append(" + ")
                    builder.Append(it.Current)
                End If
            End While
            Return builder.ToString
        End Function
    End Class

    Function GetDivisors(number As Integer) As List(Of Integer)
        Dim divisors As New List(Of Integer)
        Dim root = CType(Math.Sqrt(number), Long)
        For i = 1 To root
            If number Mod i = 0 Then
                divisors.Add(i)
                Dim div = number \ i
                If div <> i AndAlso div <> number Then
                    divisors.Add(div)
                End If
            End If
        Next
        Return divisors
    End Function

    Private ReadOnly allFactors As New Dictionary(Of Integer, Dictionary(Of Integer, Integer)) From {{2, New Dictionary(Of Integer, Integer) From {{2, 1}}}}
    Function GetFactors(number As Integer) As Dictionary(Of Integer, Integer)
        If allFactors.ContainsKey(number) Then
            Return allFactors(number)
        End If

        Dim factors As New Dictionary(Of Integer, Integer)
        If number Mod 2 = 0 Then
            Dim factorsDivTwo = GetFactors(number \ 2)
            For Each pair In factorsDivTwo
                If Not factors.ContainsKey(pair.Key) Then
                    factors.Add(pair.Key, pair.Value)
                End If
            Next
            If factors.ContainsKey(2) Then
                factors(2) += 1
            Else
                factors.Add(2, 1)
            End If
            If number < MAX_ALL_FACTORS Then
                allFactors.Add(number, factors)
            End If
            Return factors
        End If
        Dim root = CType(Math.Sqrt(number), Long)
        Dim i = 3L
        While i <= root
            If number Mod i = 0 Then
                Dim factorsDivI = GetFactors(number \ i)
                For Each pair In factorsDivI
                    If Not factors.ContainsKey(pair.Key) Then
                        factors.Add(pair.Key, pair.Value)
                    End If
                Next
                If factors.ContainsKey(i) Then
                    factors(i) += 1
                Else
                    factors.Add(i, 1)
                End If
                If number < MAX_ALL_FACTORS Then
                    allFactors.Add(number, factors)
                End If
                Return factors
            End If
            i += 2
        End While
        factors.Add(number, 1)
        If number < MAX_ALL_FACTORS Then
            allFactors.Add(number, factors)
        End If
        Return factors
    End Function

    Private ReadOnly computedPolynomials As New Dictionary(Of Integer, Polynomial)
    Function CyclotomicPolynomial(n As Integer) As Polynomial
        If computedPolynomials.ContainsKey(n) Then
            Return computedPolynomials(n)
        End If

        If n = 1 Then
            REM polynomial: x - 1
            Dim p As New Polynomial(New Term(1, 1), New Term(-1))
            computedPolynomials.Add(n, p)
            Return p
        End If

        Dim factors = GetFactors(n)
        Dim terms As New List(Of Term)
        Dim cyclo As Polynomial

        If factors.ContainsKey(n) Then
            REM n prime
            For index = 1 To n
                terms.Add(New Term(1, index - 1))
            Next

            cyclo = New Polynomial(terms)
            computedPolynomials.Add(n, cyclo)
            Return cyclo
        ElseIf factors.Count = 2 AndAlso factors.ContainsKey(2) AndAlso factors(2) = 1 AndAlso factors.ContainsKey(n / 2) AndAlso factors(n / 2) = 1 Then
            REM n = 2p
            Dim prime = n \ 2
            Dim coeff = -1

            For index = 1 To prime
                coeff *= -1
                terms.Add(New Term(coeff, index - 1))
            Next

            cyclo = New Polynomial(terms)
            computedPolynomials.Add(n, cyclo)
            Return cyclo
        ElseIf factors.Count = 1 AndAlso factors.ContainsKey(2) Then
            REM n = 2^h
            Dim h = factors(2)
            terms = New List(Of Term) From {
                New Term(1, Math.Pow(2, h - 1)),
                New Term(1)
            }

            cyclo = New Polynomial(terms)
            computedPolynomials.Add(n, cyclo)
            Return cyclo
        ElseIf factors.Count = 1 AndAlso factors.ContainsKey(n) Then
            REM n = p^k
            Dim p = 0
            Dim k = 0
            For Each it In factors
                p = it.Key
                k = it.Value
            Next
            For index = 1 To p
                terms.Add(New Term(1, (index - 1) * Math.Pow(p, k - 1)))
            Next

            cyclo = New Polynomial(terms)
            computedPolynomials.Add(n, cyclo)
            Return cyclo
        ElseIf factors.Count = 2 AndAlso factors.ContainsKey(2) Then
            REM n = 2^h * p^k
            Dim p = 0
            For Each it In factors
                If it.Key <> 2 Then
                    p = it.Key
                End If
            Next

            Dim coeff = -1
            Dim twoExp = CType(Math.Pow(2, factors(2) - 1), Long)
            Dim k = factors(p)
            For index = 1 To p
                coeff *= -1
                terms.Add(New Term(coeff, (index - 1) * twoExp * Math.Pow(p, k - 1)))
            Next

            cyclo = New Polynomial(terms)
            computedPolynomials.Add(n, cyclo)
            Return cyclo
        ElseIf factors.ContainsKey(2) AndAlso (n / 2) Mod 2 = 1 AndAlso n / 2 > 1 Then
            REM CP(2m)[x] = CP(-m)[x], n odd integer > 1
            Dim cycloDiv2 = CyclotomicPolynomial(n \ 2)
            For Each t In cycloDiv2
                If t.Exponent Mod 2 = 0 Then
                    terms.Add(t)
                Else
                    terms.Add(-t)
                End If
            Next

            cyclo = New Polynomial(terms)
            computedPolynomials.Add(n, cyclo)
            Return cyclo
        End If

#If ALGORITHM = 0 Then
        REM slow - uses basic definition
        Dim divisors = GetDivisors(n)
        REM Polynomial: (x^n - 1)
        cyclo = New Polynomial(New Term(1, n), New Term(-1))
        For Each i In divisors
            Dim p = CyclotomicPolynomial(i)
            cyclo /= p
        Next

        computedPolynomials.Add(n, cyclo)
        Return cyclo
#ElseIf ALGORITHM = 1 Then
        REM Faster.  Remove Max divisor (and all divisors of max divisor) - only one divide for all divisors of Max Divisor
        Dim divisors = GetDivisors(n)
        Dim maxDivisor = Integer.MinValue
        For Each div In divisors
            maxDivisor = Math.Max(maxDivisor, div)
        Next
        Dim divisorExceptMax As New List(Of Integer)
        For Each div In divisors
            If maxDivisor Mod div <> 0 Then
                divisorExceptMax.Add(div)
            End If
        Next

        REM Polynomial:  ( x^n - 1 ) / ( x^m - 1 ), where m is the max divisor
        cyclo = New Polynomial(New Term(1, n), New Term(-1)) / New Polynomial(New Term(1, maxDivisor), New Term(-1))
        For Each i In divisorExceptMax
            Dim p = CyclotomicPolynomial(i)
            cyclo /= p
        Next

        computedPolynomials.Add(n, cyclo)
        Return cyclo
#ElseIf ALGORITHM = 2 Then
        REM Fastest
        REM Let p ; q be primes such that p does not divide n, and q divides n
        REM Then Cp(np)[x] = CP(n)[x^p] / CP(n)[x]
        Dim m = 1
        cyclo = CyclotomicPolynomial(m)
        Dim primes As New List(Of Integer)
        For Each it In factors
            primes.Add(it.Key)
        Next
        primes.Sort()
        For Each prime In primes
            REM CP(m)[x]
            Dim cycloM = cyclo
            REM Compute CP(m)[x^p]
            terms = New List(Of Term)
            For Each t In cyclo
                terms.Add(New Term(t.Coefficient, t.Exponent * prime))
            Next
            cyclo = New Polynomial(terms) / cycloM
            m *= prime
        Next
        REM Now, m is the largest square free divisor of n
        Dim s = n \ m
        REM Compute CP(n)[x] = CP(m)[x^s]
        terms = New List(Of Term)
        For Each t In cyclo
            terms.Add(New Term(t.Coefficient, t.Exponent * s))
        Next

        cyclo = New Polynomial(terms)
        computedPolynomials.Add(n, cyclo)
        Return cyclo
#Else
        Throw New Exception("Invalid algorithm")
#End If
    End Function

    Sub Main()
        Console.WriteLine("Task 1:  cyclotomic polynomials for n <= 30:")
        For i = 1 To 30
            Dim p = CyclotomicPolynomial(i)
            Console.WriteLine("CP[{0}] = {1}", i, p)
        Next
        Console.WriteLine()

        Console.WriteLine("Task 2:  Smallest cyclotomic polynomial with n or -n as a coefficient:")
        Dim n = 0
        For i = 1 To 10
            While True
                n += 1
                Dim cyclo = CyclotomicPolynomial(n)
                If cyclo.HasCoefficentAbs(i) Then
                    Console.WriteLine("CP[{0}] has coefficient with magnitude = {1}", n, i)
                    n -= 1
                    Exit While
                End If
            End While
        Next
    End Sub

End Module
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

Wren[edit]

Translation of: Go
Library: Wren-trait
Library: Wren-sort
Library: Wren-math
Library: Wren-fmt

Second part is very slow. Limited to first 7 to finish in a reasonable time - 5 minutes on my machine.

import "/trait" for Stepped
import "/sort" for Sort
import "/math" for Int, Nums
import "/fmt" for Fmt

var algo = 2
var maxAllFactors = 1e5

class Term {
    construct new(coef, exp) {
        _coef = coef
        _exp = exp
    }

    coef { _coef }
    exp  { _exp }

    *(t) { Term.new(_coef * t.coef, _exp + t.exp) }

    +(t) {
        if (_exp != t.exp) Fiber.abort("Exponents unequal in term '+' method.")
        return Term.new(_coef + t.coef, _exp)
    }

    - { Term.new(-_coef, _exp) }

    toString {
        if (_coef == 0) return "0"
        if (_exp == 0)  return _coef.toString
        if (_coef == 1) return (_exp == 1) ? "x" : "x^%(_exp)"
        if (_exp == 1)  return "%(_coef)x"
        return "%(_coef)x^%(_exp)"
    }
}

class Poly {
    // pass coef, exp in pairs as parameters
    construct new(values) {
        var le = values.count
        if (le == 0) {
            _terms = [Term.new(0, 0)]
        } else {
            if (le%2 != 0) Fiber.abort("Odd number of parameters(%(le)) passed to Poly constructor.")
            _terms = []
            for (i in Stepped.new(0...le, 2)) _terms.add(Term.new(values[i], values[i+1]))
            tidy()
        }
    }

    terms { _terms }

    hasCoefAbs(coef) { _terms.any { |t| t.coef.abs == coef } }

    +(p2) {
        var p3 = Poly.new([])
        var le = _terms.count
        var le2 = p2.terms.count
        while (le > 0 || le2 > 0) {
            if (le == 0) {
                p3.terms.add(p2.terms[le2-1])
                le2 = le2 - 1
            } else if (le2 == 0) {
                p3.terms.add(_terms[le-1])
                le = le - 1
            } else {
                var t = _terms[le-1]
                var t2 = p2.terms[le2-1]
                if (t.exp == t2.exp) {
                    var t3 = t + t2
                    if (t3.coef != 0) p3.terms.add(t3)
                    le = le - 1
                    le2 = le2 - 1
                } else if (t.exp < t2.exp) {
                    p3.terms.add(t)
                    le = le - 1
                } else {
                    p3.terms.add(t2)
                    le2 = le2 - 1
                }
            }
        }
        p3.tidy()
        return p3
    }

    addTerm(t) {
        var q = Poly.new([])
        var added = false
        for (i in 0..._terms.count) {
            var ct = _terms[i]
            if (ct.exp == t.exp) {
                added = true
                if (ct.coef + t.coef != 0) q.terms.add(ct + t)
            } else {
                q.terms.add(ct)
            }
        }
        if (!added) q.terms.add(t)
        q.tidy()
        return q
    }

    mulTerm(t) {
        var q = Poly.new([])
        for (i in 0..._terms.count) {
            var ct = _terms[i]
            q.terms.add(ct * t)
        }
        q.tidy()
        return q
    }

    /(v) {
        var p = this
        var q = Poly.new([])
        var lcv = v.leadingCoef
        var dv = v.degree
        while (p.degree >= v.degree) {
            var lcp = p.leadingCoef
            var s = (lcp/lcv).truncate
            var t = Term.new(s, p.degree - dv)
            q = q.addTerm(t)
            p = p + v.mulTerm(-t) 
        }
        q.tidy()
        return q
    }

    leadingCoef { _terms[0].coef }

    degree { _terms[0].exp }

    toString {
        var sb = ""
        var first = true
        for (t in _terms) {
            if (first) {
                sb = sb + t.toString
                first = false
            } else {
                sb = sb + " "
                if (t.coef > 0) {
                    sb = sb + "+ "
                    sb = sb + t.toString
                } else {
                    sb = sb + "- "
                    sb = sb + (-t).toString
                }
            }
        }
        return sb
    }

    // in place descending sort by term.exp
    sortTerms() {
        var cmp  = Fn.new { |t1, t2| (t2.exp - t1.exp).sign }
        Sort.quick(_terms, 0, _terms.count-1, cmp)
    }

    // sort terms and remove any unnecesary zero terms
    tidy() {
        sortTerms()
        if (degree > 0) {
            for (i in _terms.count-1..0) {
                if (_terms[i].coef == 0) _terms.removeAt(i)
            }
            if (_terms.count == 0) _terms.add(Term.new(0, 0))
        }
    }
}

var computed = {}
var allFactors = {2: {2: 1}}

var getFactors // recursive function
getFactors = Fn.new { |n|
    var f = allFactors[n]
    if (f) return f
    var factors = {}
    if (n%2 == 0) {
        var factorsDivTwo = getFactors.call(n/2)
        for (me in factorsDivTwo) factors[me.key] = me.value
        factors[2] = factors[2] ? factors[2] + 1 : 1
        if (n < maxAllFactors) allFactors[n] = factors
        return factors
    }
    var prime = true
    var sqrt = n.sqrt.floor
    var i = 3
    while (i <= sqrt){
        if (n%i == 0) {
            prime = false
            for (me in getFactors.call(n/i)) factors[me.key] = me.value
            factors[i] = factors[i] ? factors[i] + 1 : 1
            if (n < maxAllFactors) allFactors[n] = factors
            return factors
        }
        i = i + 2
    }
    if (prime) {
        factors[n] = 1
        if (n < maxAllFactors) allFactors[n] = factors
    }
    return factors
}

var cycloPoly // recursive function
cycloPoly = Fn.new { |n|
    var p = computed[n]
    if (p) return p
    if (n == 1) {
        // polynomialL x - 1
        p = Poly.new([1, 1, -1, 0])
        computed[1] = p
        return p
    }
    var factors = getFactors.call(n)
    var cyclo = Poly.new([])
    if (factors[n]) {
        // n is prime
        for (i in 0...n) cyclo.terms.add(Term.new(1, i))
    } else if (factors.count == 2 && factors[2] == 1 && factors[n/2] == 1) {
        // n == 2p
        var prime = n / 2
        var coef = -1
        for (i in 0...prime) {
            coef = coef * (-1)
            cyclo.terms.add(Term.new<