18

In ghci, using the arithmoi package:

Math.NumberTheory.Powers.General> :set +s
Math.NumberTheory.Powers.General> integerRoot 786 ((10^32)^786)
100000000000000000000000000000000
(0.04 secs, 227,064 bytes)
Math.NumberTheory.Powers.General> integerRoot 787 ((10^32)^787)

After five minutes, it still hasn't responded. Why is it taking so long?

(From some ad-hoc testing, it appears to be slow for all choices larger than 787 and fast for all choices smaller.)

Daniel Wagner
  • 145,880
  • 9
  • 220
  • 380
  • 3
    I get a segmentation fault when I run this with the latest arithmoi version and GHCi version 8.0.1 on Mac OS 10.12.5, with the `:set +s` option. Without `:set +s` I get "Bus error: 10". Seems like this `integerRoot` function makes memory act pretty strangely. – Alex Jul 12 '17 at 22:09

1 Answers1

23

arithmoi implements integerRoot by getting an initial approximate root and refining its guess with Newton’s method. For (1032)786, the second approximation gets a really good starting point:

> appKthRoot 786 ((10^32)^786)
100000000000000005366162204393472

For (1032)787, the second approximation gets a really bad starting point. Like, really bad.

> appKthRoot 787 ((10^32)^787)   
1797693134862315907729305190789024733617976978942306572734300811577326758055009
6313270847732240753602112011387987139335765878976881441662249284743063947412437
7767893424865485276302219601246094119453082952085005768838150682342462881473913
110540827237163350510684586298239947245938479716304835356329624224137216

It actually gets this approximation for everything starting there.

> length $ nub [appKthRoot x ((10^32)^x) | x <- [787..1000]]
1

Anyway, putting in the important parts of appKthRoot, we get:

> let h = 106; k = 786; n = (10^32)^k; !(I# s) = h * k - k in floor (scaleFloat (h - 1) (fromInteger (n `shiftRInteger` s) ** (1/fromIntegral k) :: Double))
100000000000000005366162204393472

> let h = 106; k = 787; n = (10^32)^k; !(I# s) = h * k - k in floor (scaleFloat (h - 1) (fromInteger (n `shiftRInteger` s) ** (1/fromIntegral k) :: Double))
179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137216

and taking a look at what’s going into scaleFloat:

> let h = 106; k = 786; n = (10^32)^k; !(I# s) = h * k - k in fromInteger (n `shiftRInteger` s) ** (1/fromIntegral k) :: Double
2.465190328815662

> let h = 106; k = 787; n = (10^32)^k; !(I# s) = h * k - k in fromInteger (n `shiftRInteger` s) ** (1/fromIntegral k) :: Double
Infinity

Yeah, that’d do it. (1032)786 ÷ 282530 &approx; 21023.1 fits in a double, but (1032)787 ÷ 282635 &approx; 21024.4 does not.

Ry-
  • 218,210
  • 55
  • 464
  • 476
  • 2
    That is frustrating indeed. The whole reason I have these big `Integer`s lying around in the first place (and part of the reason I'm using arithmoi) is because I'm working really hard to avoid `Double`... – Daniel Wagner Jul 12 '17 at 23:03
  • @DanielWagner: Yeah, I’m trying to figure out why they’re not just scaling by more than (h − 1)k. Might have a fix in a little while. – Ry- Jul 12 '17 at 23:07
  • @DanielWagner: For example, ``let !a@(I# a#) = 1024; h = 106; k = 787; n = (10^32)^k; !(I# s) = h * k - k in floor (scaleFloat 105 (fromInteger (n `shiftRInteger` (s +# a#)) ** (1/fromIntegral k) * 2**(fromIntegral a/k) :: Double))`` results in a reasonable approximation, and you could just keep scaling `a` up as appropriate. Might be worth raising an issue for. – Ry- Jul 12 '17 at 23:25