' Visual Basic 6 functions to calculate the ModularInverse for a symmetrical leap cycle.
' by Dr. Irv Bromberg, University of Toronto, Canada
'
' 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
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