' Visual Basic 6 functions to calculate the ModularInverse for a symmetrical leap cycle.
' by Dr. Irv Bromberg, University of Toronto, Canada
' <http://www.sym454.org/leap/>
' placed in the public domain on May 12, 2008 (Gregorian)


Public Function ModularInverse(ByVal L As Long, ByVal c As Long) As Long

    ' returns the modular (multiplicative) inverse U of integers L and C
    ' where normally L=leap years per cycle, C=years per cycle
    ' returns U=modular (multiplicative) inverse for leap cycle,
    ' indicates how many years earlier the leap cycle will start for each increment of K in leap rule:
    ' year is leap if and only if (L*Y+K) mod C < L
    ' also used to calculate the earliest and latest years in Helios cycle relative to new year moment
    ' or equinox or solstice as per the Debug.Print at the end of the function
    
    ' http://mathworld.wolfram.com/ModularInverse.html
    ' the modular inverse can be computed in Mathematica using PowerMod[L, -1, C].
    
    ' see also the following, and follow its links to iterative Extended Euclidean Algorithm:
    ' http://en.wikipedia.org/wiki/Modular_multiplicative_inverse

    Dim x As Long
    Dim Y As Long
    
    ' x and y get changed by reference
    ' for valid modular (multiplicative) inverse L and c must be coprime, confirmed by testing for GCD=1
    ' take x modulo c in case it is negative (same as adding it to C if negative)
    ' if GCD<>1 then return zero, which is never a valid result
    If Extended_GCD(L, c, x, Y) = 1 Then ModularInverse = modulus(x, c) Else ModularInverse = 0

End Function


Public Function Extended_GCD(ByVal a As Long, ByVal b As Long, ByRef x As Long, ByRef Y As Long) As Long

    ' iteratively execute the extended Euclidean algorithm
    ' used by ModularInverse() instead of Extended_GCD_recursive(), has the advantage that it returns GCD and is more efficient
    ' source <http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Iterative_method>

    Dim LastX As Long
    Dim LastY As Long
    Dim t As Long   ' temporary
    Dim q As Long   ' quotient

    x = 0
    Y = 1
    LastX = 1
    LastY = 0

    Do While b <> 0
    
        t = b
        q = quotient(a, b)
        b = modulus(a, b)
        a = t
        
        t = x
        x = LastX - q * x
        LastX = t
        
        t = Y
        Y = LastY - q * Y
        LastY = t
        
        ' Remainder = Dividend - Quotient * Divisor
        
    Loop
        
    ' return x and y by reference
    x = LastX
    Y = LastY
    
    ' return GCD as function result
    Extended_GCD = a
    
End Function



Public Function quotient(ByVal numerator As Double, ByVal denominator As Double) As Double

    quotient = floor(numerator / denominator)

End Function



Public Function modulus(ByVal x As Double, ByVal n As Double) As Double
    
    modulus = x - n * floor(x / n)
    
End Function



Public Function floor(ByVal x As Double) As Double

'    floor = Int(x) ' has bug, when x is whole number, sometimes returns number below!
    
    floor = Int(CDbl(CStr(x)))  ' convert x to string, then back to double, then to integer
    ' (a lot of extra trouble, but avoids the bug above)

End Function