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.
-- demo\rosetta\Cyclotomic_Polynomial.exw
with javascript_semantics
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
            integer mi = -i
            n[mi] -= d2[mi]*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 the 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,deep_copy(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<iff(platform()=JS?5: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 and platform()!=JS then
        printf(1,"working (%d) ...\r",n)
        t1 = time()+1
    end if
end while
