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.
- task
- Find and print the first 30 cyclotomic polynomials.
- Find and print the order of the first 10 cyclotomic polynomials that have i or -i 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.
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