Bodigrim / arithmoi

Number theory: primes, arithmetic functions, modular computations, special sequences
http://hackage.haskell.org/package/arithmoi
MIT License
147 stars 40 forks source link

(WIP) Integer root fix #178

Closed b-mehta closed 4 years ago

b-mehta commented 4 years ago

Fix for issue #177.

Bodigrim commented 4 years ago

Huge thanks for investigation, now it makes sense.

It is probably worth to change otherwise case below in a similar fashion.

Also, when k >= 1024 even nshiftRInteger(h# *# k#) may appear larger than 2^1024 and thus be converted to Infinity :: Double. We can resort to a simpler approximation in this case:

appKthRoot k@(I# k#) n
  | k >= 1024 = 1 `shiftLInteger` (integerLog2# n `quotInt#` k# +# 1#)
  | otherwise = ...
b-mehta commented 4 years ago

Perhaps I'm missing something, but would

1 `shiftLInteger` (integerLog2# n `quotInt#` k#)

not be a better approximation?

b-mehta commented 4 years ago

Also, do you have suggestions for further test cases that could be added to verify these?

Bodigrim commented 4 years ago

I would probably just copy-paste Math.NumberTheory.Powers.GeneralTests.integerRootProperty to focus on huge inputs. Something along these lines:

integerRootHugeProperty :: Huge Natural -> Large Word -> Bool
integerRootHugeProperty (Huge n) (Large pow') = pow == 0 ||
  toInteger root ^ pow <= toInteger n && toInteger n < toInteger (root + 1) ^ pow
    where
      pow = pow' `mod` 2000
      root = integerRoot pow n
b-mehta commented 4 years ago

Huge inputs are already tested in testIntegral2Property, I think...

Bodigrim commented 4 years ago

Correct, but only for small powers, I believe.

b-mehta commented 4 years ago

Ah, I see!

Bodigrim commented 4 years ago

Nice, well done!