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



Task

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


See also
 Wikipedia article, Cyclotomic polynomial, showing ways to calculate them.
 The sequence A013594 with the smallest order of cyclotomic polynomial containing n or -n as a coefficient.
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
            }
        }
    }
}