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 = e2.
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
I’m not sure, if modular theory is used correctly in this example… When writing
For i = 1855 To 1 Step -1
Answer = Pow(a, Answer, 10 ^ 8)
Next i
you seem to use the relation a^b = a^(b mod n) mod n; but this relation is false in general. Please correct me, if I’m wrong…
Christian –
In spreadsheet terms, it uses MOD(a*b,n) = MOD(MOD(a,n)*MOD(b,n),n) which is valid.
…mrt
Let me explain, what irritated me. Have a look at the code
For i = 2 To 1 Step -1
Answer = Pow(12, Answer, 10 )
Next i
After the first iteration, we have Answer=Pow(12, 1, 10 )=2, and so also after the second iteration, we have Answer=Pow(12, 2, 10 )=4.
On the other hand we have MOD(12^12,10)=6.
Christian –
It’s not 12^12. It’s hyper-exponentiation. It’s 12^^12, or in this case it’s not 1777^1885, it’s 1777^^1885.
See the Wiki article referenced above.
…mrt
My point is that 12^12 is 12^^2
Christian –
My apology, I saw that right after I posted. There’s something screwy using mod 10. If I use mod 100, it seems right.
I’ll think about it and post back. On the positive side, that code does check in the right answer ;-)
…mrt
Hi Christian –
There is a logic bug in the above code when m < b, and since nothing depends on i, there is no need to count backwards. I was influenced wrongly by having to go right-to-left. Thanks, I’d guess ;-) for making me relook at that code. I never considered that case.
For Euler’s case, this works:
Answer = a
For i = 1 to 1855-1
Answer = Pow(a, Answer, 10 ^ 8)
Next i
For your case:
Answer = a
For i = 1 to 2-1
Answer = Pow(a,answer,10)
Next i
I can’t satisfy myself with an articulation of the logic bug, so I’d invite your comments.
…mrt
The reason the original code provided the right answer is that the loop settled upon it in an amazingly few number of runs. I found this by dumping the answers into a spreadsheet.
…mrt
I’ve spent some time thinking about the reason for the correct behaviour of the code. I’m not sure if it’s a 100% correct, so I’d be glad to hear your opinion:
If a and n are coprime there is an integer e called the order of a modulo n with the property that x=y mod e implies a^x=a^y mod n. In case a=1777 and n=10^8, I think the order of a mod n is 2^4*5^7 or maybe 2^5*5^7. In any case it is a divisor of 10^8.
In the general situation, this must not necessarily hold (e.g. consider n=5 and a=2)
Hi Christian –
I’m an engineer (for 3 more months) and hobbyist programmer solving math problems, not a mathematician solving programming ones ;-) so I can’t really venture a creditable mathematical opinion. There are those who hang around DDoE who can (i.e fzz), I’m just not one of them. When you have the right answer, you can check into the Euler problem and read the insights of other solvers, some of whom are math-types. All that said, my knee quivers that if the general case is bad, the understanding is bad.
In the case of Euler #188, the loop can be too large or too small by a significant number relative to 1855, and still give the right answer. And while Euler is fond of problems that support a millisec-quick answer, I think in this case it’s happenstance. I lucked up with a bad algorithm giving a good result.
…mrt