Euler Problem 188 asks:

The hyperexponentiation or tetration of a number a by a positive integer b, denoted by a^^b or

^{b}a, is recursively defined by:a^^1 = a,

a^^(k+1) = a(a^^k).Thus we have e.g. 3^^2 = 3

^{3}= 27, hence 3^^3 = 3^{27}= 7625597484987 and 3^^4 is roughly 10^{3.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 10^{8}. 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.

…mrtHi 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

not12^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