Modular exponentiation with a noninteger base



#14

I'd like to compute "floor(3.8^98) mod 100" for the latest Project Euler problem (#356). Does someone know how to do this?

Fast modular exponentiation using the binary method doesn't seem to work for noninteger bases (I may well be overlooking the clever modification that will make it work), and I'm finding nothing that can do the job. I tried replacing the power function with exp(exponent * ln base) and a Taylor series for exp(), using modular aggregation. This converges a million times too slowly, and just can't be the way to do this.

WolframAlpha will gladly do the calculation (the result is "20"), but it fails for greater exponents.

For exponent 9876543 and modulus 100000000 it still produces output. For exponent 98765432, it appears to bring down a server (!); for exponent 987654321, it switches its input interpretation over to stocks...

Any ideas or keywords that would help to research this? Thanks!

EDIT: Here's the link to PE problem 356. Warning: enormous time sink. Largest root for n=2 is '1/3*(4+16/(37+3*i*sqrt(303))^(1/3)+(37+3*i*sqrt(303))^(1/3))'.


Edited: 30 Oct 2011, 10:20 p.m.


#15

Use the "Russian Peasant" method of exponentiation. Iterate to compute b^2, b^4, b^8, b^16, etc. Then b^98 = b^64 * b^32 * b^2. Even for the problem 356 exponent like 987654321, this takes only 29 squarings (log2 of 987654321), so no more than 58 multiplications in total for each exponential.


#16

Thanks, Eric. "Russian Peasant" == "binary method". So, that's what I tried.

#17

This by itself doesn't solve the problem -- even with this approach, even with a CAS, you need to figure out how to truncate the results so you still get the correct result and so the computer doesn't have to handle billion-digit arithmetic.

#18

I'm getting 60, not 20, for that problem, with both the HP50g and TI89.

Basic idea: You can calculate 3.8^98 exactly as (38/10)^98.

You can also get the approximate value of that.

And then convert that approximate value to an exact value (eg., 6.558e50 becomes 6 followed by 50 zeroes), and subtract.

Repeat until you get a value that's small enough (less than a billion, say) that you can read all the values down to the decimal point exactly.

And then you read the last two digits.

In fact, Wolfram Alpha seems to give the same thing (...160.982...)

So why is it giving 20?


#19

Utterly bizarre.

658858978275742440319523292513429252482438987872804536320
658858978275745389910397903226911721462945575230733045160

The top is what Wolfram Alpha gives for floor of 3.8^98.

The bottom is what Wolfram Alpha gives for 3.8^98 (truncated at the same number of digits)

This seems like a bug in Wolfram Alpha.

Edited: 31 Oct 2011, 12:01 a.m.


#20

Using bc I get:

scale=150
3.8^98
658858978275745389910397903226911721462945575230733045160.9820465028\
70158400628629120383359841551012084671432225137467649553611674541077\
38631648377416187904


- Pauli

#21

Might it be the floating point number (3.8) confusing things?

Calculate it as floor( (19/5)^98 ) mod 100 and it gets the correct answer.


- Pauli

#22

Thanks, Crawl and Pauli!

I never thought for a minute that WolframAlpha could be wrong. The query "floor((38/10)^98) mod 100" computes the correct value.

EDIT: Interestingly, while I now have floor-result harmony between WolframAlpha, bc (thanks, Pauli, for the steps!), and a BigFloat result in my calculator, they each (!) go their own way when it comes to 987 as a power. Shared digits in the beginning and divergence further on.

The 50g, following the steps by Crawl (thanks!), actually accepts to EVAL '(19/5)^987', but didn't finish after a long time, so won't be providing a result to compare. (The resulting number has 500+ digits.)

Edited: 31 Oct 2011, 9:15 a.m.


#23

For the larger exponent, you'll need to increase the scale parameter significantly in bc. 1000 or 1500 ought to be enough but err cautiously and go for e.g. 2500


- Pauli

#24

Ok. Notwithstanding the WolframAlpha bug, the problem remains how to calculate the modular result, without calculating the whole number.

When I hear "exponentiation" and "8 digits", I hear modpow with 10^8.

A 50g will easily tell you the result for 4^98 mod 100:

<< 100 MODSTO 4 98 POWMOD >>
If I try this with '19/5' instead of 4, "No solution in ring" is returned. (3.8 yields "Bad Argument Type".)

It returns a negative number (?) for "4^987654321 mod 100000000" (i.e. 100000000 MODSTO 4 987654321 POWMOD) but WolframAlpha and my other calc return the correct result, in no time.

If I do modpow(3.8, 98, 100), with a standard implementation for modpow

	// credit: book "Applied Cryptography" by Bruce Schneier
result = 1;
while (exponent > 0) {
if ((exponent & 1))
result = (result * base) % modulus;
exponent >>= 1;
base = (base * base) % modulus;
}
return result;
an incorrect result is returned. (Note, for floats, "%" becomes an "fmod" and that may be the catch here.)

This works correctly for integer bases.

It also works for a small power (such as 9), suggesting that precision issues could be the problem.

(123.456 mod 100 == 23.456000000000003 in IEEE 754.)

For solving this PE problem, this is likely a wrong attack angle altogether. Brute-force hardly suffices with these problems, and the roots will likely not do, expressed as reals. Purely symbolic solutions are rare in PE, but I'm not optimistic that the task is as straightforward as it may look at first glance.

At any rate, I'd still like to know how to compute the last few digits of a fraction, or real, raised to a large power.

To spell out the PE problem, the complete calculation becomes:

floor(1/3*(4+16/(37+3*i*sqrt(303))^(1/3)+(37+3*i*sqrt(303))^(1/3)) ^ 987654321) mod 100000000

That's for the second(-simplest) root. The last term in full:

floor(1/3 * (1073741824+1152921504606846976/(1237940039285380274899123819+
9*i*sqrt(12379400392853802748991240215))^(1/3)+(1237940039285380274899123819+
9*i*sqrt(12379400392853802748991240215))^(1/3)) ^ 987654321) mod 100000000

WolframAlpha will not compute the results for these.

I'd like to do this on my calculator...

[EDIT: On reflection, as a noninteger float is multiplied the digits are spreading in both direction. So, a large exponent would cause an enormous amount of spread, and all digits needed to be retained to not cause wildly different results in the output. So... inaccurate representation reasons aside, I think it's safe to say this cannot be solved with floats. A fraction would require huge numbers, too, for the result to be accurate. The solution will be a different method, and I guess properties of those roots may well come into play.

Too bad PE doesn't let you "give up" and look at the solution thread to learn from it.]


Edited: 31 Oct 2011, 8:48 p.m.

#25

I wonder if the multinomial expansion could help here. If the term is greater than 100000000*, it can be dismissed immediately. You still have to figure out how many terms to start with.


*-that is, the power of ten that multiplies the integers is greater than that.


#26

I don't know what that is, but will read up on it. Thanks, Crawl. (Learning new stuff sure is a fun part of PE.)


Possibly Related Threads…
Thread Author Replies Views Last Post
  Interesting Base Conversions - Porting a 1975 HP 25 Program to the HP 35S Eddie W. Shore 1 1,359 10-13-2013, 07:49 PM
Last Post: BruceH
  Wp34s base mode question Cristian Arezzini 7 2,003 05-19-2011, 06:46 AM
Last Post: Marcus von Cube, Germany
  WP34 base page view Cristian Arezzini 6 2,048 04-26-2011, 05:51 AM
Last Post: Paul Dale
  HP 35s. Base conversions, fast way Pablo P (Spain) 0 764 06-30-2010, 02:00 PM
Last Post: Pablo P (Spain)
  Can the HP 38G do Binary, Hexadecimal, or Octal base kmaciej 2 1,311 02-17-2009, 07:28 PM
Last Post: BruceH
  Base of Natural Logarithms Vincze 20 4,534 08-14-2007, 08:00 PM
Last Post: Miguel Toro
  dv6258SE base hot when using Scooby Doo 5 1,564 03-02-2007, 11:47 AM
Last Post: cfh
  Log Base 2 on HP15C John 5 1,884 08-09-2005, 03:19 PM
Last Post: John
  Anyone knows what happened with the hp knowledge base? (nt) Raul Lion 2 1,169 04-08-2004, 05:41 PM
Last Post: Raul Lion
  User base estimate Jose Ignacio Gutierrez 0 595 01-26-2003, 01:34 PM
Last Post: Jose Ignacio Gutierrez

Forum Jump: