Cyclotomic polynomial: Difference between revisions

(Added C#)
(→‎{{header|jq}}: task2(10))
 
(52 intermediate revisions by 13 users not shown)
Line 14:
=={{header|C++}}==
{{trans|Java}}
<langsyntaxhighlight lang="cpp">#include <algorithm>
#include <iostream>
#include <initializer_list>
Line 20:
#include <vector>
 
const int MAX_ALL_FACTORS = 100'000100000;
const int algorithm = 2;
int divisions = 0;
Line 532:
 
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>Task 1: cyclotomic polynomials for n <= 30:
Line 578:
CP[10465] has coefficient with magnitude = 10</pre>
 
=={{header|C sharp|C#}}==
{{trans|Java}}
{{works with|C sharp|8}}
<langsyntaxhighlight lang="csharp">using System;
using System.Collections;
using System.Collections.Generic;
Line 816:
};
}
}</langsyntaxhighlight>
{{out}}
<pre>
Line 862:
CP[6545] has coefficient with magnitude = 9
CP[10465] has coefficient with magnitude = 10</pre>
 
=={{header|D}}==
{{trans|Kotlin}}
<syntaxhighlight lang="d">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;
}
}
}
}</syntaxhighlight>
{{out}}
<pre>Task 1: cyclotomic polynomials for n <= 30:
CP[1] = x - 1
CP[2] = x + 1
CP[3] = x^2 + x + 1
CP[4] = x^2 + 1
CP[5] = x^4 + x^3 + x^2 + x + 1
CP[6] = x^2 - x + 1
CP[7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[8] = x^4 + 1
CP[9] = x^6 + x^3 + 1
CP[10] = x^4 - x^3 + x^2 - x + 1
CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[12] = x^4 - x^2 + 1
CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1
CP[16] = x^8 + 1
CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[18] = x^6 - x^3 + 1
CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[20] = x^8 - x^6 + x^4 - x^2 + 1
CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1
CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[24] = x^8 - x^4 + 1
CP[25] = x^20 + x^15 + x^10 + x^5 + 1
CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[27] = x^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</pre>
 
=={{header|Fermat}}==
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.
<syntaxhighlight lang="fermat">
&(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;</syntaxhighlight>
 
=={{header|Go}}==
{{trans|Java}}
<langsyntaxhighlight lang="go">package main
 
import (
Line 1,304 ⟶ 1,824:
}
}
}</langsyntaxhighlight>
 
{{out}}
Line 1,352 ⟶ 1,872:
CP[10465] has coefficient with magnitude = 10
</pre>
<math>Insert formula here</math>
 
=={{header|Haskell}}==
Uses synthetic polynomial division and simple memoization.
 
<syntaxhighlight lang="haskell">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))</syntaxhighlight>
 
Simple examples
 
<pre>λ> 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]</pre>
 
The task solution
 
<syntaxhighlight lang="haskell">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</syntaxhighlight>
 
Result
 
<pre>"-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</pre>
 
Computations take a while...
 
=={{header|J}}==
 
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:
 
<syntaxhighlight lang="j">cyclo=: {{<.-:1+(++) p. 1;^0j2p1* y%~1+I.1=y+.1+i.y}}</syntaxhighlight>
 
This approach suggests that cyclotomic polynomial zero should be <tt>f<sub>0</sub>(x)= 1</tt>
 
Routine to find the nth cyclotomic polynomial:
 
<syntaxhighlight lang="j">{{ 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
}}</syntaxhighlight>
 
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, x<sup>2</sup>+x+1, x<sup>2</sup>+1, x<sup>2</sup>-x+1, and x<sup>4</sup>-x<sup>2</sup>+1) you get x<sup>n</sup>-1 (for 12, that would of course be x<sup>12</sup>-1).
 
Notes:
 
* the coefficients of cyclotomic polynomials after 1 form a palindrome (that's the <tt>q#1</tt> 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 <tt>,(y%*/q) {."0</tt> ... 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 <tt>(* 1 _1 $~ #)</tt> phrase in the implementation.
* To deal with the general case, we use polynomial division, x<sup>n</sup>-1 divided by the polynomial product the cyclotomic polynomials of the proper divisors of number we're looking for.
* <tt>+//.@(*/)</tt> is polynomial product in J.
 
Task examples:
 
<syntaxhighlight lang="j">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</syntaxhighlight>
 
=== Another approach ===
 
As noted in the [http://jsoftware.com/pipermail/programming/2022-March/060209.html J programming forum], we can improve the big-O character of this algorithm by using the [[Fast Fourier transform#J|fast fourier transform]] for polynomial multiplication and division.
 
<syntaxhighlight lang="j">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.</syntaxhighlight>
 
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 [[j:Essays/FFT|fft for polynomial multiplication]]
 
This approach gave slightly over a 16x speedup for <tt>taskorder 10</tt>, 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.)
 
=={{header|Java}}==
<langsyntaxhighlight lang="java">
import java.util.ArrayList;
import java.util.Collections;
Line 1,881 ⟶ 2,703:
 
}
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 1,927 ⟶ 2,749:
CP[6545] has coefficient with magnitude = 9
CP[10465] has coefficient with magnitude = 10
</pre>
 
===Alternative Version===
An alternative example using the algorithm from "Matters Computational" by Jorg Arndt, pages 704 - 705.
It completes the task in less than 2 seconds.
<syntaxhighlight lang="java">
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
 
public class CyclotomicPolynomials {
 
public static void main(String[] args) {
System.out.println("Task 1: Cyclotomic polynomials for n <= 30:");
for ( int cpIndex = 1; cpIndex <= 30; cpIndex++ ) {
System.out.println("CP[" + cpIndex + "] = " + toString(cyclotomicPolynomial(cpIndex)));
}
System.out.println();
System.out.println("Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:");
System.out.println("CP[1] has a coefficient with magnitude 1");
int cpIndex = 2;
for ( int coefficient = 2; coefficient <= 10; coefficient++ ) {
while ( BigInteger.valueOf(cpIndex).isProbablePrime(PRIME_CERTAINTY)
|| ! hasHeight(cyclotomicPolynomial(cpIndex), coefficient) ) {
cpIndex += 1;
}
System.out.println("CP[" + cpIndex + "] has a coefficient with magnitude " + coefficient);
}
}
// Return the Cyclotomic Polynomial of order 'cpIndex' as an array of coefficients,
// where, for example, the polynomial 3x^2 - 1 is represented by the array [3, 0, -1].
private static int[] cyclotomicPolynomial(int cpIndex) {
int[] polynomial = new int[] { 1, -1 };
if ( cpIndex == 1 ) {
return polynomial;
}
if ( BigInteger.valueOf(cpIndex).isProbablePrime(PRIME_CERTAINTY) ) {
int[] result = new int[cpIndex];
Arrays.fill(result, 1);
return result;
}
List<Integer> primes = distinctPrimeFactors(cpIndex);
int product = 1;
for ( int prime : primes ) {
int[] numerator = substituteExponent(polynomial, prime);
polynomial = exactDivision(numerator, polynomial);
product *= prime;
}
polynomial = substituteExponent(polynomial, cpIndex / product);
return polynomial;
}
// Return the Cyclotomic Polynomial obtained from 'polynomial' by replacing x with x^'exponent'.
private static int[] substituteExponent(int[] polynomial, int exponent) {
int[] result = new int[exponent * ( polynomial.length - 1 ) + 1];
for ( int i = polynomial.length - 1; i >= 0; i-- ) {
result[i * exponent] = polynomial[i];
}
return result;
}
// Return the Cyclotomic Polynomial equal to 'dividend' / 'divisor'. The division is always exact.
private static int[] exactDivision(int[] dividend, int[] divisor) {
int[] result = Arrays.copyOf(dividend, dividend.length);
for ( int i = 0; i < dividend.length - divisor.length + 1; i++ ) {
if ( result[i] != 0 ) {
for ( int j = 1; j < divisor.length; j ++ ) {
result[i + j] += -divisor[j] * result[i];
}
}
}
result = Arrays.copyOf(result, result.length - divisor.length + 1);
return result;
}
 
// Return whether 'polynomial' has a coefficient of equal magnitude to 'coefficient'.
private static boolean hasHeight(int[] polynomial, int coefficient) {
for ( int i = 0; i <= ( polynomial.length + 1 ) / 2; i++ ) {
if ( Math.abs(polynomial[i]) == coefficient ) {
return true;
}
}
return false;
}
// Return a string representation of 'polynomial'.
private static String toString(int[] polynomial) {
StringBuilder text = new StringBuilder();
for ( int i = 0; i < polynomial.length; i++ ) {
if ( polynomial[i] == 0 ) {
continue;
}
text.append(( polynomial[i] < 0 ) ? ( ( i == 0 ) ? "-" : " - " ) : ( ( i == 0 ) ? "" : " + " ));
final int exponent = polynomial.length - 1 - i;
if ( exponent > 0 && Math.abs(polynomial[i]) > 1 ) {
text.append(Math.abs(polynomial[i]));
}
text.append(( exponent > 1 ) ?
( "x^" + String.valueOf(exponent) ) : ( ( exponent == 1 ) ? "x" : Math.abs(polynomial[i]) ));
}
return text.toString();
}
// Return a list of the distinct prime factors of 'number'.
private static List<Integer> distinctPrimeFactors(int number) {
List<Integer> primeFactors = new ArrayList<Integer>();
for ( int divisor = 2; divisor * divisor <= number; divisor++ ) {
if ( number % divisor == 0 ) {
primeFactors.add(divisor);
}
while ( number % divisor == 0 ) {
number = number / divisor;
}
}
if ( number > 1 ) {
primeFactors.add(number);
}
return primeFactors;
}
private static final int PRIME_CERTAINTY = 15;
 
}
</syntaxhighlight>
<pre>
Task 1: Cyclotomic polynomials for n <= 30:
CP[1] = x - 1
CP[2] = x + 1
CP[3] = x^2 + x + 1
CP[4] = x^2 + 1
CP[5] = x^4 + x^3 + x^2 + x + 1
CP[6] = x^2 - x + 1
CP[7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[8] = x^4 + 1
CP[9] = x^6 + x^3 + 1
CP[10] = x^4 - x^3 + x^2 - x + 1
CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[12] = x^4 - x^2 + 1
CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1
CP[16] = x^8 + 1
CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[18] = x^6 - x^3 + 1
CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[20] = x^8 - x^6 + x^4 - x^2 + 1
CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1
CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[24] = x^8 - x^4 + 1
CP[25] = x^20 + x^15 + x^10 + x^5 + 1
CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[27] = x^18 + x^9 + 1
CP[28] = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1
CP[29] = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[30] = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1
 
Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:
CP[1] has 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
</pre>
 
=={{header|jq}}==
'''Adapted from the [[#Wren|Wren]] implementation of the Arndt algorithm'''
 
'''Works with jq, the C implementation of jq'''
 
'''Works with gojq, the Go implementation of jq'''
 
'''Works with jaq, the Rust implementation of jq'''
 
In this entry, a polynomial of degree n is represented by a JSON
array of length n+1 beginning with the leading coefficient.
 
For the second task, besides exploiting the fact that
CP[1] has height 1, the following program only assumes that the
cyclotomic polynomials are palindromic, but avoids recomputing them.
<syntaxhighlight lang="jq">
### Generic utilities
 
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;
 
# Return the maximum item in the stream assuming it is not empty:
def max(s): reduce s as $s (null; if . == null then $s elif $s > . then $s else . end);
 
# Truncated integer division (consistent with % operator).
# `round` is used for the sake of jaq.
def quo($x; $y): ($x - ($x%$y)) / $y | round;
 
### Primes
def is_prime:
. as $n
| if ($n < 2) then false
elif ($n % 2 == 0) then $n == 2
elif ($n % 3 == 0) then $n == 3
elif ($n % 5 == 0) then $n == 5
elif ($n % 7 == 0) then $n == 7
elif ($n % 11 == 0) then $n == 11
elif ($n % 13 == 0) then $n == 13
elif ($n % 17 == 0) then $n == 17
elif ($n % 19 == 0) then $n == 19
else sqrt as $s
| 23
| until( . > $s or ($n % . == 0); . + 2)
| . > $s
end;
 
# Emit an array of the distinct prime factors of 'n' in order using a wheel
# with basis [2, 3, 5], e.g. 44 | distinctPrimeFactors #=> [2,11]
def distinctPrimeFactors:
def augment($x): if .[-1] == $x then . else . + [$x] end;
def out($i):
if (.n % $i) == 0
then .factors += [$i]
| until (.n % $i != 0; .n = ((.n/$i)|floor) )
else .
end;
if . < 2 then []
else [4, 2, 4, 2, 4, 6, 2, 6] as $inc
| { n: .,
factors: [] }
| out(2)
| out(3)
| out(5)
| .k = 7
| .i = 0
| until(.k * .k > .n;
if .n % .k == 0
then .k as $k | .factors |= augment($k)
| .n = ((.n/.k)|floor)
else .k += $inc[.i]
| .i = ((.i + 1) % 8)
end)
| if .n > 1 then .n as $n | .factors |= augment($n) else . end
| .factors
end;
 
### Polynomials
def canonical:
if length == 0 then .
elif .[-1] == 0 then .[:-1]|canonical
else .
end;
 
# For pretty-printing the input array as the polynomial it represents
# e.g. [1,-1] => x-1
def pp:
def digits: tostring | explode[] | [.] | implode | tonumber;
def superscript:
if . <= 1 then ""
else ["\u2070", "\u00b9", "\u00b2", "\u00b3", "\u2074", "\u2075", "\u2076", "\u2077", "\u2078", "\u2079"] as $ss
| reduce digits as $d (""; . + $ss[$d] )
end;
 
if length == 1 then .[0] | tostring
else reverse as $p
| reduce range(length-1; -1; -1) as $i ([];
if $p[$i] != 0
then (if $i > 0 then "x" else "" end) as $x
| ( if $i > 0 and ($p[$i]|length) == 1
then (if $p[$i] == 1 then "" else "-" end)
else ($p[$i]|tostring)
end ) as $c
| . + ["\($c)\($x)\($i|superscript)"]
else . end )
| join("+")
| gsub("\\+-"; "-")
end ;
 
def polynomialDivide($divisor):
. as $in
| ($divisor|canonical) as $divisor
| { curr: canonical}
| .base = ((.curr|length) - ($divisor|length))
| until( .base < 0;
(.curr[-1] / $divisor[-1]) as $res
| .result += [$res]
| .curr |= .[:-1]
| reduce range (0;$divisor|length-1) as $i (.;
.curr[.base + $i] += (- $res * $divisor[$i]) )
| .base += -1 )
| (.result | reverse) as $quot
| (.curr | canonical) as $rem
| [$quot, $rem];
 
# Call `round` for the sake of jaq
def exactDivision($numerator; $denominator):
($numerator | polynomialDivide($denominator))
| .[0]
| map(round);
 
def init($n; $value): [range(0;$n)|$value];
 
### Cyclotomic Polynomials
 
# The Cyclotomic Polynomial obtained from $polynomial
# by replacing x with x^$exponent
def substituteExponent($polynomial; $exponent):
init( ($polynomial|length - 1) * $exponent + 1; 0)
| reduce range(0; $polynomial|length) as $i (.; .[$i*$exponent] = $polynomial[$i]);
 
# Return the Cyclotomic Polynomial of order 'cpIndex' as a JSON array of coefficients,
# where, for example, the polynomial 3x^2 - 1 is represented by [3, 0, -1].
def cycloPoly($cpIndex):
{ polynomial: [1, -1] }
| if $cpIndex == 1 then .polynomial
elif ($cpIndex|is_prime) then [range(0; $cpIndex) | 1 ]
else .product = 1
| reduce ($cpIndex | distinctPrimeFactors[]) as $prime (.;
substituteExponent(.polynomial; $prime) as $numerator
| .polynomial = exactDivision($numerator; .polynomial)
| .product *= $prime )
| substituteExponent(.polynomial; quo($cpIndex; .product) )
end;
 
# The Cyclotomic Polynomial equal to $dividend / $divisor
def exactDivision($dividend; $divisor):
reduce range(0; 1 + ($dividend|length) - ($divisor|length)) as $i ($dividend;
if .[$i] != 0
then reduce range(1; $divisor|length) as $j (.;
.[$i+$j] = .[$i+$j] - $divisor[$j] * .[$i] )
else .
end)
| .[0: 1 + length - ($divisor|length)];
 
### The tasks
def task1($n):
"Task 1: Cyclotomic polynomials for n <= \($n):",
( range(1;$n+1) | "CP[\(lpad(2))]: \(cycloPoly(.)|pp)" );
 
# For range(1;$n+1) as $c, report the first cpIndex which has a coefficient
# equal in magnitude to $c, possibly reporting others as well.
def task2($n):
def height: max(.[]|length); # i.e. abs
# update .summary and .todo
def register($cpIndex):
cycloPoly($cpIndex) as $poly
| if ($poly|height) < .todo[0] then .
else # it is a palindrome so we can halve the checks
reduce ($poly | .[0: quo(length + 1; 2)][]|length|select(.>1)) as $c (.;
if .summary[$c|tostring] then .
else .summary[$c|tostring] = $cpIndex
| .todo -= [$c]
| debug
end)
end;
 
{cpIndex:1, summary: {"1": 1}, todo: [range(2; $n + 1)]}
| until(.todo|length == 0;
if .cpIndex|is_prime then . else register(.cpIndex) end
| .cpIndex += 1)
| .summary
| (keys | sort_by(tonumber)[]) as $key
| "CP[\(.[$key]|lpad(5))] has a coefficient with magnitude \($key)"
;
 
task1(30),
"",
task2(10)
</syntaxhighlight>
{{output}}
<pre>
Task 1: Cyclotomic polynomials for n <= 30:
CP[ 1]: x-1
CP[ 2]: x+1
CP[ 3]: x²+x+1
CP[ 4]: x²+1
CP[ 5]: x⁴+x³+x²+x+1
CP[ 6]: x²-x+1
CP[ 7]: x⁶+x⁵+x⁴+x³+x²+x+1
CP[ 8]: x⁴+1
CP[ 9]: x⁶+x³+1
CP[10]: x⁴-x³+x²-x+1
CP[11]: x¹⁰+x⁹+x⁸+x⁷+x⁶+x⁵+x⁴+x³+x²+x+1
CP[12]: x⁴-x²+1
CP[13]: x¹²+x¹¹+x¹⁰+x⁹+x⁸+x⁷+x⁶+x⁵+x⁴+x³+x²+x+1
CP[14]: x⁶-x⁵+x⁴-x³+x²-x+1
CP[15]: x⁸-x⁷+x⁵-x⁴+x³-x+1
CP[16]: x⁸+1
CP[17]: x¹⁶+x¹⁵+x¹⁴+x¹³+x¹²+x¹¹+x¹⁰+x⁹+x⁸+x⁷+x⁶+x⁵+x⁴+x³+x²+x+1
CP[18]: x⁶-x³+1
CP[19]: x¹⁸+x¹⁷+x¹⁶+x¹⁵+x¹⁴+x¹³+x¹²+x¹¹+x¹⁰+x⁹+x⁸+x⁷+x⁶+x⁵+x⁴+x³+x²+x+1
CP[20]: x⁸-x⁶+x⁴-x²+1
CP[21]: x¹²-x¹¹+x⁹-x⁸+x⁶-x⁴+x³-x+1
CP[22]: x¹⁰-x⁹+x⁸-x⁷+x⁶-x⁵+x⁴-x³+x²-x+1
CP[23]: x²²+x²¹+x²⁰+x¹⁹+x¹⁸+x¹⁷+x¹⁶+x¹⁵+x¹⁴+x¹³+x¹²+x¹¹+x¹⁰+x⁹+x⁸+x⁷+x⁶+x⁵+x⁴+x³+x²+x+1
CP[24]: x⁸-x⁴+1
CP[25]: x²⁰+x¹⁵+x¹⁰+x⁵+1
CP[26]: x¹²-x¹¹+x¹⁰-x⁹+x⁸-x⁷+x⁶-x⁵+x⁴-x³+x²-x+1
CP[27]: x¹⁸+x⁹+1
CP[28]: x¹²-x¹⁰+x⁸-x⁶+x⁴-x²+1
CP[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
CP[30]: x⁸+x⁷-x⁵-x⁴-x³+x+1
 
Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:
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
CP[10465] has a coefficient with magnitude 11
CP[10465] has a coefficient with magnitude 12
CP[10465] has a coefficient with magnitude 13
CP[10465] has a coefficient with magnitude 14
</pre>
 
=={{header|Julia}}==
<langsyntaxhighlight lang="julia">using Primes, Polynomials
# memoize cache for recursive calls
Line 1,987 ⟶ 3,242:
println("The cyclotomic polynomial Φ(", abs(n), ") has a coefficient that is ", n < 0 ? -i : i)
end
</langsyntaxhighlight>{{out}}
<pre>
First 30 cyclotomic polynomials:
Line 2,034 ⟶ 3,289:
=={{header|Kotlin}}==
{{trans|Java}}
<langsyntaxhighlight lang="scala">import java.util.TreeMap
import kotlin.math.abs
import kotlin.math.pow
Line 2,490 ⟶ 3,745:
} else coefficient.toString() + "x^" + exponent
}
}</langsyntaxhighlight>
{{out}}
<pre>Task 1: cyclotomic polynomials for n <= 30:
Line 2,535 ⟶ 3,790:
CP[6545] has coefficient with magnitude = 9
CP[10465] has coefficient with magnitude = 10</pre>
 
=={{header|Maple}}==
<syntaxhighlight lang="maple">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]</syntaxhighlight>
 
=={{header|Mathematica}} / {{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">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++;
]</syntaxhighlight>
{{out}}
<pre>-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</pre>
 
=={{header|Nim}}==
{{trans|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.
 
<syntaxhighlight lang="nim">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</syntaxhighlight>
 
{{out}}
The program runs in 41 seconds on our reasonably performing laptop.
 
<pre>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</pre>
 
=={{header|PARI/GP}}==
Cyclotomic polynomials are a built-in function.
<syntaxhighlight lang="parigp">
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))
</syntaxhighlight>
 
{{out}}<pre>
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
</pre>
 
=={{header|Perl}}==
Conveniently, the module <code>Math::Polynomial::Cyclotomic</code> exists to do all the work. An <code>exponent too large</code> error prevents reaching the 10th step of the 2nd part of the task.
<langsyntaxhighlight lang="perl">use feature 'say';
use List::Util qw(first);
use Math::Polynomial::Cyclotomic qw(cyclo_poly_iterate);
Line 2,555 ⟶ 4,357:
$n++;
}
}</langsyntaxhighlight>
{{out}}
<pre>First 30 cyclotomic polynomials:
Line 2,603 ⟶ 4,405:
{{trans|Julia}}
Uses several routines from [[Polynomial_long_division#Phix]], tweaked slightly to check remainder is zero and trim the quotient.
<!--<syntaxhighlight lang="phix">(phixonline)-->
<lang Phix>-- demo\rosetta\Cyclotomic_Polynomial.exw
<span style="color: #000080;font-style:italic;">-- demo\rosetta\Cyclotomic_Polynomial.exw</span>
function degree(sequence p)
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
for i=length(p) to 1 by -1 do
<span style="color: #008080;">function</span> <span style="color: #000000;">degree</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
if p[i]!=0 then return i end if
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
end for
<span style="color: #008080;">if</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]!=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">i</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
return -1
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<span style="color: #008080;">return</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function poly_div(sequence n, d)
while length(d)<length(n) do d &=0 end while
<span style="color: #008080;">function</span> <span style="color: #000000;">poly_div</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
integer dn = degree(n),
<span style="color: #008080;">while</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)<</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> <span style="color: #000000;">d</span> <span style="color: #0000FF;">&=</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
dd = degree(d)
<span style="color: #004080;">integer</span> <span style="color: #000000;">dn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">degree</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span>
if dd<0 then throw("divide by zero") end if
<span style="color: #000000;">dd</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">degree</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
sequence quot = repeat(0,dn)
<span style="color: #008080;">if</span> <span style="color: #000000;">dd</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">throw</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"divide by zero"</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
while dn>=dd do
<span style="color: #004080;">sequence</span> <span style="color: #000000;">quot</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">dn</span><span style="color: #0000FF;">)</span>
integer k = dn-dd
<span style="color: #008080;">while</span> <span style="color: #000000;">dn</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">dd</span> <span style="color: #008080;">do</span>
integer qk = n[dn]/d[dd]
<span style="color: #004080;">integer</span> <span style="color: #000000;">k</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">dn</span><span style="color: #0000FF;">-</span><span style="color: #000000;">dd</span>
quot[k+1] = qk
<span style="color: #004080;">integer</span> <span style="color: #000000;">qk</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">[</span><span style="color: #000000;">dn</span><span style="color: #0000FF;">]/</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">dd</span><span style="color: #0000FF;">]</span>
sequence d2 = d[1..length(d)-k]
<span style="color: #000000;">quot</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">qk</span>
for i=1 to length(d2) do
<span style="color: #004080;">sequence</span> <span style="color: #000000;">d2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)-</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span>
n[-i] -= d2[-i]*qk
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d2</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
end for
<span style="color: #004080;">integer</span> <span style="color: #000000;">mi</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">i</span>
dn = degree(n)
<span style="color: #000000;">n</span><span style="color: #0000FF;">[</span><span style="color: #000000;">mi</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">d2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">mi</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">qk</span>
end while
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
-- return {quot,n} -- (n is now the remainder)
<span style="color: #000000;">dn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">degree</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
if n!=repeat(0,length(n)) then ?9/0 end if
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
while quot[$]=0 do quot = quot[1..$-1] end while
<span style="color: #000080;font-style:italic;">-- return {quot,n} -- (n is now the remainder)</span>
return quot
<span style="color: #008080;">if</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">!=</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
end function
<span style="color: #008080;">while</span> <span style="color: #000000;">quot</span><span style="color: #0000FF;">[$]=</span><span style="color: #000000;">0</span> <span style="color: #008080;">do</span> <span style="color: #000000;">quot</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">quot</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..$-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">quot</span>
function poly(sequence si)
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
-- display helper
string r = ""
<span style="color: #008080;">function</span> <span style="color: #000000;">poly</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">si</span><span style="color: #0000FF;">)</span>
for t=length(si) to 1 by -1 do
<span style="color: #000080;font-style:italic;">-- display helper</span>
integer sit = si[t]
<span style="color: #004080;">string</span> <span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">""</span>
if sit!=0 then
<span style="color: #008080;">for</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">si</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
if sit=1 and t>1 then
<span style="color: #004080;">integer</span> <span style="color: #000000;">sit</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">si</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">]</span>
r &= iff(r=""? "":" + ")
<span style="color: #008080;">if</span> <span style="color: #000000;">sit</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
elsif sit=-1 and t>1 then
<span style="color: #008080;">if</span> <span style="color: #000000;">sit</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">and</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">></span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span>
r &= iff(r=""?"-":" - ")
<span style="color: #000000;">r</span> <span style="color: #0000FF;">&=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #008000;">""</span><span style="color: #0000FF;">?</span> <span style="color: #008000;">""</span><span style="color: #0000FF;">:</span><span style="color: #008000;">" + "</span><span style="color: #0000FF;">)</span>
else
<span style="color: #008080;">elsif</span> <span style="color: #000000;">sit</span><span style="color: #0000FF;">=-</span><span style="color: #000000;">1</span> <span style="color: #008080;">and</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">></span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span>
if r!="" then
<span style="color: #000000;">r</span> <span style="color: #0000FF;">&=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #008000;">""</span><span style="color: #0000FF;">?</span><span style="color: #008000;">"-"</span><span style="color: #0000FF;">:</span><span style="color: #008000;">" - "</span><span style="color: #0000FF;">)</span>
r &= iff(sit<0?" - ":" + ")
<span sit style="color: abs(sit)#008080;">else</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">!=</span><span style="color: #008000;">""</span> <span style="color: #008080;">then</span>
end if
<span style="color: #000000;">r</span> <span style="color: #0000FF;">&=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sit</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span><span style="color: #0000FF;">?</span><span style="color: #008000;">" - "</span><span style="color: #0000FF;">:</span><span style="color: #008000;">" + "</span><span style="color: #0000FF;">)</span>
r &= sprintf("%d",sit)
<span style="color: #000000;">sit</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sit</span><span style="color: #0000FF;">)</span>
end if
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
r &= iff(t>1?"x"&iff(t>2?sprintf("^%d",t-1):""):"")
<span style="color: #000000;">r</span> <span style="color: #0000FF;">&=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"%d"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sit</span><span style="color: #0000FF;">)</span>
end if
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
end for
<span style="color: #000000;">r</span> <span style="color: #0000FF;">&=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">></span><span style="color: #000000;">1</span><span style="color: #0000FF;">?</span><span style="color: #008000;">"x"</span><span style="color: #0000FF;">&</span><span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">></span><span style="color: #000000;">2</span><span style="color: #0000FF;">?</span><span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"^%d"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">):</span><span style="color: #008000;">""</span><span style="color: #0000FF;">):</span><span style="color: #008000;">""</span><span style="color: #0000FF;">)</span>
if r="" then r="0" end if
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
return r
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<span style="color: #008080;">if</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #008000;">""</span> <span style="color: #008080;">then</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #008000;">"0"</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
--</Polynomial_long_division.exw>
<span style="color: #008080;">return</span> <span style="color: #000000;">r</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
--# memoize cache for recursive calls
<span style="color: #000080;font-style:italic;">--&lt;/Polynomial_long_division.exw&gt;
constant cyclotomics = new_dict({{1,{-1,1}},{2,{1,1}}})
--# memoize cache for recursive calls</span>
function cyclotomic(integer n)
<span style="color: #008080;">constant</span> <span style="color: #000000;">cyclotomics</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">new_dict</span><span style="color: #0000FF;">({{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,{-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}},{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}}})</span>
--
-- Calculate nth cyclotomic polynomial.
<span style="color: #008080;">function</span> <span style="color: #000000;">cyclotomic</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
-- See wikipedia article at bottom of section /wiki/Cyclotomic_polynomial#Fundamental_tools
<span style="color: #000080;font-style:italic;">--
-- The algorithm is reliable but slow for large n > 1000.
-- Calculate the nth cyclotomic polynomial.
--
-- See wikipedia article at bottom of section /wiki/Cyclotomic_polynomial#Fundamental_tools
sequence c
-- The algorithm is reliable but slow for large n &gt; 1000.
if getd_index(n,cyclotomics)!=NULL then
--</span>
c = getd(n,cyclotomics)
<span style="color: #004080;">sequence</span> <span style="color: #000000;">c</span>
else
<span style="color: #008080;">if</span> <span style="color: #7060A8;">getd_index</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cyclotomics</span><span style="color: #0000FF;">)!=</span><span style="color: #004600;">NULL</span> <span style="color: #008080;">then</span>
if is_prime(n) then
<span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">getd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cyclotomics</span><span style="color: #0000FF;">)</span>
c = repeat(1,n)
<span style="color: #008080;">else</span>
else -- recursive formula seen in wikipedia article
<span style="color: #008080;">if</span> <span style="color: #7060A8;">is_prime</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
c = -1&repeat(0,n-1)&1
<span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
sequence f = factors(n,-1)
<span style="color: #008080;">else</span> <span style="color: #000080;font-style:italic;">-- recursive formula seen in wikipedia article</span>
for i=1 to length(f) do
<span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">&</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)&</span><span style="color: #000000;">1</span>
c = poly_div(c,cyclotomic(f[i]))
<span style="color: #004080;">sequence</span> <span style="color: #000000;">f</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">factors</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
end for
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
end if
<span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">poly_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">c</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">cyclotomic</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])))</span>
setd(n,c,cyclotomics)
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end if
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
return c
<span style="color: #7060A8;">setd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cyclotomics</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">c</span>
for i=1 to 30 do
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
sequence z = cyclotomic(i)
string s = poly(z)
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">30</span> <span style="color: #008080;">do</span>
printf(1,"cp(%2d) = %s\n",{i,s})
<span style="color: #004080;">sequence</span> <span style="color: #000000;">z</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">cyclotomic</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span>
if i>1 and z!=reverse(z) then ?9/0 end if -- sanity check
<span style="color: #004080;">string</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">poly</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
end for
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"cp(%2d) = %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">></span><span style="color: #000000;">1</span> <span style="color: #008080;">and</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">!=</span><span style="color: #7060A8;">reverse</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #000080;font-style:italic;">-- sanity check</span>
integer found = 0, n = 1, cheat = 0
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
sequence fn = repeat(false,10),
nxt = {105,385,1365,1785,2805,3135,6545,6545,10465,10465}
<span style="color: #004080;">integer</span> <span style="color: #000000;">found</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">cheat</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
atom t1 = time()+1
<span style="color: #004080;">sequence</span> <span style="color: #000000;">fn</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #004600;">false</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">),</span>
puts(1,"\n")
<span style="color: #000000;">nxt</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">105</span><span style="color: #0000FF;">,</span><span style="color: #000000;">385</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1365</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1785</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2805</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3135</span><span style="color: #0000FF;">,</span><span style="color: #000000;">6545</span><span style="color: #0000FF;">,</span><span style="color: #000000;">6545</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10465</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10465</span><span style="color: #0000FF;">}</span>
while found<10 do
<span style="color: #004080;">atom</span> <span style="color: #000000;">t1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()+</span><span style="color: #000000;">1</span>
sequence z = cyclotomic(n)
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span>
for i=1 to length(z) do
<span style="color: #008080;">while</span> <span style="color: #000000;">found</span><span style="color: #0000FF;"><</span><span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">platform</span><span style="color: #0000FF;">()=</span><span style="color: #004600;">JS</span><span style="color: #0000FF;">?</span><span style="color: #000000;">5</span><span style="color: #0000FF;">:</span><span style="color: #000000;">10</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
atom azi = abs(z[i])
<span style="color: #004080;">sequence</span> <span style="color: #000000;">z</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">cyclotomic</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
if azi>=1 and azi<=10 and fn[azi]=0 then
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
printf(1,"cp(%d) has a coefficient with magnitude %d\n",{n,azi})
<span style="color: #004080;">atom</span> <span style="color: #000000;">azi</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
cheat = azi -- (comment this out to prevent cheating!)
<span style="color: #008080;">if</span> <span style="color: #000000;">azi</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">1</span> <span style="color: #008080;">and</span> <span style="color: #000000;">azi</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">10</span> <span style="color: #008080;">and</span> <span style="color: #000000;">fn</span><span style="color: #0000FF;">[</span><span style="color: #000000;">azi</span><span style="color: #0000FF;">]=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
found += 1
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"cp(%d) has a coefficient with magnitude %d\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">azi</span><span style="color: #0000FF;">})</span>
fn[azi] = true
<span style="color: #000000;">cheat</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">azi</span> <span style="color: #000080;font-style:italic;">-- (comment this out to prevent cheating!)</span>
t1 = time()+1
<span style="color: #000000;">found</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
end if
<span style="color: #000000;">fn</span><span style="color: #0000FF;">[</span><span style="color: #000000;">azi</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #004600;">true</span>
end for
<span style="color: #000000;">t1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()+</span><span style="color: #000000;">1</span>
if cheat then {n,cheat} = {nxt[cheat],0} else n += iff(n=1?4:10) end if
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
if time()>t1 then
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
printf(1,"working (%d) ...\r",n)
<span style="color: #008080;">if</span> <span style="color: #000000;">cheat</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cheat</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">nxt</span><span style="color: #0000FF;">[</span><span style="color: #000000;">cheat</span><span style="color: #0000FF;">],</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">else</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">+=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span><span style="color: #0000FF;">?</span><span style="color: #000000;">4</span><span style="color: #0000FF;">:</span><span style="color: #000000;">10</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
t1 = time()+1
<span style="color: #008080;">if</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()></span><span style="color: #000000;">t1</span> <span style="color: #008080;">and</span> <span style="color: #7060A8;">platform</span><span style="color: #0000FF;">()!=</span><span style="color: #004600;">JS</span> <span style="color: #008080;">then</span>
end if
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"working (%d) ...\r"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
end while</lang>
<span style="color: #000000;">t1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()+</span><span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<!--</syntaxhighlight>-->
{{out}}
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.<br>
Line 2,763 ⟶ 4,569:
 
=={{header|Python}}==
<langsyntaxhighlight lang="python">from itertools import count, chain
from collections import deque
 
Line 2,884 ⟶ 4,690:
while want in c or -want in c:
print(f'C[{want}]: {n}')
want += 1</langsyntaxhighlight>
{{out}}
Only showing first 10 polynomials to avoid clutter.
Line 2,930 ⟶ 4,736:
 
Uses the same library as Perl, so comes with the same caveats.
<syntaxhighlight lang="raku" perl6line>use Math::Polynomial::Cyclotomic:from<Perl5> <cyclo_poly_iterate cyclo_poly>;
 
say 'First 30 cyclotomic polynomials:';
Line 2,953 ⟶ 4,759:
sub super ($str) {
$str.subst( / '^' (\d+) /, { $0.trans([<0123456789>.comb] => [<⁰¹²³⁴⁵⁶⁷⁸⁹>.comb]) }, :g)
}</langsyntaxhighlight>
<pre>First 30 cyclotomic polynomials:
Φ(1) = (x - 1)
Line 2,998 ⟶ 4,804:
 
=={{header|Sidef}}==
Built-in:
Solution based on polynomial interpolation (slow).
<syntaxhighlight lang="ruby">say "First 30 cyclotomic polynomials:"
<lang ruby>var Poly = require('Math::Polynomial')
Poly.string_config(Hash(fold_sign => true, prefix => "", suffix => ""))
 
func poly_interpolation(v) {
v.len.of {|n| v.len.of {|k| n**k } }.msolve(v)
}
 
say "First 30 cyclotomic polynomials:"
for k in (1..30) {
var a =say ("Φ(#{k+1}).of {= ", cyclotomic(k, _) })
var Φ = poly_interpolation(a)
say ("Φ(#{k}) = ", Poly.new(Φ...))
}
 
say "\nSmallest cyclotomic polynomial with n or -n as a coefficient:"
for n in (1..10) { # very slow
var k = (1..Inf -> first {|k|
poly_interpolation(cyclotomic(k+1).ofcoeffs.any { cyclotomic(k, _) }).first { tail.abs == n }
})
say "Φ(#{k}) has coefficient with magnitude #{n}"
}</langsyntaxhighlight>
 
Slightly faster solution, using the '''Math::Polynomial::Cyclotomic''' Perl module.
<lang ruby>var Poly = require('Math::Polynomial')
require('Math::Polynomial::Cyclotomic')
 
Poly.string_config(Hash(fold_sign => true, prefix => "", suffix => ""))
 
say "First 30 cyclotomic polynomials:"
for k in (1..30) {
say ("Φ(#{k}) = ", Poly.new.cyclotomic(k))
}
 
say "\nSmallest cyclotomic polynomial with n or -n as a coefficient:"
for n in (1..10) {
var p = Poly.new
var k = (1..Inf -> first {|k|
[p.cyclotomic(k).coeff].first { .abs == n }
})
say "Φ(#{k}) has coefficient with magnitude = #{n}"
}</lang>
 
{{out}}
Line 3,084 ⟶ 4,861:
Φ(3135) has coefficient with magnitude = 7
^C
</pre>
 
=={{header|Visual Basic .NET}}==
{{trans|C++}}
<syntaxhighlight lang="vbnet">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</syntaxhighlight>
{{out}}
<pre>Task 1: cyclotomic polynomials for n <= 30:
CP[1] = x - 1
CP[2] = x + 1
CP[3] = x^2 + x + 1
CP[4] = x^2 + 1
CP[5] = x^4 + x^3 + x^2 + x + 1
CP[6] = x^2 - x + 1
CP[7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[8] = x^4 + 1
CP[9] = x^6 + x^3 + 1
CP[10] = x^4 - x^3 + x^2 - x + 1
CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[12] = x^4 - x^2 + 1
CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1
CP[16] = x^8 + 1
CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[18] = x^6 - x^3 + 1
CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[20] = x^8 - x^6 + x^4 - x^2 + 1
CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1
CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[24] = x^8 - x^4 + 1
CP[25] = x^20 + x^15 + x^10 + x^5 + 1
CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[27] = x^18 + x^9 + 1
CP[28] = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1
CP[29] = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[30] = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1
 
Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:
CP[1] has coefficient with magnitude = 1
CP[105] has coefficient with magnitude = 2
CP[385] has coefficient with magnitude = 3
CP[1365] has coefficient with magnitude = 4
CP[1785] has coefficient with magnitude = 5
CP[2805] has coefficient with magnitude = 6
CP[3135] has coefficient with magnitude = 7
CP[6545] has coefficient with magnitude = 8
CP[6545] has coefficient with magnitude = 9
CP[10465] has coefficient with magnitude = 10</pre>
 
=={{header|Wren}}==
===Version 1===
{{trans|Go}}
{{libheader|Wren-iterate}}
{{libheader|Wren-sort}}
{{libheader|Wren-math}}
{{libheader|Wren-fmt}}
Second part is very slow. Limited to first 7 to finish in a reasonable time - 5 minutes on my machine.
<syntaxhighlight lang="wren">import "./iterate" 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(coef, i))
}
} else if (factors.count == 1) {
var h = factors[2]
if (h) { // n == 2^h
cyclo.terms.addAll([Term.new(1, 1 << (h-1)), Term.new(1, 0)])
} else if (!factors[n]) {
// n == p ^ k
var p = 0
for (prime in factors.keys) p = prime
var k = factors[p]
for (i in 0...p) {
var pk = p.pow(k-1).floor
cyclo.terms.add(Term.new(1, i * pk))
}
}
} else if (factors.count == 2 && factors[2]) {
// n = 2^h * p^k
var p = 0
for (prime in factors.keys) if (prime != 2) p = prime
var coef = -1
var twoExp = 1 << (factors[2] - 1)
var k = factors[p]
for (i in 0...p) {
coef = coef * (-1)
var pk = p.pow(k-1).floor
cyclo.terms.add(Term.new(coef, i * twoExp * pk))
}
} else if (factors[2] && (n/2) % 2 == 1 && (n/2) > 1) {
// CP(2m)[x] == CP(-m)[x], n odd integer > 1
var cycloDiv2 = cycloPoly.call(n/2)
for (t in cycloDiv2.terms) {
var t2 = t
if (t.exp % 2 != 0) t2 = -t
cyclo.terms.add(t2)
}
} else if (algo == 0) {
// slow - uses basic definition
var divs = Int.properDivisors(n)
// polynomial: x^n - 1
var cyclo = Poly.new([1, n, -1, 0])
for (i in divs) {
var p = cycloPoly.call(i)
cyclo = cyclo / p
}
} else if (algo == 1) {
// faster - remove max divisor (and all divisors of max divisor)
// only one divide for all divisors of max divisor
var divs = Int.properDivisors(n)
var maxDiv = Nums.max(divs)
var divsExceptMax = divs.where { |d| maxDiv % d != 0 }.toList
// polynomial: ( x^n - 1 ) / ( x^m - 1 ), where m is the max divisor
cyclo = Poly.new([1, n, -1, 0])
cyclo = cyclo / Poly.new([1, maxDiv, -1, 0])
for (i in divsExceptMax) {
var p = cycloPoly.call(i)
cyclo = cyclo / 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]
var m = 1
cyclo = cycloPoly.call(m)
var primes = []
for (prime in factors.keys) primes.add(prime)
Sort.quick(primes)
for (prime in primes) {
// CP(m)[x]
var cycloM = cyclo
// compute CP(m)[x^p]
var terms = []
for (t in cycloM.terms) terms.add(Term.new(t.coef, t.exp * prime))
cyclo = Poly.new([])
cyclo.terms.addAll(terms)
cyclo.tidy()
cyclo = cyclo / cycloM
m = m * prime
}
// now, m is the largest square free divisor of n
var s = n / m
// Compute CP(n)[x] = CP(m)[x^s]
var terms = []
for (t in cyclo.terms) terms.add(Term.new(t.coef, t.exp * s))
cyclo = Poly.new([])
cyclo.terms.addAll(terms)
} else {
Fiber.abort("Invalid algorithm.")
}
cyclo.tidy()
computed[n] = cyclo
return cyclo
}
 
System.print("Task 1: cyclotomic polynomials for n <= 30:")
for (i in 1..30) {
var p = cycloPoly.call(i)
Fmt.print("CP[$2d] = $s", i, p)
}
 
System.print("\nTask 2: Smallest cyclotomic polynomial with n or -n as a coefficient:")
var n = 0
for (i in 1..7) {
while(true) {
n = n + 1
var cyclo = cycloPoly.call(n)
if (cyclo.hasCoefAbs(i)) {
Fmt.print("CP[$d] has coefficient with magnitude = $d", n, i)
n = n - 1
break
}
}
}</syntaxhighlight>
 
{{out}}
<pre>
Task 1: cyclotomic polynomials for n <= 30:
CP[ 1] = x - 1
CP[ 2] = x + 1
CP[ 3] = x^2 + x + 1
CP[ 4] = x^2 + 1
CP[ 5] = x^4 + x^3 + x^2 + x + 1
CP[ 6] = x^2 - x + 1
CP[ 7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[ 8] = x^4 + 1
CP[ 9] = x^6 + x^3 + 1
CP[10] = x^4 - x^3 + x^2 - x + 1
CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[12] = x^4 - x^2 + 1
CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1
CP[16] = x^8 + 1
CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[18] = x^6 - x^3 + 1
CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[20] = x^8 - x^6 + x^4 - x^2 + 1
CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1
CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[24] = x^8 - x^4 + 1
CP[25] = x^20 + x^15 + x^10 + x^5 + 1
CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
CP[27] = x^18 + x^9 + 1
CP[28] = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1
CP[29] = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
CP[30] = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1
 
Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:
CP[1] has coefficient with magnitude = 1
CP[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
</pre>
 
===Version 2===
{{trans|Java}}
A translation of the alternative version which completes the second part in about 33 seconds.
<syntaxhighlight lang="wren">import "./math" for Int
import "./fmt" for Fmt
 
class CP {
// Return the Cyclotomic Polynomial of order 'cpIndex' as a list of coefficients,
// where, for example, the polynomial 3x^2 - 1 is represented by the list [3, 0, -1].
static cycloPoly(cpIndex) {
var polynomial = [1, -1]
if (cpIndex == 1) return polynomial
if (Int.isPrime(cpIndex)) return List.filled(cpIndex, 1)
var primes = Int.distinctPrimeFactors(cpIndex)
var product = 1
for (prime in primes) {
var numerator = substituteExponent(polynomial, prime)
polynomial = exactDivision(numerator, polynomial)
product = product * prime
}
return substituteExponent(polynomial, Int.quo(cpIndex, product))
}
 
// Return the Cyclotomic Polynomial obtained from 'polynomial'
// by replacing x with x^'exponent'.
static substituteExponent(polynomial, exponent) {
var result = List.filled(exponent * (polynomial.count - 1) + 1, 0)
for (i in polynomial.count-1..0) result[i*exponent] = polynomial[i]
return result
}
 
// Return the Cyclotomic Polynomial equal to 'dividend' / 'divisor'.
// The division is always exact.
static exactDivision(dividend, divisor) {
var result = dividend.toList
for (i in 0..dividend.count - divisor.count) {
if (result[i] != 0) {
for (j in 1...divisor.count) {
result[i+j] = result[i+j] - divisor[j] * result[i]
}
}
}
return result[0..result.count - divisor.count]
}
 
// Return whether 'polynomial' has a coefficient of equal magnitude
// to 'coefficient'.
static hasHeight(polynomial, coefficient) {
for (i in 0..Int.quo(polynomial.count + 1, 2)) {
if (polynomial[i].abs == coefficient) return true
}
return false
}
}
 
System.print("Task 1: Cyclotomic polynomials for n <= 30:")
for (cpIndex in 1..30) {
Fmt.write("CP[$2d] = ", cpIndex)
Fmt.pprint("$d", CP.cycloPoly(cpIndex), "", "x")
}
 
System.print("\nTask 2: Smallest cyclotomic polynomial with n or -n as a coefficient:")
System.print("CP[ 1] has a coefficient with magnitude 1")
var cpIndex = 2
for (coeff in 2..10) {
while (Int.isPrime(cpIndex) || !CP.hasHeight(CP.cycloPoly(cpIndex), coeff)) {
cpIndex = cpIndex + 1
}
Fmt.print("CP[$5d] has a coefficient with magnitude $d", cpIndex, coeff)
}</syntaxhighlight>
 
{{out}}
<pre>
Task 1: Cyclotomic polynomials for n <= 30:
CP[ 1] = x - 1
CP[ 2] = x + 1
CP[ 3] = x² + x + 1
CP[ 4] = x² + 1
CP[ 5] = x⁴ + x³ + x² + x + 1
CP[ 6] = x² - x + 1
CP[ 7] = x⁶ + x⁵ + x⁴ + x³ + x² + x + 1
CP[ 8] = x⁴ + 1
CP[ 9] = x⁶ + x³ + 1
CP[10] = x⁴ - x³ + x² - x + 1
CP[11] = x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1
CP[12] = x⁴ - x² + 1
CP[13] = x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1
CP[14] = x⁶ - x⁵ + x⁴ - x³ + x² - x + 1
CP[15] = x⁸ - x⁷ + x⁵ - x⁴ + x³ - x + 1
CP[16] = x⁸ + 1
CP[17] = x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1
CP[18] = x⁶ - x³ + 1
CP[19] = x¹⁸ + x¹⁷ + x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1
CP[20] = x⁸ - x⁶ + x⁴ - x² + 1
CP[21] = x¹² - x¹¹ + x⁹ - x⁸ + x⁶ - x⁴ + x³ - x + 1
CP[22] = x¹⁰ - x⁹ + x⁸ - x⁷ + x⁶ - x⁵ + x⁴ - x³ + x² - x + 1
CP[23] = x²² + x²¹ + x²⁰ + x¹⁹ + x¹⁸ + x¹⁷ + x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1
CP[24] = x⁸ - x⁴ + 1
CP[25] = x²⁰ + x¹⁵ + x¹⁰ + x⁵ + 1
CP[26] = x¹² - x¹¹ + x¹⁰ - x⁹ + x⁸ - x⁷ + x⁶ - x⁵ + x⁴ - x³ + x² - x + 1
CP[27] = x¹⁸ + x⁹ + 1
CP[28] = x¹² - x¹⁰ + x⁸ - x⁶ + x⁴ - x² + 1
CP[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
CP[30] = x⁸ + x⁷ - x⁵ - x⁴ - x³ + x + 1
 
Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:
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
</pre>
2,472

edits