Anonymous user
Cyclotomic polynomial: Difference between revisions
→{{header|Kotlin}}
Thundergnat (talk | contribs) m (Thundergnat moved page Cyclotomic Polynomial to Cyclotomic polynomial: Conway's Game of Life) |
|||
Line 3,022:
CP[6545] has coefficient with magnitude = 9
CP[10465] has coefficient with magnitude = 10</pre>
=={{header|Nim}}==
{{trans|Java}}
We use Java algorithm with ideas from C#, D, Go and Kotlin. We have kept only algorithm number 2 as other algorithms are much less efficient. We have also done some Nim specific improvements in order to get better performances.
<lang Nim>import algorithm, math, sequtils, strformat, tables
type
Term = tuple[coeff: int; exp: Natural]
Polynomial = seq[Term]
# Table used to represent the list of factors of a number.
# If, for a number "n", "k" is present in the table "f" of its factors,
# "f[k]" contains the exponent of "k" in the prime factor decomposition.
Factors = Table[int, int]
####################################################################################################
# Miscellaneous.
## Parity tests.
template isOdd(n: int): bool = (n and 1) != 0
template isEven(n: int): bool = (n and 1) == 0
#---------------------------------------------------------------------------------------------------
proc sort(poly: var Polynomial) {.inline.} =
## Sort procedure for the terms of a polynomial (high degree first).
algorithm.sort(poly, proc(x, y: Term): int = cmp(x.exp, y.exp), Descending)
####################################################################################################
# Superscripts and subscripts.
const
Superscripts: array['0'..'9', string] = ["⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹"]
Subscripts: array['0'..'9', string] = ["₀", "₁", "₂", "₃", "₄", "₅", "₆", "₇", "₈", "₉"]
func superscript(n: Natural): string =
## Return the Unicode string to use to represent an exponent.
if n == 1:
return ""
for d in $n:
result.add(Superscripts[d])
#---------------------------------------------------------------------------------------------------
func subscript(n: Natural): string =
## Return the Unicode string to use to represent an index.
for d in $n:
result.add(Subscripts[d])
####################################################################################################
# Term operations.
func term(coeff, exp: int): Term =
## Create a term.
if exp < 0:
raise newException(ValueError, "term exponent cannot be negative")
(coeff, Natural exp)
#---------------------------------------------------------------------------------------------------
func `*`(a, b: Term): Term =
## Multiply two terms.
(a.coeff * b.coeff, Natural a.exp + b.exp)
#---------------------------------------------------------------------------------------------------
func `+`(a, b: Term): Term =
## Add two terms.
if a.exp != b.exp:
raise newException(ValueError, "addition of terms with unequal exponents")
(a.coeff + b.coeff, a.exp)
#---------------------------------------------------------------------------------------------------
func `-`(a: Term): Term =
## Return the opposite of a term.
(-a.coeff, a.exp)
#---------------------------------------------------------------------------------------------------
func `$`(a: Term): string =
## Return the string representation of a term.
if a.coeff == 0: "0"
elif a.exp == 0: $a.coeff
elif a.coeff == 1: 'x' & superscript(a.exp)
elif a.coeff == -1: "-x" & superscript(a.exp)
else: $a.coeff & 'x' & superscript(a.exp)
####################################################################################################
# Polynomial.
func polynomial(terms: varargs[Term]): Polynomial =
## Create a polynomial described by its terms.
for t in terms:
if t.coeff != 0:
result.add(t)
if result.len == 0:
return @[term(0, 0)]
sort(result)
#---------------------------------------------------------------------------------------------------
func hasCoeffAbs(poly: Polynomial; coeff: int): bool =
## Return true if the polynomial contains a given coefficient.
for t in poly:
if abs(t.coeff) == coeff:
return true
#---------------------------------------------------------------------------------------------------
func leadingCoeff(poly: Polynomial): int {.inline.} =
## Return the coefficient of the term with the highest degree.
poly[0].coeff
#---------------------------------------------------------------------------------------------------
func degree(poly: Polynomial): int {.inline.} =
## Return the degree of the polynomial.
if poly.len == 0: -1
else: poly[0].exp
#---------------------------------------------------------------------------------------------------
func `+`(poly: Polynomial; someTerm: Term): Polynomial =
## Add a term to a polynomial.
var added = false
for currTerm in poly:
if currterm.exp == someTerm.exp:
added = true
if currTerm.coeff + someTerm.coeff != 0:
result.add(currTerm + someTerm)
else:
result.add(currTerm)
if not added:
result.add(someTerm)
#---------------------------------------------------------------------------------------------------
func `+`(a, b: Polynomial): Polynomial =
## Add two polynomials.
var aIndex = a.high
var bIndex = b.high
while aIndex >= 0 or bIndex >= 0:
if aIndex < 0:
result &= b[bIndex]
dec bIndex
elif bIndex < 0:
result &= a[aIndex]
dec aIndex
else:
let t1 = a[aIndex]
let t2 = b[bIndex]
if t1.exp == t2.exp:
let t3 = t1 + t2
if t3.coeff != 0:
result.add(t3)
dec aIndex
dec bIndex
elif t1.exp < t2.exp:
result.add(t1)
dec aIndex
else:
result.add(t2)
dec bIndex
sort(result)
#---------------------------------------------------------------------------------------------------
func `*`(poly: Polynomial; someTerm: Term): Polynomial =
## Multiply a polynomial by a term.
for currTerm in poly:
result.add(currTerm * someTerm)
#---------------------------------------------------------------------------------------------------
func `/`(a, b: Polynomial): Polynomial =
## Divide a polynomial by another polynomial.
var a = a
let lcb = b.leadingCoeff
let db = b.degree
while a.degree >= b.degree:
let lca = a.leadingCoeff
let s = lca div lcb
let t = term(s, a.degree - db)
result = result + t
a = a + b * -t
#---------------------------------------------------------------------------------------------------
func `$`(poly: Polynomial): string =
## Return the string representation of a polynomial.
for t in poly:
if result.len == 0:
result.add($t)
else:
if t.coeff > 0:
result.add('+')
result.add($t)
else:
result.add('-')
result.add($(-t))
####################################################################################################
# Cyclotomic polynomial.
var
# Cache of list of factors.
factorCache: Table[int, Factors] = {2: {2: 1}.toTable}.toTable
# Cache of cyclotomic polynomials. Initialized with 1 -> x - 1.
polyCache: Table[int, Polynomial] = {1: polynomial(term(1, 1), term(-1, 0))}.toTable
#---------------------------------------------------------------------------------------------------
proc getFactors(n: int): Factors =
## Return the list of factors of a number.
if n in factorCache:
return factorCache[n]
if n.isEven:
result = getFactors(n div 2)
result[2] = result.getOrDefault(2) + 1
factorCache[n] = result
return
var i = 3
while i * i <= n:
if n mod i == 0:
result = getFactors( n div i)
result[i] = result.getOrDefault(i) + 1
factorCache[n] = result
return
inc i, 2
result[n] = 1
factorCache[n] = result
#---------------------------------------------------------------------------------------------------
proc cycloPoly(n: int): Polynomial =
## Find the nth cyclotomic polynomial.
if n in polyCache:
return polyCache[n]
let factors = getFactors(n)
if n in factors:
# n is prime.
for i in countdown(n - 1, 0): # Add the terms by decreasing degrees.
result.add(term(1, i))
elif factors.len == 2 and factors.getOrDefault(2) == 1 and factors.getOrDefault(n div 2) == 1:
# n = 2 x prime.
let prime = n div 2
var coeff = -1
for i in countdown(prime - 1, 0): # Add the terms by decreasing degrees.
coeff *= -1
result.add(term(coeff, i))
elif factors.len == 1 and 2 in factors:
# n = 2 ^ h.
let h = factors[2]
result.add([term(1, 1 shl (h - 1)), term(1, 0)])
elif factors.len == 1 and n notin factors:
# n = prime ^ k.
var p, k = 0
for prime, v in factors.pairs:
if prime > p:
p = prime
k = v
for i in countdown(p - 1, 0): # Add the terms by decreasing degrees.
result.add(term(1, i * p^(k-1)))
elif factors.len == 2 and 2 in factors:
# n = 2 ^ h x prime ^ k.
var p, k = 0
for prime, v in factors.pairs:
if prime != 2 and prime > p:
p = prime
k = v
var coeff = -1
let twoExp = 1 shl (factors[2] - 1)
for i in countdown(p - 1, 0): # Add the terms by decreasing degrees.
coeff *= -1
result.add(term(coeff, i * twoExp * p^(k-1)))
elif 2 in factors and isOdd(n div 2) and n div 2 > 1:
# CP(2m)[x] = CP(-m)[x], n odd integer > 1.
let cycloDiv2 = cycloPoly(n div 2)
for t in cycloDiv2:
result.add(if t.exp.isEven: t else: -t)
else:
# 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].
var m = 1
var cyclo = cycloPoly(m)
let primes = sorted(toSeq(factors.keys))
for prime in primes:
# Compute CP(m)[x^p].
var terms: Polynomial
for t in cyclo:
terms.add(term(t.coeff, t.exp * prime))
cyclo = terms / cyclo
m *= prime
# Now, m is the largest square free divisor of n.
let s = n div m
# Compute CP(n)[x] = CP(m)[x^s].
for t in cyclo:
result.add(term(t.coeff, t.exp * s))
polyCache[n] = result
#———————————————————————————————————————————————————————————————————————————————————————————————————
echo "Cyclotomic polynomials for n ⩽ 30:"
for i in 1..30:
echo &"Φ{subscript(i):<2} = {cycloPoly(i)}"
echo ""
echo "Smallest cyclotomic polynomial with n or -n as a coefficient:"
var n = 0
for i in 1..10:
while true:
inc n
let cyclo = cycloPoly(n)
if cyclo.hasCoeffAbs(i):
echo &"Φ{subscript(n):6} has coefficient with magnitude = {i}"
dec n
break</lang>
{{out}}
The program runs in 41 seconds on our reasonably performing laptop.
<pre>Cyclotomic polynomials for n ⩽ 30:
Φ₁ = x-1
Φ₂ = x+1
Φ₃ = x²+x+1
Φ₄ = x²+1
Φ₅ = x⁴+x³+x²+x+1
Φ₆ = x²-x+1
Φ₇ = x⁶+x⁵+x⁴+x³+x²+x+1
Φ₈ = x⁴+1
Φ₉ = x⁶+x³+1
Φ₁₀ = x⁴-x³+x²-x+1
Φ₁₁ = x¹⁰+x⁹+x⁸+x⁷+x⁶+x⁵+x⁴+x³+x²+x+1
Φ₁₂ = x⁴-x²+1
Φ₁₃ = x¹²+x¹¹+x¹⁰+x⁹+x⁸+x⁷+x⁶+x⁵+x⁴+x³+x²+x+1
Φ₁₄ = x⁶-x⁵+x⁴-x³+x²-x+1
Φ₁₅ = x⁸-x⁷+x⁵-x⁴+x³-x+1
Φ₁₆ = x⁸+1
Φ₁₇ = x¹⁶+x¹⁵+x¹⁴+x¹³+x¹²+x¹¹+x¹⁰+x⁹+x⁸+x⁷+x⁶+x⁵+x⁴+x³+x²+x+1
Φ₁₈ = x⁶-x³+1
Φ₁₉ = x¹⁸+x¹⁷+x¹⁶+x¹⁵+x¹⁴+x¹³+x¹²+x¹¹+x¹⁰+x⁹+x⁸+x⁷+x⁶+x⁵+x⁴+x³+x²+x+1
Φ₂₀ = x⁸-x⁶+x⁴-x²+1
Φ₂₁ = x¹²-x¹¹+x⁹-x⁸+x⁶-x⁴+x³-x+1
Φ₂₂ = x¹⁰-x⁹+x⁸-x⁷+x⁶-x⁵+x⁴-x³+x²-x+1
Φ₂₃ = x²²+x²¹+x²⁰+x¹⁹+x¹⁸+x¹⁷+x¹⁶+x¹⁵+x¹⁴+x¹³+x¹²+x¹¹+x¹⁰+x⁹+x⁸+x⁷+x⁶+x⁵+x⁴+x³+x²+x+1
Φ₂₄ = x⁸-x⁴+1
Φ₂₅ = x²⁰+x¹⁵+x¹⁰+x⁵+1
Φ₂₆ = x¹²-x¹¹+x¹⁰-x⁹+x⁸-x⁷+x⁶-x⁵+x⁴-x³+x²-x+1
Φ₂₇ = x¹⁸+x⁹+1
Φ₂₈ = x¹²-x¹⁰+x⁸-x⁶+x⁴-x²+1
Φ₂₉ = 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
Φ₃₀ = x⁸+x⁷-x⁵-x⁴-x³+x+1
Smallest cyclotomic polynomial with n or -n as a coefficient:
Φ₁ has coefficient with magnitude = 1
Φ₁₀₅ has coefficient with magnitude = 2
Φ₃₈₅ has coefficient with magnitude = 3
Φ₁₃₆₅ has coefficient with magnitude = 4
Φ₁₇₈₅ has coefficient with magnitude = 5
Φ₂₈₀₅ has coefficient with magnitude = 6
Φ₃₁₃₅ has coefficient with magnitude = 7
Φ₆₅₄₅ has coefficient with magnitude = 8
Φ₆₅₄₅ has coefficient with magnitude = 9
Φ₁₀₄₆₅ has coefficient with magnitude = 10</pre>
=={{header|Perl}}==
| |||