484 lines
16 KiB
VB.net
484 lines
16 KiB
VB.net
Imports System.Text
|
|
|
|
Module Module1
|
|
Private ReadOnly MAX_ALL_FACTORS As Integer = 100_000
|
|
#Const ALGORITHM = 2
|
|
|
|
Class Term
|
|
Implements IComparable(Of Term)
|
|
|
|
Public ReadOnly Property Coefficient As Long
|
|
Public ReadOnly Property Exponent As Long
|
|
|
|
Public Sub New(c As Long, Optional e As Long = 0)
|
|
Coefficient = c
|
|
Exponent = e
|
|
End Sub
|
|
|
|
Public Shared Operator -(t As Term) As Term
|
|
Return New Term(-t.Coefficient, t.Exponent)
|
|
End Operator
|
|
|
|
Public Shared Operator +(lhs As Term, rhs As Term) As Term
|
|
If lhs.Exponent <> rhs.Exponent Then
|
|
Throw New ArgumentException("Exponents not equal")
|
|
End If
|
|
Return New Term(lhs.Coefficient + rhs.Coefficient, lhs.Exponent)
|
|
End Operator
|
|
|
|
Public Shared Operator *(lhs As Term, rhs As Term) As Term
|
|
Return New Term(lhs.Coefficient * rhs.Coefficient, lhs.Exponent + rhs.Exponent)
|
|
End Operator
|
|
|
|
Public Function CompareTo(other As Term) As Integer Implements IComparable(Of Term).CompareTo
|
|
Return -Exponent.CompareTo(other.Exponent)
|
|
End Function
|
|
|
|
Public Overrides Function ToString() As String
|
|
If Coefficient = 0 Then
|
|
Return "0"
|
|
End If
|
|
If Exponent = 0 Then
|
|
Return Coefficient.ToString
|
|
End If
|
|
If Coefficient = 1 Then
|
|
If Exponent = 1 Then
|
|
Return "x"
|
|
End If
|
|
Return String.Format("x^{0}", Exponent)
|
|
End If
|
|
If Coefficient = -1 Then
|
|
If Exponent = 1 Then
|
|
Return "-x"
|
|
End If
|
|
Return String.Format("-x^{0}", Exponent)
|
|
End If
|
|
If Exponent = 1 Then
|
|
Return String.Format("{0}x", Coefficient)
|
|
End If
|
|
Return String.Format("{0}x^{1}", Coefficient, Exponent)
|
|
End Function
|
|
End Class
|
|
|
|
Class Polynomial
|
|
Implements IEnumerable(Of Term)
|
|
|
|
Private ReadOnly polyTerms As New List(Of Term)
|
|
|
|
Public Sub New()
|
|
polyTerms.Add(New Term(0))
|
|
End Sub
|
|
|
|
Public Sub New(ParamArray values() As Term)
|
|
If values.Length = 0 Then
|
|
polyTerms.Add(New Term(0))
|
|
Else
|
|
polyTerms.AddRange(values)
|
|
End If
|
|
Normalize()
|
|
End Sub
|
|
|
|
Public Sub New(values As IEnumerable(Of Term))
|
|
polyTerms.AddRange(values)
|
|
If polyTerms.Count = 0 Then
|
|
polyTerms.Add(New Term(0))
|
|
End If
|
|
Normalize()
|
|
End Sub
|
|
|
|
Public Function LeadingCoeficient() As Long
|
|
Return polyTerms(0).Coefficient
|
|
End Function
|
|
|
|
Public Function Degree() As Long
|
|
Return polyTerms(0).Exponent
|
|
End Function
|
|
|
|
Public Function HasCoefficentAbs(coeff As Long) As Boolean
|
|
For Each t In polyTerms
|
|
If Math.Abs(t.Coefficient) = coeff Then
|
|
Return True
|
|
End If
|
|
Next
|
|
Return False
|
|
End Function
|
|
|
|
Public Function GetEnumerator() As IEnumerator(Of Term) Implements IEnumerable(Of Term).GetEnumerator
|
|
Return polyTerms.GetEnumerator
|
|
End Function
|
|
|
|
Private Function IEnumerable_GetEnumerator() As IEnumerator Implements IEnumerable.GetEnumerator
|
|
Return polyTerms.GetEnumerator
|
|
End Function
|
|
|
|
Private Sub Normalize()
|
|
polyTerms.Sort(Function(a As Term, b As Term) a.CompareTo(b))
|
|
End Sub
|
|
|
|
Public Shared Operator +(lhs As Polynomial, rhs As Term) As Polynomial
|
|
Dim terms As New List(Of Term)
|
|
Dim added = False
|
|
For Each ct In lhs
|
|
If ct.Exponent = rhs.Exponent Then
|
|
added = True
|
|
If ct.Coefficient + rhs.Coefficient <> 0 Then
|
|
terms.Add(ct + rhs)
|
|
End If
|
|
Else
|
|
terms.Add(ct)
|
|
End If
|
|
Next
|
|
If Not added Then
|
|
terms.Add(rhs)
|
|
End If
|
|
Return New Polynomial(terms)
|
|
End Operator
|
|
|
|
Public Shared Operator *(lhs As Polynomial, rhs As Term) As Polynomial
|
|
Dim terms As New List(Of Term)
|
|
For Each ct In lhs
|
|
terms.Add(ct * rhs)
|
|
Next
|
|
Return New Polynomial(terms)
|
|
End Operator
|
|
|
|
Public Shared Operator +(lhs As Polynomial, rhs As Polynomial) As Polynomial
|
|
Dim terms As New List(Of Term)
|
|
Dim thisCount = lhs.polyTerms.Count
|
|
Dim polyCount = rhs.polyTerms.Count
|
|
While thisCount > 0 OrElse polyCount > 0
|
|
If thisCount = 0 Then
|
|
Dim polyTerm = rhs.polyTerms(polyCount - 1)
|
|
terms.Add(polyTerm)
|
|
polyCount -= 1
|
|
ElseIf polyCount = 0 Then
|
|
Dim thisTerm = lhs.polyTerms(thisCount - 1)
|
|
terms.Add(thisTerm)
|
|
thisCount -= 1
|
|
Else
|
|
Dim polyTerm = rhs.polyTerms(polyCount - 1)
|
|
Dim thisTerm = lhs.polyTerms(thisCount - 1)
|
|
If thisTerm.Exponent = polyTerm.Exponent Then
|
|
Dim t = thisTerm + polyTerm
|
|
If t.Coefficient <> 0 Then
|
|
terms.Add(t)
|
|
End If
|
|
thisCount -= 1
|
|
polyCount -= 1
|
|
ElseIf thisTerm.Exponent < polyTerm.Exponent Then
|
|
terms.Add(thisTerm)
|
|
thisCount -= 1
|
|
Else
|
|
terms.Add(polyTerm)
|
|
polyCount -= 1
|
|
End If
|
|
End If
|
|
End While
|
|
Return New Polynomial(terms)
|
|
End Operator
|
|
|
|
Public Shared Operator *(lhs As Polynomial, rhs As Polynomial) As Polynomial
|
|
Throw New Exception("Not implemented")
|
|
End Operator
|
|
|
|
Public Shared Operator /(lhs As Polynomial, rhs As Polynomial) As Polynomial
|
|
Dim q As New Polynomial
|
|
Dim r = lhs
|
|
Dim lcv = rhs.LeadingCoeficient
|
|
Dim dv = rhs.Degree
|
|
While r.Degree >= rhs.Degree
|
|
Dim lcr = r.LeadingCoeficient
|
|
Dim s = lcr \ lcv
|
|
Dim t As New Term(s, r.Degree() - dv)
|
|
q += t
|
|
r += rhs * -t
|
|
End While
|
|
Return q
|
|
End Operator
|
|
|
|
Public Overrides Function ToString() As String
|
|
Dim builder As New StringBuilder
|
|
Dim it = polyTerms.GetEnumerator()
|
|
If it.MoveNext Then
|
|
builder.Append(it.Current)
|
|
End If
|
|
While it.MoveNext
|
|
If it.Current.Coefficient < 0 Then
|
|
builder.Append(" - ")
|
|
builder.Append(-it.Current)
|
|
Else
|
|
builder.Append(" + ")
|
|
builder.Append(it.Current)
|
|
End If
|
|
End While
|
|
Return builder.ToString
|
|
End Function
|
|
End Class
|
|
|
|
Function GetDivisors(number As Integer) As List(Of Integer)
|
|
Dim divisors As New List(Of Integer)
|
|
Dim root = CType(Math.Sqrt(number), Long)
|
|
For i = 1 To root
|
|
If number Mod i = 0 Then
|
|
divisors.Add(i)
|
|
Dim div = number \ i
|
|
If div <> i AndAlso div <> number Then
|
|
divisors.Add(div)
|
|
End If
|
|
End If
|
|
Next
|
|
Return divisors
|
|
End Function
|
|
|
|
Private ReadOnly allFactors As New Dictionary(Of Integer, Dictionary(Of Integer, Integer)) From {{2, New Dictionary(Of Integer, Integer) From {{2, 1}}}}
|
|
Function GetFactors(number As Integer) As Dictionary(Of Integer, Integer)
|
|
If allFactors.ContainsKey(number) Then
|
|
Return allFactors(number)
|
|
End If
|
|
|
|
Dim factors As New Dictionary(Of Integer, Integer)
|
|
If number Mod 2 = 0 Then
|
|
Dim factorsDivTwo = GetFactors(number \ 2)
|
|
For Each pair In factorsDivTwo
|
|
If Not factors.ContainsKey(pair.Key) Then
|
|
factors.Add(pair.Key, pair.Value)
|
|
End If
|
|
Next
|
|
If factors.ContainsKey(2) Then
|
|
factors(2) += 1
|
|
Else
|
|
factors.Add(2, 1)
|
|
End If
|
|
If number < MAX_ALL_FACTORS Then
|
|
allFactors.Add(number, factors)
|
|
End If
|
|
Return factors
|
|
End If
|
|
Dim root = CType(Math.Sqrt(number), Long)
|
|
Dim i = 3L
|
|
While i <= root
|
|
If number Mod i = 0 Then
|
|
Dim factorsDivI = GetFactors(number \ i)
|
|
For Each pair In factorsDivI
|
|
If Not factors.ContainsKey(pair.Key) Then
|
|
factors.Add(pair.Key, pair.Value)
|
|
End If
|
|
Next
|
|
If factors.ContainsKey(i) Then
|
|
factors(i) += 1
|
|
Else
|
|
factors.Add(i, 1)
|
|
End If
|
|
If number < MAX_ALL_FACTORS Then
|
|
allFactors.Add(number, factors)
|
|
End If
|
|
Return factors
|
|
End If
|
|
i += 2
|
|
End While
|
|
factors.Add(number, 1)
|
|
If number < MAX_ALL_FACTORS Then
|
|
allFactors.Add(number, factors)
|
|
End If
|
|
Return factors
|
|
End Function
|
|
|
|
Private ReadOnly computedPolynomials As New Dictionary(Of Integer, Polynomial)
|
|
Function CyclotomicPolynomial(n As Integer) As Polynomial
|
|
If computedPolynomials.ContainsKey(n) Then
|
|
Return computedPolynomials(n)
|
|
End If
|
|
|
|
If n = 1 Then
|
|
REM polynomial: x - 1
|
|
Dim p As New Polynomial(New Term(1, 1), New Term(-1))
|
|
computedPolynomials.Add(n, p)
|
|
Return p
|
|
End If
|
|
|
|
Dim factors = GetFactors(n)
|
|
Dim terms As New List(Of Term)
|
|
Dim cyclo As Polynomial
|
|
|
|
If factors.ContainsKey(n) Then
|
|
REM n prime
|
|
For index = 1 To n
|
|
terms.Add(New Term(1, index - 1))
|
|
Next
|
|
|
|
cyclo = New Polynomial(terms)
|
|
computedPolynomials.Add(n, cyclo)
|
|
Return cyclo
|
|
ElseIf factors.Count = 2 AndAlso factors.ContainsKey(2) AndAlso factors(2) = 1 AndAlso factors.ContainsKey(n / 2) AndAlso factors(n / 2) = 1 Then
|
|
REM n = 2p
|
|
Dim prime = n \ 2
|
|
Dim coeff = -1
|
|
|
|
For index = 1 To prime
|
|
coeff *= -1
|
|
terms.Add(New Term(coeff, index - 1))
|
|
Next
|
|
|
|
cyclo = New Polynomial(terms)
|
|
computedPolynomials.Add(n, cyclo)
|
|
Return cyclo
|
|
ElseIf factors.Count = 1 AndAlso factors.ContainsKey(2) Then
|
|
REM n = 2^h
|
|
Dim h = factors(2)
|
|
terms = New List(Of Term) From {
|
|
New Term(1, Math.Pow(2, h - 1)),
|
|
New Term(1)
|
|
}
|
|
|
|
cyclo = New Polynomial(terms)
|
|
computedPolynomials.Add(n, cyclo)
|
|
Return cyclo
|
|
ElseIf factors.Count = 1 AndAlso factors.ContainsKey(n) Then
|
|
REM n = p^k
|
|
Dim p = 0
|
|
Dim k = 0
|
|
For Each it In factors
|
|
p = it.Key
|
|
k = it.Value
|
|
Next
|
|
For index = 1 To p
|
|
terms.Add(New Term(1, (index - 1) * Math.Pow(p, k - 1)))
|
|
Next
|
|
|
|
cyclo = New Polynomial(terms)
|
|
computedPolynomials.Add(n, cyclo)
|
|
Return cyclo
|
|
ElseIf factors.Count = 2 AndAlso factors.ContainsKey(2) Then
|
|
REM n = 2^h * p^k
|
|
Dim p = 0
|
|
For Each it In factors
|
|
If it.Key <> 2 Then
|
|
p = it.Key
|
|
End If
|
|
Next
|
|
|
|
Dim coeff = -1
|
|
Dim twoExp = CType(Math.Pow(2, factors(2) - 1), Long)
|
|
Dim k = factors(p)
|
|
For index = 1 To p
|
|
coeff *= -1
|
|
terms.Add(New Term(coeff, (index - 1) * twoExp * Math.Pow(p, k - 1)))
|
|
Next
|
|
|
|
cyclo = New Polynomial(terms)
|
|
computedPolynomials.Add(n, cyclo)
|
|
Return cyclo
|
|
ElseIf factors.ContainsKey(2) AndAlso (n / 2) Mod 2 = 1 AndAlso n / 2 > 1 Then
|
|
REM CP(2m)[x] = CP(-m)[x], n odd integer > 1
|
|
Dim cycloDiv2 = CyclotomicPolynomial(n \ 2)
|
|
For Each t In cycloDiv2
|
|
If t.Exponent Mod 2 = 0 Then
|
|
terms.Add(t)
|
|
Else
|
|
terms.Add(-t)
|
|
End If
|
|
Next
|
|
|
|
cyclo = New Polynomial(terms)
|
|
computedPolynomials.Add(n, cyclo)
|
|
Return cyclo
|
|
End If
|
|
|
|
#If ALGORITHM = 0 Then
|
|
REM slow - uses basic definition
|
|
Dim divisors = GetDivisors(n)
|
|
REM Polynomial: (x^n - 1)
|
|
cyclo = New Polynomial(New Term(1, n), New Term(-1))
|
|
For Each i In divisors
|
|
Dim p = CyclotomicPolynomial(i)
|
|
cyclo /= p
|
|
Next
|
|
|
|
computedPolynomials.Add(n, cyclo)
|
|
Return cyclo
|
|
#ElseIf ALGORITHM = 1 Then
|
|
REM Faster. Remove Max divisor (and all divisors of max divisor) - only one divide for all divisors of Max Divisor
|
|
Dim divisors = GetDivisors(n)
|
|
Dim maxDivisor = Integer.MinValue
|
|
For Each div In divisors
|
|
maxDivisor = Math.Max(maxDivisor, div)
|
|
Next
|
|
Dim divisorExceptMax As New List(Of Integer)
|
|
For Each div In divisors
|
|
If maxDivisor Mod div <> 0 Then
|
|
divisorExceptMax.Add(div)
|
|
End If
|
|
Next
|
|
|
|
REM Polynomial: ( x^n - 1 ) / ( x^m - 1 ), where m is the max divisor
|
|
cyclo = New Polynomial(New Term(1, n), New Term(-1)) / New Polynomial(New Term(1, maxDivisor), New Term(-1))
|
|
For Each i In divisorExceptMax
|
|
Dim p = CyclotomicPolynomial(i)
|
|
cyclo /= p
|
|
Next
|
|
|
|
computedPolynomials.Add(n, cyclo)
|
|
Return cyclo
|
|
#ElseIf ALGORITHM = 2 Then
|
|
REM Fastest
|
|
REM Let p ; q be primes such that p does not divide n, and q divides n
|
|
REM Then Cp(np)[x] = CP(n)[x^p] / CP(n)[x]
|
|
Dim m = 1
|
|
cyclo = CyclotomicPolynomial(m)
|
|
Dim primes As New List(Of Integer)
|
|
For Each it In factors
|
|
primes.Add(it.Key)
|
|
Next
|
|
primes.Sort()
|
|
For Each prime In primes
|
|
REM CP(m)[x]
|
|
Dim cycloM = cyclo
|
|
REM Compute CP(m)[x^p]
|
|
terms = New List(Of Term)
|
|
For Each t In cyclo
|
|
terms.Add(New Term(t.Coefficient, t.Exponent * prime))
|
|
Next
|
|
cyclo = New Polynomial(terms) / cycloM
|
|
m *= prime
|
|
Next
|
|
REM Now, m is the largest square free divisor of n
|
|
Dim s = n \ m
|
|
REM Compute CP(n)[x] = CP(m)[x^s]
|
|
terms = New List(Of Term)
|
|
For Each t In cyclo
|
|
terms.Add(New Term(t.Coefficient, t.Exponent * s))
|
|
Next
|
|
|
|
cyclo = New Polynomial(terms)
|
|
computedPolynomials.Add(n, cyclo)
|
|
Return cyclo
|
|
#Else
|
|
Throw New Exception("Invalid algorithm")
|
|
#End If
|
|
End Function
|
|
|
|
Sub Main()
|
|
Console.WriteLine("Task 1: cyclotomic polynomials for n <= 30:")
|
|
For i = 1 To 30
|
|
Dim p = CyclotomicPolynomial(i)
|
|
Console.WriteLine("CP[{0}] = {1}", i, p)
|
|
Next
|
|
Console.WriteLine()
|
|
|
|
Console.WriteLine("Task 2: Smallest cyclotomic polynomial with n or -n as a coefficient:")
|
|
Dim n = 0
|
|
For i = 1 To 10
|
|
While True
|
|
n += 1
|
|
Dim cyclo = CyclotomicPolynomial(n)
|
|
If cyclo.HasCoefficentAbs(i) Then
|
|
Console.WriteLine("CP[{0}] has coefficient with magnitude = {1}", n, i)
|
|
n -= 1
|
|
Exit While
|
|
End If
|
|
End While
|
|
Next
|
|
End Sub
|
|
|
|
End Module
|