# Euler Problem 188

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:

Function HyperExp(a As Double, k As Double) As Double
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):

Public Function Pow(b As Variant, e As Variant, m As Variant) As Long
‘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:

Public Function BitShift(ByVal value As Long, ByVal shift As Integer) As Long
‘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 &amp;H40000000   ‘ save 30th bit
shl = (shl And &amp;H3FFFFFFF)   ‘ clear 30th and 31st bits
shl = shl * 2   ‘ multiply by 2
If M  0 Then
shl = shl Or &amp;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:

Sub Problem_188()

Const a As Long = 1777
Dim i As Long
Dim Answer As Long, T As Single

T = Timer

For i = 1855 To 1 Step -1
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

Posted in Uncategorized

## 18 thoughts on “Euler Problem 188”

1. Matthew Pfluger says:

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

2. JoshG says:

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

3. Michael says:

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

4. Michael says:

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

5. JoshG says:

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

6. Andrew says:

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.

7. Michael says:

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

and

http://www.petesqbsite.com/sections/express/issue12/TUTALGOS.TXT

…mrt

8. JoshG says:

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

9. Christian says:

I’m not sure, if modular theory is used correctly in this example… When writing

For i = 1855 To 1 Step -1
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…

10. Michael says:

Christian –

In spreadsheet terms, it uses MOD(a*b,n) = MOD(MOD(a,n)*MOD(b,n),n) which is valid.

…mrt

11. Christian says:

Let me explain, what irritated me. Have a look at the code

For i = 2 To 1 Step -1
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.

12. Michael says:

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

13. Christian says:

My point is that 12^12 is 12^^2

14. Michael says:

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

15. Michael says:

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:

a = 1777
For i = 1 to 1855-1
Next i

a = 12
For i = 1 to 2-1
Next i

I can’t satisfy myself with an articulation of the logic bug, so I’d invite your comments.

…mrt

16. Michael says:

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

17. Christian says:

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)

18. Michael says:

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

Posting code? Use <pre> tags for VBA and <code> tags for inline.