Cyclotomic polynomial

Cyclotomic polynomial
You are encouraged to solve this task according to the task description, using any language you may know.

The nth Cyclotomic polynomial, for any positive integer n, is the unique irreducible polynomial of largest degree with integer coefficients that is a divisor of x^n − 1, and is not a divisor of x^k − 1 for any k < n.

task

• Find and print the first 30 cyclotomic polynomials.
• Find and print the order of the first 10 cyclotomic polynomials that have i or -i as a coefficient.

See also

• wikipedia Cyclotomic_polynomial showing ways to calculate
• 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