Cyclotomic polynomial
The nth Cyclotomic polynomial, for any positive integer n, is the unique irreducible polynomial of largest degree with integer coefficients that is a divisor of x^n − 1, and is not a divisor of x^k − 1 for any k < n.
You are encouraged to solve this task according to the task description, using any language you may know.
- Task
- Find and print the first 30 cyclotomic polynomials.
- Find and print the order of the first 10 cyclotomic polynomials that have n or -n as a coefficient.
- See also
- Wikipedia article, Cyclotomic polynomial, showing ways to calculate them.
- The sequence A013594 with the smallest order of cyclotomic polynomial containing n or -n as a coefficient.
Go
<lang go>package main
import (
"fmt" "log" "math" "sort" "strings"
)
const (
algo = 2 maxAllFactors = 100000
)
func iabs(i int) int {
if i < 0 {
return -i
}
return i
}
type term struct{ coef, exp int }
func (t term) mul(t2 term) term {
return term{t.coef * t2.coef, t.exp + t2.exp}
}
func (t term) add(t2 term) term {
if t.exp != t2.exp {
log.Fatal("exponents unequal in term.add method")
}
return term{t.coef + t2.coef, t.exp}
}
func (t term) negate() term { return term{-t.coef, t.exp} }
func (t term) String() string {
switch {
case t.coef == 0:
return "0"
case t.exp == 0:
return fmt.Sprintf("%d", t.coef)
case t.coef == 1:
if t.exp == 1 {
return "x"
} else {
return fmt.Sprintf("x^%d", t.exp)
}
case t.exp == 1:
return fmt.Sprintf("%dx", t.coef)
}
return fmt.Sprintf("%dx^%d", t.coef, t.exp)
}
type poly struct{ terms []term }
// pass coef, exp in pairs as parameters func newPoly(values ...int) poly {
le := len(values)
if le == 0 {
return poly{[]term{term{0, 0}}}
}
if le%2 != 0 {
log.Fatalf("odd number of parameters (%d) passed to newPoly function", le)
}
var terms []term
for i := 0; i < le; i += 2 {
terms = append(terms, term{values[i], values[i+1]})
}
p := poly{terms}.tidy()
return p
}
func (p poly) hasCoefAbs(coef int) bool {
for _, t := range p.terms {
if iabs(t.coef) == coef {
return true
}
}
return false
}
func (p poly) add(p2 poly) poly {
p3 := newPoly()
le, le2 := len(p.terms), len(p2.terms)
for le > 0 || le2 > 0 {
if le == 0 {
p3.terms = append(p3.terms, p2.terms[le2-1])
le2--
} else if le2 == 0 {
p3.terms = append(p3.terms, p.terms[le-1])
le--
} else {
t := p.terms[le-1]
t2 := p2.terms[le2-1]
if t.exp == t2.exp {
t3 := t.add(t2)
if t3.coef != 0 {
p3.terms = append(p3.terms, t3)
}
le--
le2--
} else if t.exp < t2.exp {
p3.terms = append(p3.terms, t)
le--
} else {
p3.terms = append(p3.terms, t2)
le2--
}
}
}
return p3.tidy()
}
func (p poly) addTerm(t term) poly {
q := newPoly()
added := false
for i := 0; i < len(p.terms); i++ {
ct := p.terms[i]
if ct.exp == t.exp {
added = true
if ct.coef+t.coef != 0 {
q.terms = append(q.terms, ct.add(t))
}
} else {
q.terms = append(q.terms, ct)
}
}
if !added {
q.terms = append(q.terms, t)
}
return q.tidy()
}
func (p poly) mulTerm(t term) poly {
q := newPoly()
for i := 0; i < len(p.terms); i++ {
ct := p.terms[i]
q.terms = append(q.terms, ct.mul(t))
}
return q.tidy()
}
func (p poly) div(v poly) poly {
q := newPoly()
lcv := v.leadingCoef()
dv := v.degree()
for p.degree() >= v.degree() {
lcp := p.leadingCoef()
s := lcp / lcv
t := term{s, p.degree() - dv}
q = q.addTerm(t)
p = p.add(v.mulTerm(t.negate()))
}
return q.tidy()
}
func (p poly) leadingCoef() int {
return p.terms[0].coef
}
func (p poly) degree() int {
return p.terms[0].exp
}
func (p poly) String() string {
var sb strings.Builder
first := true
for _, t := range p.terms {
if first {
sb.WriteString(t.String())
first = false
} else {
sb.WriteString(" ")
if t.coef > 0 {
sb.WriteString("+ ")
sb.WriteString(t.String())
} else {
sb.WriteString("- ")
sb.WriteString(t.negate().String())
}
}
}
return sb.String()
}
// in place descending sort by term.exp func (p poly) sortTerms() {
sort.Slice(p.terms, func(i, j int) bool {
return p.terms[i].exp > p.terms[j].exp
})
}
// sort terms and remove any unnecesary zero terms func (p poly) tidy() poly {
p.sortTerms()
if p.degree() == 0 {
return p
}
for i := len(p.terms) - 1; i >= 0; i-- {
if p.terms[i].coef == 0 {
copy(p.terms[i:], p.terms[i+1:])
p.terms[len(p.terms)-1] = term{0, 0}
p.terms = p.terms[:len(p.terms)-1]
}
}
if len(p.terms) == 0 {
p.terms = append(p.terms, term{0, 0})
}
return p
}
func getDivisors(n int) []int {
var divs []int
sqrt := int(math.Sqrt(float64(n)))
for i := 1; i <= sqrt; i++ {
if n%i == 0 {
divs = append(divs, i)
d := n / i
if d != i && d != n {
divs = append(divs, d)
}
}
}
return divs
}
var (
computed = make(map[int]poly) allFactors = make(map[int]map[int]int)
)
func init() {
f := map[int]int{2: 1}
allFactors[2] = f
}
func getFactors(n int) map[int]int {
if f, ok := allFactors[n]; ok {
return f
}
factors := make(map[int]int)
if n%2 == 0 {
factorsDivTwo := getFactors(n / 2)
for k, v := range factorsDivTwo {
factors[k] = v
}
factors[2]++
if n < maxAllFactors {
allFactors[n] = factors
}
return factors
}
prime := true
sqrt := int(math.Sqrt(float64(n)))
for i := 3; i <= sqrt; i += 2 {
if n%i == 0 {
prime = false
for k, v := range getFactors(n / i) {
factors[k] = v
}
factors[i]++
if n < maxAllFactors {
allFactors[n] = factors
}
return factors
}
}
if prime {
factors[n] = 1
if n < maxAllFactors {
allFactors[n] = factors
}
}
return factors
}
func cycloPoly(n int) poly {
if p, ok := computed[n]; ok {
return p
}
if n == 1 {
// polynomial: x - 1
p := newPoly(1, 1, -1, 0)
computed[1] = p
return p
}
factors := getFactors(n)
cyclo := newPoly()
if _, ok := factors[n]; ok {
// n is prime
for i := 0; i < n; i++ {
cyclo.terms = append(cyclo.terms, term{1, i})
}
} else if len(factors) == 2 && factors[2] == 1 && factors[n/2] == 1 {
// n == 2p
prime := n / 2
coef := -1
for i := 0; i < prime; i++ {
coef *= -1
cyclo.terms = append(cyclo.terms, term{coef, i})
}
} else if len(factors) == 1 {
if h, ok := factors[2]; ok {
// n == 2^h
cyclo.terms = append(cyclo.terms, term{1, 1 << (h - 1)}, term{1, 0})
} else if _, ok := factors[n]; !ok {
// n == p ^ k
p := 0
for prime := range factors {
p = prime
}
k := factors[p]
for i := 0; i < p; i++ {
pk := int(math.Pow(float64(p), float64(k-1)))
cyclo.terms = append(cyclo.terms, term{1, i * pk})
}
}
} else if len(factors) == 2 && factors[2] != 0 {
// n = 2^h * p^k
p := 0
for prime := range factors {
if prime != 2 {
p = prime
}
}
coef := -1
twoExp := 1 << (factors[2] - 1)
k := factors[p]
for i := 0; i < p; i++ {
coef *= -1
pk := int(math.Pow(float64(p), float64(k-1)))
cyclo.terms = append(cyclo.terms, term{coef, i * twoExp * pk})
}
} else if factors[2] != 0 && ((n/2)%2 == 1) && (n/2) > 1 {
// CP(2m)[x] == CP(-m)[x], n odd integer > 1
cycloDiv2 := cycloPoly(n / 2)
for _, t := range cycloDiv2.terms {
t2 := t
if t.exp%2 != 0 {
t2 = t.negate()
}
cyclo.terms = append(cyclo.terms, t2)
}
} else if algo == 0 {
// slow - uses basic definition
divs := getDivisors(n)
// polynomial: x^n - 1
cyclo = newPoly(1, n, -1, 0)
for _, i := range divs {
p := cycloPoly(i)
cyclo = cyclo.div(p)
}
} else if algo == 1 {
// faster - remove max divisor (and all divisors of max divisor)
// only one divide for all divisors of max divisor
divs := getDivisors(n)
maxDiv := math.MinInt32
for _, d := range divs {
if d > maxDiv {
maxDiv = d
}
}
var divsExceptMax []int
for _, d := range divs {
if maxDiv%d != 0 {
divsExceptMax = append(divsExceptMax, d)
}
}
// polynomial: ( x^n - 1 ) / ( x^m - 1 ), where m is the max divisor
cyclo = newPoly(1, n, -1, 0)
cyclo = cyclo.div(newPoly(1, maxDiv, -1, 0))
for _, i := range divsExceptMax {
p := cycloPoly(i)
cyclo = cyclo.div(p)
}
} else if algo == 2 {
// fastest
// let p, q be primes such that p does not divide n, and q divides n
// then CP(np)[x] = CP(n)[x^p] / CP(n)[x]
m := 1
cyclo = cycloPoly(m)
var primes []int
for prime := range factors {
primes = append(primes, prime)
}
sort.Ints(primes)
for _, prime := range primes {
// CP(m)[x]
cycloM := cyclo
// compute CP(m)[x^p]
var terms []term
for _, t := range cycloM.terms {
terms = append(terms, term{t.coef, t.exp * prime})
}
cyclo = newPoly()
cyclo.terms = append(cyclo.terms, terms...)
cyclo = cyclo.tidy()
cyclo = cyclo.div(cycloM)
m *= prime
}
// now, m is the largest square free divisor of n
s := n / m
// Compute CP(n)[x] = CP(m)[x^s]
var terms []term
for _, t := range cyclo.terms {
terms = append(terms, term{t.coef, t.exp * s})
}
cyclo = newPoly()
cyclo.terms = append(cyclo.terms, terms...)
} else {
log.Fatal("invalid algorithm")
}
cyclo = cyclo.tidy()
computed[n] = cyclo
return cyclo
}
func main() {
fmt.Println("Task 1: cyclotomic polynomials for n <= 30:")
for i := 1; i <= 30; i++ {
p := cycloPoly(i)
fmt.Printf("CP[%2d] = %s\n", i, p)
}
fmt.Println("\nTask 2: Smallest cyclotomic polynomial with n or -n as a coefficient:")
n := 0
for i := 1; i <= 10; i++ {
for {
n++
cyclo := cycloPoly(n)
if cyclo.hasCoefAbs(i) {
fmt.Printf("CP[%d] has coefficient with magnitude = %d\n", n, i)
n--
break
}
}
}
}</lang>
- Output:
Task 1: cyclotomic polynomials for n <= 30: CP[ 1] = x - 1 CP[ 2] = x + 1 CP[ 3] = x^2 + x + 1 CP[ 4] = x^2 + 1 CP[ 5] = x^4 + x^3 + x^2 + x + 1 CP[ 6] = x^2 - x + 1 CP[ 7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[ 8] = x^4 + 1 CP[ 9] = x^6 + x^3 + 1 CP[10] = x^4 - x^3 + x^2 - x + 1 CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[12] = x^4 - x^2 + 1 CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 CP[16] = x^8 + 1 CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[18] = x^6 - x^3 + 1 CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[20] = x^8 - x^6 + x^4 - x^2 + 1 CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[24] = x^8 - x^4 + 1 CP[25] = x^20 + x^15 + x^10 + x^5 + 1 CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[27] = x^18 + x^9 + 1 CP[28] = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 CP[29] = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[30] = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient: CP[1] has coefficient with magnitude = 1 CP[105] has coefficient with magnitude = 2 CP[385] has coefficient with magnitude = 3 CP[1365] has coefficient with magnitude = 4 CP[1785] has coefficient with magnitude = 5 CP[2805] has coefficient with magnitude = 6 CP[3135] has coefficient with magnitude = 7 CP[6545] has coefficient with magnitude = 8 CP[6545] has coefficient with magnitude = 9 CP[10465] has coefficient with magnitude = 10
Java
<lang java> import java.util.ArrayList; import java.util.Collections; import java.util.Comparator; import java.util.HashMap; import java.util.List; import java.util.Map; import java.util.TreeMap;
public class CyclotomicPolynomial {
@SuppressWarnings("unused")
private static int divisions = 0;
private static int algorithm = 2;
public static void main(String[] args) throws Exception {
System.out.println("Task 1: cyclotomic polynomials for n <= 30:");
for ( int i = 1 ; i <= 30 ; i++ ) {
Polynomial p = cyclotomicPolynomial(i);
System.out.printf("CP[%d] = %s%n", i, p);
}
System.out.println("Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:");
int n = 0;
for ( int i = 1 ; i <= 10 ; i++ ) {
while ( true ) {
n++;
Polynomial cyclo = cyclotomicPolynomial(n);
if ( cyclo.hasCoefficientAbs(i) ) {
System.out.printf("CP[%d] has coefficient with magnitude = %d%n", n, i);
n--;
break;
}
}
}
}
private static final Map<Integer, Polynomial> COMPUTED = new HashMap<>();
private static Polynomial cyclotomicPolynomial(int n) {
if ( COMPUTED.containsKey(n) ) {
return COMPUTED.get(n);
}
//System.out.println("COMPUTE: n = " + n);
if ( n == 1 ) {
// Polynomial: x - 1
Polynomial p = new Polynomial(1, 1, -1, 0);
COMPUTED.put(1, p);
return p;
}
Map<Integer,Integer> factors = getFactors(n);
if ( factors.containsKey(n) ) {
// n prime
List<Term> termList = new ArrayList<>();
for ( int index = 0 ; index < n ; index++ ) {
termList.add(new Term(1, index));
}
Polynomial cyclo = new Polynomial(termList);
COMPUTED.put(n, cyclo);
return cyclo;
}
else if ( factors.size() == 2 && factors.containsKey(2) && factors.get(2) == 1 && factors.containsKey(n/2) && factors.get(n/2) == 1 ) {
// n = 2p
int prime = n/2;
List<Term> termList = new ArrayList<>();
int coeff = -1;
for ( int index = 0 ; index < prime ; index++ ) {
coeff *= -1;
termList.add(new Term(coeff, index));
}
Polynomial cyclo = new Polynomial(termList);
COMPUTED.put(n, cyclo);
return cyclo;
}
else if ( factors.size() == 1 && factors.containsKey(2) ) {
// n = 2^h
int h = factors.get(2);
List<Term> termList = new ArrayList<>();
termList.add(new Term(1, (int) Math.pow(2, h-1)));
termList.add(new Term(1, 0));
Polynomial cyclo = new Polynomial(termList);
COMPUTED.put(n, cyclo);
return cyclo;
}
else if ( factors.size() == 1 && ! factors.containsKey(n) ) {
// n = p^k
int p = 0;
for ( int prime : factors.keySet() ) {
p = prime;
}
int k = factors.get(p);
List<Term> termList = new ArrayList<>();
for ( int index = 0 ; index < p ; index++ ) {
termList.add(new Term(1, index * (int) Math.pow(p, k-1)));
}
Polynomial cyclo = new Polynomial(termList);
COMPUTED.put(n, cyclo);
return cyclo;
}
else if ( factors.size() == 2 && factors.containsKey(2) ) {
// n = 2^h * p^k
int p = 0;
for ( int prime : factors.keySet() ) {
if ( prime != 2 ) {
p = prime;
}
}
List<Term> termList = new ArrayList<>();
int coeff = -1;
int twoExp = (int) Math.pow(2, factors.get(2)-1);
int k = factors.get(p);
for ( int index = 0 ; index < p ; index++ ) {
coeff *= -1;
termList.add(new Term(coeff, index * twoExp * (int) Math.pow(p, k-1)));
}
Polynomial cyclo = new Polynomial(termList);
COMPUTED.put(n, cyclo);
return cyclo;
}
else if ( factors.containsKey(2) && ((n/2) % 2 == 1) && (n/2) > 1 ) {
// CP(2m)[x] = CP(-m)[x], n odd integer > 1
Polynomial cycloDiv2 = cyclotomicPolynomial(n/2);
List<Term> termList = new ArrayList<>();
for ( Term term : cycloDiv2.polynomialTerms ) {
termList.add(term.exponent % 2 == 0 ? term : term.negate());
}
Polynomial cyclo = new Polynomial(termList);
COMPUTED.put(n, cyclo);
return cyclo;
}
// General Case
if ( algorithm == 0 ) {
// Slow - uses basic definition.
List<Integer> divisors = getDivisors(n);
// Polynomial: ( x^n - 1 )
Polynomial cyclo = new Polynomial(1, n, -1, 0);
for ( int i : divisors ) {
Polynomial p = cyclotomicPolynomial(i);
cyclo = cyclo.divide(p);
}
COMPUTED.put(n, cyclo);
return cyclo;
}
else if ( algorithm == 1 ) {
// Faster. Remove Max divisor (and all divisors of max divisor) - only one divide for all divisors of Max Divisor
List<Integer> divisors = getDivisors(n);
int maxDivisor = Integer.MIN_VALUE;
for ( int div : divisors ) {
maxDivisor = Math.max(maxDivisor, div);
}
List<Integer> divisorsExceptMax = new ArrayList<Integer>();
for ( int div : divisors ) {
if ( maxDivisor % div != 0 ) {
divisorsExceptMax.add(div);
}
}
// Polynomial: ( x^n - 1 ) / ( x^m - 1 ), where m is the max divisor
Polynomial cyclo = new Polynomial(1, n, -1, 0).divide(new Polynomial(1, maxDivisor, -1, 0));
for ( int i : divisorsExceptMax ) {
Polynomial p = cyclotomicPolynomial(i);
cyclo = cyclo.divide(p);
}
COMPUTED.put(n, cyclo);
return cyclo;
}
else if ( algorithm == 2 ) {
// Fastest
// Let p ; q be primes such that p does not divide n, and q q divides n.
// Then CP(np)[x] = CP(n)[x^p] / CP(n)[x]
int m = 1;
Polynomial cyclo = cyclotomicPolynomial(m);
List<Integer> primes = new ArrayList<>(factors.keySet());
Collections.sort(primes);
for ( int prime : primes ) {
// CP(m)[x]
Polynomial cycloM = cyclo;
// Compute CP(m)[x^p].
List<Term> termList = new ArrayList<>();
for ( Term t : cycloM.polynomialTerms ) {
termList.add(new Term(t.coefficient, t.exponent * prime));
}
cyclo = new Polynomial(termList).divide(cycloM);
m = m * prime;
}
// Now, m is the largest square free divisor of n
int s = n / m;
// Compute CP(n)[x] = CP(m)[x^s]
List<Term> termList = new ArrayList<>();
for ( Term t : cyclo.polynomialTerms ) {
termList.add(new Term(t.coefficient, t.exponent * s));
}
cyclo = new Polynomial(termList);
COMPUTED.put(n, cyclo);
return cyclo;
}
else {
throw new RuntimeException("ERROR 103: Invalid algorithm.");
}
}
private static final List<Integer> getDivisors(int number) {
List<Integer> divisors = new ArrayList<Integer>();
long sqrt = (long) Math.sqrt(number);
for ( int i = 1 ; i <= sqrt ; i++ ) {
if ( number % i == 0 ) {
divisors.add(i);
int div = number / i;
if ( div != i && div != number ) {
divisors.add(div);
}
}
}
return divisors;
}
private static final Map<Integer,Map<Integer,Integer>> allFactors = new TreeMap<Integer,Map<Integer,Integer>>();
static {
Map<Integer,Integer> factors = new TreeMap<Integer,Integer>();
factors.put(2, 1);
allFactors.put(2, factors);
}
public static Integer MAX_ALL_FACTORS = 100000;
public static final Map<Integer,Integer> getFactors(Integer number) {
if ( allFactors.containsKey(number) ) {
return allFactors.get(number);
}
Map<Integer,Integer> factors = new TreeMap<Integer,Integer>();
if ( number % 2 == 0 ) {
Map<Integer,Integer> factorsdDivTwo = getFactors(number/2);
factors.putAll(factorsdDivTwo);
factors.merge(2, 1, (v1, v2) -> v1 + v2);
if ( number < MAX_ALL_FACTORS )
allFactors.put(number, factors);
return factors;
}
boolean prime = true;
long sqrt = (long) Math.sqrt(number);
for ( int i = 3 ; i <= sqrt ; i += 2 ) {
if ( number % i == 0 ) {
prime = false;
factors.putAll(getFactors(number/i));
factors.merge(i, 1, (v1, v2) -> v1 + v2);
if ( number < MAX_ALL_FACTORS )
allFactors.put(number, factors);
return factors;
}
}
if ( prime ) {
factors.put(number, 1);
if ( number < MAX_ALL_FACTORS )
allFactors.put(number, factors);
}
return factors;
}
private static final class Polynomial {
private List<Term> polynomialTerms;
// Format - coeff, exp, coeff, exp, (repeating in pairs) . . .
public Polynomial(int ... values) {
if ( values.length % 2 != 0 ) {
throw new IllegalArgumentException("ERROR 102: Polynomial constructor. Length must be even. Length = " + values.length);
}
polynomialTerms = new ArrayList<>();
for ( int i = 0 ; i < values.length ; i += 2 ) {
Term t = new Term(values[i], values[i+1]);
polynomialTerms.add(t);
}
Collections.sort(polynomialTerms, new TermSorter());
}
public Polynomial() {
// zero
polynomialTerms = new ArrayList<>();
polynomialTerms.add(new Term(0,0));
}
private boolean hasCoefficientAbs(int coeff) {
for ( Term term : polynomialTerms ) {
if ( Math.abs(term.coefficient) == coeff ) {
return true;
}
}
return false;
}
private Polynomial(List<Term> termList) {
if ( termList.size() == 0 ) {
// zero
termList.add(new Term(0,0));
}
else {
// Remove zero terms if needed
for ( int i = 0 ; i < termList.size() ; i++ ) {
if ( termList.get(i).coefficient == 0 ) {
termList.remove(i);
}
}
}
if ( termList.size() == 0 ) {
// zero
termList.add(new Term(0,0));
}
polynomialTerms = termList;
Collections.sort(polynomialTerms, new TermSorter());
}
public Polynomial divide(Polynomial v) {
//long start = System.currentTimeMillis();
divisions++;
Polynomial q = new Polynomial();
Polynomial r = this;
long lcv = v.leadingCoefficient();
long dv = v.degree();
while ( r.degree() >= v.degree() ) {
long lcr = r.leadingCoefficient();
long s = lcr / lcv; // Integer division
Term term = new Term(s, r.degree() - dv);
q = q.add(term);
r = r.add(v.multiply(term.negate()));
}
//long end = System.currentTimeMillis();
//System.out.printf("Divide: Elapsed = %d, Degree 1 = %d, Degree 2 = %d%n", (end-start), this.degree(), v.degree());
return q;
}
public Polynomial add(Polynomial polynomial) {
List<Term> termList = new ArrayList<>();
int thisCount = polynomialTerms.size();
int polyCount = polynomial.polynomialTerms.size();
while ( thisCount > 0 || polyCount > 0 ) {
Term thisTerm = thisCount == 0 ? null : polynomialTerms.get(thisCount-1);
Term polyTerm = polyCount == 0 ? null : polynomial.polynomialTerms.get(polyCount-1);
if ( thisTerm == null ) {
termList.add(polyTerm.clone());
polyCount--;
}
else if (polyTerm == null ) {
termList.add(thisTerm.clone());
thisCount--;
}
else if ( thisTerm.degree() == polyTerm.degree() ) {
Term t = thisTerm.add(polyTerm);
if ( t.coefficient != 0 ) {
termList.add(t);
}
thisCount--;
polyCount--;
}
else if ( thisTerm.degree() < polyTerm.degree() ) {
termList.add(thisTerm.clone());
thisCount--;
}
else {
termList.add(polyTerm.clone());
polyCount--;
}
}
return new Polynomial(termList);
}
public Polynomial add(Term term) {
List<Term> termList = new ArrayList<>();
boolean added = false;
for ( int index = 0 ; index < polynomialTerms.size() ; index++ ) {
Term currentTerm = polynomialTerms.get(index);
if ( currentTerm.exponent == term.exponent ) {
added = true;
if ( currentTerm.coefficient + term.coefficient != 0 ) {
termList.add(currentTerm.add(term));
}
}
else {
termList.add(currentTerm.clone());
}
}
if ( ! added ) {
termList.add(term.clone());
}
return new Polynomial(termList);
}
public Polynomial multiply(Term term) {
List<Term> termList = new ArrayList<>();
for ( int index = 0 ; index < polynomialTerms.size() ; index++ ) {
Term currentTerm = polynomialTerms.get(index);
termList.add(currentTerm.clone().multiply(term));
}
return new Polynomial(termList);
}
public Polynomial clone() {
List<Term> clone = new ArrayList<>();
for ( Term t : polynomialTerms ) {
clone.add(new Term(t.coefficient, t.exponent));
}
return new Polynomial(clone);
}
public long leadingCoefficient() {
// long coefficient = 0; // long degree = Integer.MIN_VALUE; // for ( Term t : polynomialTerms ) { // if ( t.degree() > degree ) { // coefficient = t.coefficient; // degree = t.degree(); // } // }
return polynomialTerms.get(0).coefficient;
}
public long degree() {
// long degree = Integer.MIN_VALUE; // for ( Term t : polynomialTerms ) { // if ( t.degree() > degree ) { // degree = t.degree(); // } // }
return polynomialTerms.get(0).exponent;
}
@Override
public String toString() {
StringBuilder sb = new StringBuilder();
boolean first = true;
for ( Term term : polynomialTerms ) {
if ( first ) {
sb.append(term);
first = false;
}
else {
sb.append(" ");
if ( term.coefficient > 0 ) {
sb.append("+ ");
sb.append(term);
}
else {
sb.append("- ");
sb.append(term.negate());
}
}
}
return sb.toString();
}
}
private static final class TermSorter implements Comparator<Term> {
@Override
public int compare(Term o1, Term o2) {
return (int) (o2.exponent - o1.exponent);
}
}
// Note: Cyclotomic Polynomials have small coefficients. Not appropriate for general polynomial usage.
private static final class Term {
long coefficient;
long exponent;
public Term(long c, long e) {
coefficient = c;
exponent = e;
}
public Term clone() {
return new Term(coefficient, exponent);
}
public Term multiply(Term term) {
return new Term(coefficient * term.coefficient, exponent + term.exponent);
}
public Term add(Term term) {
if ( exponent != term.exponent ) {
throw new RuntimeException("ERROR 102: Exponents not equal.");
}
return new Term(coefficient + term.coefficient, exponent);
}
public Term negate() {
return new Term(-coefficient, exponent);
}
public long degree() {
return exponent;
}
@Override
public String toString() {
if ( coefficient == 0 ) {
return "0";
}
if ( exponent == 0 ) {
return "" + coefficient;
}
if ( coefficient == 1 ) {
if ( exponent == 1 ) {
return "x";
}
else {
return "x^" + exponent;
}
}
if ( exponent == 1 ) {
return coefficient + "x";
}
return coefficient + "x^" + exponent;
}
}
} </lang>
- Output:
Task 1: cyclotomic polynomials for n <= 30: CP[1] = x - 1 CP[2] = x + 1 CP[3] = x^2 + x + 1 CP[4] = x^2 + 1 CP[5] = x^4 + x^3 + x^2 + x + 1 CP[6] = x^2 - x + 1 CP[7] = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[8] = x^4 + 1 CP[9] = x^6 + x^3 + 1 CP[10] = x^4 - x^3 + x^2 - x + 1 CP[11] = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[12] = x^4 - x^2 + 1 CP[13] = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[14] = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[15] = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 CP[16] = x^8 + 1 CP[17] = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[18] = x^6 - x^3 + 1 CP[19] = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[20] = x^8 - x^6 + x^4 - x^2 + 1 CP[21] = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 CP[22] = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[23] = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[24] = x^8 - x^4 + 1 CP[25] = x^20 + x^15 + x^10 + x^5 + 1 CP[26] = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 CP[27] = x^18 + x^9 + 1 CP[28] = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 CP[29] = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 CP[30] = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient: CP[1] has coefficient with magnitude = 1 CP[105] has coefficient with magnitude = 2 CP[385] has coefficient with magnitude = 3 CP[1365] has coefficient with magnitude = 4 CP[1785] has coefficient with magnitude = 5 CP[2805] has coefficient with magnitude = 6 CP[3135] has coefficient with magnitude = 7 CP[6545] has coefficient with magnitude = 8 CP[6545] has coefficient with magnitude = 9 CP[10465] has coefficient with magnitude = 10
Julia
<lang julia>using Primes, Polynomials
- memoize cache for recursive calls
const cyclotomics = Dict([1 => Poly([big"-1", big"1"]), 2 => Poly([big"1", big"1"])])
- get all integer divisors of an integer, except itself
function divisors(n::Integer)
f = [one(n)]
for (p, e) in factor(n)
f = reduce(vcat, [f * p^j for j in 1:e], init=f)
end
return resize!(f, length(f) - 1)
end
"""
cyclotomic(n::Integer)
Calculate the n -th cyclotomic polynomial. See wikipedia article at bottom of section /wiki/Cyclotomic_polynomial#Fundamental_tools The algorithm is reliable but slow for large n > 1000. """ function cyclotomic(n::Integer)
if haskey(cyclotomics, n)
c = cyclotomics[n]
elseif isprime(n)
c = Poly(ones(BigInt, n))
cyclotomics[n] = c
else # recursive formula seen in wikipedia article
c = Poly([big"-1"; zeros(BigInt, n - 1); big"1"])
for d in divisors(n)
c ÷= cyclotomic(d)
end
cyclotomics[n] = c
end
return c
end
println("First 30 cyclotomic polynomials:") for i in 1:30
println(rpad("$i: ", 5), cyclotomic(BigInt(i)))
end
const dig = zeros(BigInt, 10) for i in 1:1000000
if all(x -> x != 0, dig)
break
end
for coef in coeffs(cyclotomic(i))
x = abs(coef)
if 0 < x < 11 && dig[Int(x)] == 0
dig[Int(x)] = coef < 0 ? -i : i
end
end
end for (i, n) in enumerate(dig)
println("The cyclotomic polynomial Φ(", abs(n), ") has a coefficient that is ", n < 0 ? -i : i)
end
</lang>
- Output:
First 30 cyclotomic polynomials: 1: Poly(-1 + x) 2: Poly(1 + x) 3: Poly(1 + x + x^2) 4: Poly(1.0 + 1.0*x^2) 5: Poly(1 + x + x^2 + x^3 + x^4) 6: Poly(1.0 - 1.0*x + 1.0*x^2) 7: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6) 8: Poly(1.0 + 1.0*x^4) 9: Poly(1.0 + 1.0*x^3 + 1.0*x^6) 10: Poly(1.0 - 1.0*x + 1.0*x^2 - 1.0*x^3 + 1.0*x^4) 11: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10) 12: Poly(1.0 - 1.0*x^2 + 1.0*x^4) 13: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12) 14: Poly(1.0 - 1.0*x + 1.0*x^2 - 1.0*x^3 + 1.0*x^4 - 1.0*x^5 + 1.0*x^6) 15: Poly(1.0 - 1.0*x + 1.0*x^3 - 1.0*x^4 + 1.0*x^5 - 1.0*x^7 + 1.0*x^8) 16: Poly(1.0 + 1.0*x^8) 17: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12 + x^13 + x^14 + x^15 + x^16) 18: Poly(1.0 - 1.0*x^3 + 1.0*x^6) 19: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12 + x^13 + x^14 + x^15 + x^16 + x^17 + x^18) 20: Poly(1.0 - 1.0*x^2 + 1.0*x^4 - 1.0*x^6 + 1.0*x^8) 21: Poly(1.0 - 1.0*x + 1.0*x^3 - 1.0*x^4 + 1.0*x^6 - 1.0*x^8 + 1.0*x^9 - 1.0*x^11 + 1.0*x^12) 22: Poly(1.0 - 1.0*x + 1.0*x^2 - 1.0*x^3 + 1.0*x^4 - 1.0*x^5 + 1.0*x^6 - 1.0*x^7 + 1.0*x^8 - 1.0*x^9 + 1.0*x^10) 23: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12 + x^13 + x^14 + x^15 + x^16 + x^17 + x^18 + x^19 + x^20 + x^21 + x^22) 24: Poly(1.0 - 1.0*x^4 + 1.0*x^8) 25: Poly(1.0 + 1.0*x^5 + 1.0*x^10 + 1.0*x^15 + 1.0*x^20) 26: Poly(1.0 - 1.0*x + 1.0*x^2 - 1.0*x^3 + 1.0*x^4 - 1.0*x^5 + 1.0*x^6 - 1.0*x^7 + 1.0*x^8 - 1.0*x^9 + 1.0*x^10 - 1.0*x^11 + 1.0*x^12) 27: Poly(1.0 + 1.0*x^9 + 1.0*x^18) 28: Poly(1.0 - 1.0*x^2 + 1.0*x^4 - 1.0*x^6 + 1.0*x^8 - 1.0*x^10 + 1.0*x^12) 29: Poly(1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12 + x^13 + x^14 + x^15 + x^16 + x^17 + x^18 + x^19 + x^20 + x^21 + x^22 + x^23 + x^24 + x^25 + x^26 + x^27 + x^28) 30: Poly(1.0 + 1.0*x - 1.0*x^3 - 1.0*x^4 - 1.0*x^5 + 1.0*x^7 + 1.0*x^8) The cyclotomic polynomial Φ(1) has a coefficient that is -1 The cyclotomic polynomial Φ(105) has a coefficient that is -2 The cyclotomic polynomial Φ(385) has a coefficient that is -3 The cyclotomic polynomial Φ(1365) has a coefficient that is -4 The cyclotomic polynomial Φ(1785) has a coefficient that is 5 The cyclotomic polynomial Φ(2805) has a coefficient that is -6 The cyclotomic polynomial Φ(3135) has a coefficient that is 7 The cyclotomic polynomial Φ(6545) has a coefficient that is -8 The cyclotomic polynomial Φ(6545) has a coefficient that is 9 The cyclotomic polynomial Φ(10465) has a coefficient that is 10
Perl
Conveniently, the module Math::Polynomial::Cyclotomic exists to do all the work. An exponent too large error prevents reaching the 10th step of the 2nd part of the task.
<lang perl>use feature 'say';
use List::Util qw(first);
use Math::Polynomial::Cyclotomic qw(cyclo_poly_iterate);
say 'First 30 cyclotomic polynomials:'; my $it = cyclo_poly_iterate(1); say "$_: " . $it->() for 1 .. 30;
say "\nSmallest cyclotomic polynomial with n or -n as a coefficient:"; $it = cyclo_poly_iterate(1);
for (my ($n, $k) = (1, 1) ; $n <= 10 ; ++$k) {
my $poly = $it->();
while (my $c = first { abs($_) == $n } $poly->coeff) {
say "CP $k has coefficient with magnitude = $n";
$n++;
}
}</lang>
- Output:
First 30 cyclotomic polynomials: 1: (x - 1) 2: (x + 1) 3: (x^2 + x + 1) 4: (x^2 + 1) 5: (x^4 + x^3 + x^2 + x + 1) 6: (x^2 - x + 1) 7: (x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 8: (x^4 + 1) 9: (x^6 + x^3 + 1) 10: (x^4 - x^3 + x^2 - x + 1) 11: (x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 12: (x^4 - x^2 + 1) 13: (x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 14: (x^6 - x^5 + x^4 - x^3 + x^2 - x + 1) 15: (x^8 - x^7 + x^5 - x^4 + x^3 - x + 1) 16: (x^8 + 1) 17: (x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 18: (x^6 - x^3 + 1) 19: (x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 20: (x^8 - x^6 + x^4 - x^2 + 1) 21: (x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1) 22: (x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1) 23: (x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 24: (x^8 - x^4 + 1) 25: (x^20 + x^15 + x^10 + x^5 + 1) 26: (x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1) 27: (x^18 + x^9 + 1) 28: (x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1) 29: (x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) 30: (x^8 + x^7 - x^5 - x^4 - x^3 + x + 1) Smallest cyclotomic polynomial with n or -n as a coefficient: CP 1 has coefficient with magnitude = 1 CP 105 has coefficient with magnitude = 2 CP 385 has coefficient with magnitude = 3 CP 1365 has coefficient with magnitude = 4 CP 1785 has coefficient with magnitude = 5 CP 2805 has coefficient with magnitude = 6 CP 3135 has coefficient with magnitude = 7 CP 6545 has coefficient with magnitude = 8 CP 6545 has coefficient with magnitude = 9
Perl 6
Uses the same library as Perl, so comes with the same caveats. <lang perl6>use Math::Polynomial::Cyclotomic:from<Perl5> <cyclo_poly_iterate cyclo_poly>;
say 'First 30 cyclotomic polynomials:'; my $iterator = cyclo_poly_iterate(1); say "Φ($_) = " ~ super $iterator().Str for 1..30;
say "\nSmallest cyclotomic polynomial with |n| as a coefficient:"; say "Φ(1) has a coefficient magnitude: 1";
my $index = 0; for 2..9 -> $coefficient {
loop {
$index += 5;
my \Φ = cyclo_poly($index);
next unless Φ ~~ / $coefficient\* /;
say "Φ($index) has a coefficient magnitude: $coefficient";
$index -= 5;
last;
}
}
sub super ($str) {
$str.subst( / '^' (\d+) /, { $0.trans([<0123456789>.comb] => [<⁰¹²³⁴⁵⁶⁷⁸⁹>.comb]) }, :g)
}</lang>
First 30 cyclotomic polynomials: Φ(1) = (x - 1) Φ(2) = (x + 1) Φ(3) = (x² + x + 1) Φ(4) = (x² + 1) Φ(5) = (x⁴ + x³ + x² + x + 1) Φ(6) = (x² - x + 1) Φ(7) = (x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(8) = (x⁴ + 1) Φ(9) = (x⁶ + x³ + 1) Φ(10) = (x⁴ - x³ + x² - x + 1) Φ(11) = (x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(12) = (x⁴ - x² + 1) Φ(13) = (x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(14) = (x⁶ - x⁵ + x⁴ - x³ + x² - x + 1) Φ(15) = (x⁸ - x⁷ + x⁵ - x⁴ + x³ - x + 1) Φ(16) = (x⁸ + 1) Φ(17) = (x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(18) = (x⁶ - x³ + 1) Φ(19) = (x¹⁸ + x¹⁷ + x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(20) = (x⁸ - x⁶ + x⁴ - x² + 1) Φ(21) = (x¹² - x¹¹ + x⁹ - x⁸ + x⁶ - x⁴ + x³ - x + 1) Φ(22) = (x¹⁰ - x⁹ + x⁸ - x⁷ + x⁶ - x⁵ + x⁴ - x³ + x² - x + 1) Φ(23) = (x²² + x²¹ + x²⁰ + x¹⁹ + x¹⁸ + x¹⁷ + x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(24) = (x⁸ - x⁴ + 1) Φ(25) = (x²⁰ + x¹⁵ + x¹⁰ + x⁵ + 1) Φ(26) = (x¹² - x¹¹ + x¹⁰ - x⁹ + x⁸ - x⁷ + x⁶ - x⁵ + x⁴ - x³ + x² - x + 1) Φ(27) = (x¹⁸ + x⁹ + 1) Φ(28) = (x¹² - x¹⁰ + x⁸ - x⁶ + x⁴ - x² + 1) Φ(29) = (x²⁸ + x²⁷ + x²⁶ + x²⁵ + x²⁴ + x²³ + x²² + x²¹ + x²⁰ + x¹⁹ + x¹⁸ + x¹⁷ + x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² + x¹¹ + x¹⁰ + x⁹ + x⁸ + x⁷ + x⁶ + x⁵ + x⁴ + x³ + x² + x + 1) Φ(30) = (x⁸ + x⁷ - x⁵ - x⁴ - x³ + x + 1) Smallest cyclotomic polynomial with |n| as a coefficient: Φ(1) has a coefficient magnitude: 1 Φ(105) has a coefficient magnitude: 2 Φ(385) has a coefficient magnitude: 3 Φ(1365) has a coefficient magnitude: 4 Φ(1785) has a coefficient magnitude: 5 Φ(2805) has a coefficient magnitude: 6 Φ(3135) has a coefficient magnitude: 7 Φ(6545) has a coefficient magnitude: 8 Φ(6545) has a coefficient magnitude: 9
Phix
Uses several routines from Polynomial_long_division#Phix, tweaked slightly to check remainder is zero and trim the quotient. <lang Phix>-- demo\rosetta\Cyclotomic_Polynomial.exw function degree(sequence p)
for i=length(p) to 1 by -1 do
if p[i]!=0 then return i end if
end for
return -1
end function
function poly_div(sequence n, d)
while length(d)<length(n) do d &=0 end while
integer dn = degree(n),
dd = degree(d)
if dd<0 then throw("divide by zero") end if
sequence quot = repeat(0,dn)
while dn>=dd do
integer k = dn-dd
integer qk = n[dn]/d[dd]
quot[k+1] = qk
sequence d2 = d[1..length(d)-k]
for i=1 to length(d2) do
n[-i] -= d2[-i]*qk
end for
dn = degree(n)
end while
-- return {quot,n} -- (n is now the remainder)
if n!=repeat(0,length(n)) then ?9/0 end if while quot[$]=0 do quot = quot[1..$-1] end while return quot
end function
function poly(sequence si) -- display helper
string r = ""
for t=length(si) to 1 by -1 do
integer sit = si[t]
if sit!=0 then
if sit=1 and t>1 then
r &= iff(r=""? "":" + ")
elsif sit=-1 and t>1 then
r &= iff(r=""?"-":" - ")
else
if r!="" then
r &= iff(sit<0?" - ":" + ")
sit = abs(sit)
end if
r &= sprintf("%d",sit)
end if
r &= iff(t>1?"x"&iff(t>2?sprintf("^%d",t-1):""):"")
end if
end for
if r="" then r="0" end if
return r
end function --</Polynomial_long_division.exw>
--# memoize cache for recursive calls constant cyclotomics = new_dict({{1,{-1,1}},{2,{1,1}}})
function cyclotomic(integer n) -- -- Calculate nth cyclotomic polynomial. -- See wikipedia article at bottom of section /wiki/Cyclotomic_polynomial#Fundamental_tools -- The algorithm is reliable but slow for large n > 1000. --
sequence c
if getd_index(n,cyclotomics)!=NULL then
c = getd(n,cyclotomics)
else
if is_prime(n) then
c = repeat(1,n)
else -- recursive formula seen in wikipedia article
c = -1&repeat(0,n-1)&1
sequence f = factors(n,-1)
for i=1 to length(f) do
c = poly_div(c,cyclotomic(f[i]))
end for
end if
setd(n,c,cyclotomics)
end if
return c
end function
for i=1 to 30 do
sequence z = cyclotomic(i)
string s = poly(z)
printf(1,"cp(%2d) = %s\n",{i,s})
if i>1 and z!=reverse(z) then ?9/0 end if -- sanity check
end for
integer found = 0, n = 1, cheat = 0 sequence fn = repeat(false,10),
nxt = {105,385,1365,1785,2805,3135,6545,6545,10465,10465}
atom t1 = time()+1 puts(1,"\n") while found<10 do
sequence z = cyclotomic(n)
for i=1 to length(z) do
atom azi = abs(z[i])
if azi>=1 and azi<=10 and fn[azi]=0 then
printf(1,"cp(%d) has a coefficient with magnitude %d\n",{n,azi})
cheat = azi -- (comment this out to prevent cheating!)
found += 1
fn[azi] = true
t1 = time()+1
end if
end for
if cheat then {n,cheat} = {nxt[cheat],0} else n += iff(n=1?4:10) end if
if time()>t1 then
printf(1,"working (%d) ...\r",n)
t1 = time()+1
end if
end while</lang>
- Output:
If you disable the cheating, and if in a particularly self harming mood replace it with n+=1, you will get exactly the same output, eventually.
(The distributed version contains simple instrumentation showing cp(1260) executes the line in the heart of poly_div() that subtracts a multiple of qk over 15 million times.)
cp( 1) = x - 1 cp( 2) = x + 1 cp( 3) = x^2 + x + 1 cp( 4) = x^2 + 1 cp( 5) = x^4 + x^3 + x^2 + x + 1 cp( 6) = x^2 - x + 1 cp( 7) = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp( 8) = x^4 + 1 cp( 9) = x^6 + x^3 + 1 cp(10) = x^4 - x^3 + x^2 - x + 1 cp(11) = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(12) = x^4 - x^2 + 1 cp(13) = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(14) = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 cp(15) = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 cp(16) = x^8 + 1 cp(17) = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(18) = x^6 - x^3 + 1 cp(19) = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(20) = x^8 - x^6 + x^4 - x^2 + 1 cp(21) = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 cp(22) = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 cp(23) = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(24) = x^8 - x^4 + 1 cp(25) = x^20 + x^15 + x^10 + x^5 + 1 cp(26) = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 cp(27) = x^18 + x^9 + 1 cp(28) = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 cp(29) = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 cp(30) = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 cp(1) has a coefficient with magnitude 1 cp(105) has a coefficient with magnitude 2 cp(385) has a coefficient with magnitude 3 cp(1365) has a coefficient with magnitude 4 cp(1785) has a coefficient with magnitude 5 cp(2805) has a coefficient with magnitude 6 cp(3135) has a coefficient with magnitude 7 cp(6545) has a coefficient with magnitude 8 cp(6545) has a coefficient with magnitude 9 cp(10465) has a coefficient with magnitude 10
Python
<lang python>from itertools import count from collections import deque
def primes(_cache=[2, 3]):
yield from _cache
for n in count(_cache[-1]+2, 2):
if isprime(n):
_cache.append(n)
yield n
def isprime(n):
for p in primes():
if n%p == 0:
return False
if p*p > n:
return True
def factors(n):
for p in primes():
if p*p > n:
if n > 1:
yield(n, 1, 1)
break
if n%p == 0:
cnt = 0
while True:
n, cnt = n//p, cnt+1
if n%p != 0: break
yield p, cnt, n
- ^^ not the most sophiscated prime number routines, because no need
def cyclotomic(n):
# Returns (list1, list2) representing the division between
# two polimoials. A list p of integers means the product
# (x^p[0] - 1) * (x^p[1] - 1) * ...
if n == 0:
return ([], [])
if n == 1:
return ([1], [])
p, m, r = next(factors(n)) poly = cyclotomic(r) return elevate(poly_div(elevate(poly, p), poly), p**(m-1))
def poly_div(num, den):
return (num[0] + den[1], num[1] + den[0])
def to_text(poly):
return .join(['%+d x^%d' % x for x in poly])
def elevate(poly, n): # replace poly p(x) with p(x**n)
def powerup(p, n):
return [a*n for a in p]
return poly if n == 1 else (powerup(poly[0], n), powerup(poly[1], n))
def terms(poly):
# convert above representation of divisions to (coef, exponent) pairs
def merge(a, b):
# a, b should be deques. They may change during the course.
while a or b:
l = a[0] if a else (0, -1) # sentinel value
r = b[0] if b else (0, -1)
if l[1] > r[1]:
a.popleft()
elif l[1] < r[1]:
b.popleft()
l = r
else:
a.popleft()
b.popleft()
l = (l[0] + r[0], l[1])
yield l
def mul(poly, p): # p means polynomial x^p - 1
poly = list(poly)
return merge(deque((c, e+p) for c,e in poly),
deque((-c, e) for c,e in poly))
def div(poly, p): # p means polynomial x^p - 1
q = deque()
for c,e in merge(deque(poly), q):
if c:
q.append((c, e - p))
yield (c, e - p)
if e == p: break
p = [(1, 0)] # 1*x^0, i.e. 1
for x in poly[0]: # numerator
p = mul(p, x)
for x in sorted(poly[1], reverse=True): # denominator
p = div(p, x)
return p
for n in range(11):
print(f'{n}: {to_text(terms(cyclotomic(n)))}')
want = 1 for n in count():
c = [c for c,_ in terms(cyclotomic(n))]
while want in c or -want in c:
print(f'C[{want}]: {n}')
want += 1</lang>
- Output:
Only showing first 10 polynomials to avoid clutter.
0: +1 x^0 1: +1 x^1-1 x^0 2: +1 x^1+1 x^0 3: +1 x^2+1 x^1+1 x^0 4: +1 x^2+1 x^0 5: +1 x^4+1 x^3+1 x^2+1 x^1+1 x^0 6: +1 x^2-1 x^1+1 x^0 7: +1 x^6+1 x^5+1 x^4+1 x^3+1 x^2+1 x^1+1 x^0 8: +1 x^4+1 x^0 9: +1 x^6+1 x^3+1 x^0 10: +1 x^4-1 x^3+1 x^2-1 x^1+1 x^0 C[1]: 0 C[2]: 105 C[3]: 385 C[4]: 1365 C[5]: 1785 C[6]: 2805 C[7]: 3135 C[8]: 6545 C[9]: 6545 C[10]: 10465 C[11]: 10465 C[12]: 10465 C[13]: 10465 C[14]: 10465 C[15]: 11305 C[16]: 11305 C[17]: 11305 C[18]: 11305 C[19]: 11305 C[20]: 11305 C[21]: 11305 C[22]: 15015 C[23]: 15015 ...
Sidef
Solution based on polynomial interpolation (slow). <lang ruby>var Poly = require('Math::Polynomial') Poly.string_config(Hash(fold_sign => true, prefix => "", suffix => ""))
func poly_interpolation(v) {
v.len.of {|n| v.len.of {|k| n**k } }.msolve(v)
}
say "First 30 cyclotomic polynomials:" for k in (1..30) {
var a = (k+1).of { cyclotomic(k, _) }
var Φ = poly_interpolation(a)
say ("Φ(#{k}) = ", Poly.new(Φ...))
}
say "\nSmallest cyclotomic polynomial with n or -n as a coefficient:" for n in (1..10) { # very slow
var k = (1..Inf -> first {|k|
poly_interpolation((k+1).of { cyclotomic(k, _) }).first { .abs == n }
})
say "Φ(#{k}) has coefficient with magnitude #{n}"
}</lang>
Slightly faster solution, using the Math::Polynomial::Cyclotomic Perl module. <lang ruby>var Poly = require('Math::Polynomial')
require('Math::Polynomial::Cyclotomic')
Poly.string_config(Hash(fold_sign => true, prefix => "", suffix => ""))
say "First 30 cyclotomic polynomials:" for k in (1..30) {
say ("Φ(#{k}) = ", Poly.new.cyclotomic(k))
}
say "\nSmallest cyclotomic polynomial with n or -n as a coefficient:" for n in (1..10) {
var p = Poly.new
var k = (1..Inf -> first {|k|
[p.cyclotomic(k).coeff].first { .abs == n }
})
say "Φ(#{k}) has coefficient with magnitude = #{n}"
}</lang>
- Output:
First 30 cyclotomic polynomials: Φ(1) = x - 1 Φ(2) = x + 1 Φ(3) = x^2 + x + 1 Φ(4) = x^2 + 1 Φ(5) = x^4 + x^3 + x^2 + x + 1 Φ(6) = x^2 - x + 1 Φ(7) = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(8) = x^4 + 1 Φ(9) = x^6 + x^3 + 1 Φ(10) = x^4 - x^3 + x^2 - x + 1 Φ(11) = x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(12) = x^4 - x^2 + 1 Φ(13) = x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(14) = x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 Φ(15) = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1 Φ(16) = x^8 + 1 Φ(17) = x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(18) = x^6 - x^3 + 1 Φ(19) = x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(20) = x^8 - x^6 + x^4 - x^2 + 1 Φ(21) = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1 Φ(22) = x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 Φ(23) = x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(24) = x^8 - x^4 + 1 Φ(25) = x^20 + x^15 + x^10 + x^5 + 1 Φ(26) = x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 Φ(27) = x^18 + x^9 + 1 Φ(28) = x^12 - x^10 + x^8 - x^6 + x^4 - x^2 + 1 Φ(29) = x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 Φ(30) = x^8 + x^7 - x^5 - x^4 - x^3 + x + 1 Smallest cyclotomic polynomial with n or -n as a coefficient: Φ(1) has coefficient with magnitude = 1 Φ(105) has coefficient with magnitude = 2 Φ(385) has coefficient with magnitude = 3 Φ(1365) has coefficient with magnitude = 4 Φ(1785) has coefficient with magnitude = 5 Φ(2805) has coefficient with magnitude = 6 Φ(3135) has coefficient with magnitude = 7 ^C