Euler Problem 188
Euler Problem 188 asks:
The hyperexponentiation or tetration of a number a by a positive integer b, denoted by a^^b or ba, is recursively defined by:
a^^1 = a,
a^^(k+1) = a(a^^k).Thus we have e.g. 3^^2 = 33 = 27, hence 3^^3 = 327 = 7625597484987 and 3^^4 is roughly 103.6383346400240996*10^12.
Find the last 8 digits of 1777^^1855.
Euler uses double up-arrows for hyperexponentiation. I substituted double carets as a “reasonable facsimile.” Tetration is covered by this Wikipedia article. A key point is to note that tetration is not associative, and we must evaluate the expression from right to left (top to bottom).
This is the recursive version:
If k = 0 Then
HyperExp = 1
Exit Function
ElseIf k = 1 Then
HyperExp = a
Exit Function
End If
HyperExp = a ^ HyperExp(a, k - 1)
End Function
This works fine, but it won’t handle 1777^^1885. Python has a Pow(b,e,m) function that returns base b raised to exponent e modulo m.
This is what we want to duplicate, particularly since returning the last 8 digits in to the same as modulo 108. Here is the VBA translation of Pow(b,e,m):
'pow(base,exponent,modulus): b^e mod m
'That works as long as (m-1)^2 fits into your integer type.
Dim a As Variant, x As Variant
If e = 0 Then
Pow = 1
Exit Function
End If
a = CDec(1)
x = CDec(b - m * Int(b / m)) 'b Mod m
While (e GT 1)
If e And 1 Then a = a * x - m * Int(a * x / m) 'If odd e then ax Mod m
x = x * x - M * Int(x * x / M) 'x^2 Mod m
e = BitShift(e, 1)
Wend
Pow = a * x - m * Int(a * x / m) 'ax Mod m
End Function
I used decimal variants, so this will work for m-1 up to the square root of ~7.92e29, or about ~8.9e14. Big enough. BitShift in this case is integer division by 32. Here are those functions:
'Right shift positive, left shift negative
If shift GT 0 Then
BitShift = shr(value, shift)
Else
BitShift = shl(value, -shift)
End If
End Function
Public Function shr(ByVal value As Long, ByVal shift As Byte) As Long
'http://www.excely.com/excel-vba/bit-shifting-function/
'Right shifting is equal to dividing Value by 2^Shift.
Dim i As Byte
shr = value
If shift GT 0 Then
shr = Int(shr / (2 ^ shift))
End If
End Function
Public Function shl(ByVal value As Long, ByVal shift As Byte) As Long
'http://www.excely.com/excel-vba/bit-shifting-function/
'Left shifting is equal to multiplying Value by 2^Shift. But to avoid an overflow error we use small trick:
shl = value
If shift GT 0 Then
Dim i As Byte
Dim M As Long
For i = 1 To shift
M = shl And &H40000000 ' save 30th bit
shl = (shl And &H3FFFFFFF) ' clear 30th and 31st bits
shl = shl * 2 ' multiply by 2
If M 0 Then
shl = shl Or &H80000000 ' set 31st bit
End If
Next i
End If
End Function
The usual angle brackets substitutions are above. Altogether then, this is the code for Problem 188:
Const a As Long = 1777
Dim i As Long
Dim Answer As Long, T As Single
T = Timer
Answer = 1
For i = 1855 To 1 Step -1
Answer = Pow(a, Answer, 10 ^ 8)
Next i
Debug.Print Answer; " Time:"; Timer - T
End Sub
Simple enough, but a lot of homework for this one. It put some more tools in the toolbox, and runs in less than 1/10 of a second.
..mrt
This is really neat stuff. I love reading about your Euler exploits… Keep up the great work!
I tried the functions myself, and I’m confused about two things:
1. What is the purpose of BitShifting in this problem?
2. If Pow(b,e,10^8) returns only the last 8 digits of the answer, and you need to recursively use this value, how do you achieve the correct answer when you never have all the digits?
Thanks,
Matthew
Michael,
Another good one. I’m going to have to look closer at this bit-shifting business - I simply used the binary representation of the exponent. I had some trouble fitting everything into the right variable, so I had to do some manipulating. On the bright side, I’ve used this modular exponentiation for several problems, so it’s a nice tool to have.
Matthew,
I believe (and Michael can correct me if I’m wrong) that:
1) The bit shifting arises from the fact that, for example, 7^6 = (7^2)^3. Thus if you break a number into its binary representation, you have a much more efficient way to exponentiate a number. I believe the bit-shifting is simply a more efficient means of multiplying / dividing by 2.
2) It turns out that if you take the modulus at every step and use that number in your calculations, then your final result will remain the same. Take for example 756^4 mod 1000. Doing out the multiplication gives a result of 296. Now if we do it stepwise, 756*756 = 571536. Mod it by 1000 and we get 536. 536*756 = 405216, mod by 1000 = 216. Finally, 216*756 = 163296, mod by 1000 gives 296. See http://en.wikipedia.org/wiki/Modular_exponentiation.
Josh
Hi Matthew -
Thank you. There is a textual error above that I’ll fix soon. It’s not integer division by 32, rather it’s by 2.
The purpose of bit-shifting is to get the problem into manageable chunks that fit into small variable types. It uses the identity that MOD(a*b,n) = MOD(MOD(a,n)*MOD(b,n),n) (in spreadsheet terms) repeatedly. I could have just said e = e\2.
The final code is not recursive. It’s a loop and relies on the fact that if you’re only interested in the last eight digits of the product, you only need to multiply the last eight digits of the multiplicands. Much harder problem if Euler wanted the left eight digits.
…mrt
Hi Josh -
Thanks. One neat thing I learned was
If n and 1 then … ‘ tests for odd n
and
If not n and 1 then … ‘ tests for even n
…mrt
Michael,
Well that is a neat trick. Of course, it looks like gibberish to me, but it apparently does work. It seemed to be about 30-40% faster for testing Longs than doing (n mod 2).
Josh
The second line of the definition is missing an operation. In your notation, it should read:
a^^(k+1) = a^^(a^^k).
otherwise, it is just the definition of exponentiation.
Hi Andrew -
Thanks. I had too few, you have too many. According to Wiki and Euler, it’s actually a^^(k+1) = a^(a^^k), which is how HyperExp() executes it.
Josh -
I’ve googled several of these logic operations. See
http://academic.evergreen.edu/projects/biophysics/technotes/program/logic.htm
and
http://www.petesqbsite.com/sections/express/issue12/TUTALGOS.TXT
…mrt
Michael,
Thanks for the links. I’ve seen various programs that make use of these logical operations and I always have trouble deciphering what they mean. This will be a good help.
Josh