137 lines
3.8 KiB
Text
137 lines
3.8 KiB
Text
#include "isprime.bas"
|
|
|
|
Type IntArray
|
|
Dim values(Any) As Integer
|
|
Dim length As Integer
|
|
End Type
|
|
|
|
Function distinctPrimeFactors(n As Integer) As IntArray
|
|
Dim result As IntArray
|
|
Redim result.values(0)
|
|
result.length = 0
|
|
|
|
For i As Integer = 2 To n
|
|
If n Mod i = 0 Andalso isPrime(i) Then
|
|
result.length += 1
|
|
Redim Preserve result.values(result.length - 1)
|
|
result.values(result.length - 1) = i
|
|
While n Mod i = 0
|
|
n \= i
|
|
Wend
|
|
End If
|
|
Next
|
|
Return result
|
|
End Function
|
|
|
|
Function substituteExponent(polynomial As IntArray, exponent As Integer) As IntArray
|
|
Dim result As IntArray
|
|
result.length = exponent * (polynomial.length - 1) + 1
|
|
Redim result.values(result.length - 1)
|
|
|
|
For i As Integer = polynomial.length - 1 To 0 Step -1
|
|
result.values(i * exponent) = polynomial.values(i)
|
|
Next
|
|
|
|
Return result
|
|
End Function
|
|
|
|
Function exactDivision(dividend As IntArray, divisor As IntArray) As IntArray
|
|
Dim As Integer i, j
|
|
Dim result As IntArray
|
|
result.length = dividend.length - divisor.length + 1
|
|
Redim result.values(result.length - 1)
|
|
|
|
Dim temp(dividend.length - 1) As Integer
|
|
For i = 0 To dividend.length - 1
|
|
temp(i) = dividend.values(i)
|
|
Next
|
|
|
|
For i = 0 To dividend.length - divisor.length
|
|
result.values(i) = temp(i)
|
|
If temp(i) <> 0 Then
|
|
For j = 1 To divisor.length - 1
|
|
temp(i + j) -= divisor.values(j) * temp(i)
|
|
Next
|
|
End If
|
|
Next
|
|
|
|
Return result
|
|
End Function
|
|
|
|
Function cycloPoly(cpIndex As Integer) As IntArray
|
|
Dim i As Integer
|
|
Dim polynomial As IntArray
|
|
polynomial.length = 2
|
|
Redim polynomial.values(1)
|
|
polynomial.values(0) = 1
|
|
polynomial.values(1) = -1
|
|
|
|
If cpIndex = 1 Then Return polynomial
|
|
|
|
If isPrime(cpIndex) Then
|
|
Dim result As IntArray
|
|
result.length = cpIndex
|
|
Redim result.values(cpIndex - 1)
|
|
For i = 0 To cpIndex - 1
|
|
result.values(i) = 1
|
|
Next
|
|
Return result
|
|
End If
|
|
|
|
Dim primes As IntArray = distinctPrimeFactors(cpIndex)
|
|
Dim product As Integer = 1
|
|
|
|
For i = 0 To primes.length - 1
|
|
Dim numerator As IntArray = substituteExponent(polynomial, primes.values(i))
|
|
polynomial = exactDivision(numerator, polynomial)
|
|
product *= primes.values(i)
|
|
Next
|
|
|
|
Return substituteExponent(polynomial, cpIndex \ product)
|
|
End Function
|
|
|
|
Function hasHeight(polynomial As IntArray, coefficient As Integer) As Boolean
|
|
For i As Integer = 0 To (polynomial.length + 1) \ 2 - 1
|
|
If Abs(polynomial.values(i)) = coefficient Then Return True
|
|
Next
|
|
Return False
|
|
End Function
|
|
|
|
' Main program
|
|
Print "Task 1: Cyclotomic polynomials for n <= 30:"
|
|
Print "CP( 1) = x - 1"
|
|
|
|
For cpIndex As Integer = 2 To 30
|
|
Print Using "CP(##) = "; cpIndex;
|
|
Dim poly As IntArray = cycloPoly(cpIndex)
|
|
|
|
Dim first As Boolean = True
|
|
For i As Integer = poly.length - 1 To 0 Step -1
|
|
If poly.values(i) <> 0 Then
|
|
If Not first Then Print Iif(poly.values(i) > 0, " + ", " ");
|
|
If poly.values(i) <> 1 Or i = 0 Then
|
|
Print Iif(poly.values(i) = -1 And i > 0, "- ", Str(poly.values(i)));
|
|
End If
|
|
|
|
If i > 0 Then
|
|
Print "x";
|
|
If i > 1 Then Print "^" & i;
|
|
End If
|
|
first = False
|
|
End If
|
|
Next
|
|
Print
|
|
Next
|
|
|
|
Print !"\nTask 2: Smallest cyclotomic polynomial with n or -n as a coefficient:"
|
|
Print "CP( 1) has a coefficient with magnitude 1"
|
|
|
|
Dim cpIndex As Integer = 2
|
|
For coeff As Integer = 2 To 10
|
|
While isPrime(cpIndex) Or Not hasHeight(cycloPoly(cpIndex), coeff)
|
|
cpIndex += 1
|
|
Wend
|
|
Print Using "CP(#####) has a coefficient with magnitude &"; cpIndex; coeff
|
|
Next
|
|
|
|
Sleep
|