' 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