#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