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:

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 &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:

Sub 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

Posted in Uncategorized

18 thoughts on “Euler Problem 188

  1. 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. 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. 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. 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. 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. 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,
    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

  8. 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…

  9. Christian –

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

    …mrt

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

    Answer=1
    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.

  11. 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

  12. 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

  13. 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
    Answer = a
    For i = 1 to 1855-1
       Answer = Pow(a, Answer, 10 ^ 8)
    Next i

    For your case:

    a = 12
    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

  14. 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

  15. 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)

  16. 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.

Leave a Reply

Your email address will not be published. Required fields are marked *