Monday, June 25, 2007

Number Theory Hacking

Every now and again you I see a sum that I try and code up into a scripting language to see if my intuition about the problem is right.
So on page 237 of Pickover's "A passion for mathematics" is see an interesting problem.
Find x where y is an integer and x squared times 1597 plus one is an integer.
So I coded this up in Ruby as

puts "trying to find y=sqrt(1597xˆ2 + 1) problem"
y=0.1
x=1
while !y.integer?
y=Math.sqrt((1597*(x*x))+1)
x=x+1
end
puts y

The answer is a very big number. So big that using a static language is likely to lead to errors as saying int x or somesuch will lead to overflows. However ruby just keeps churning away without yet reaching an answer.
So my question is
1. How do you rapidly code up these number theory problems?
2. What language has the right mixture of rapid programming and speed of execution for doing so?

4 comments:

Christophe said...

Here is my Haskell solution (yes, a strongly typed language...however with infinite integers). Suffice to say that I used exactly your method, and it occurred to me that doing so is dangerous, cause you're using floating point arithmethic, which is unexact, during the calculation.

Anyways, your solution recoded in Haskell is:

head $ dropWhile (not . (\x -> let y = sqrt (1597*((fromIntegral x)**2+1)) in fromIntegral (floor y) == y)) [1..]

If I were to use

Mike said...

I agree that haskell makes a good language for cases like this. A good lisp implementation works too if you like a lot of (). Here's my implementation, which takes a slightly different approach but avoids floating point math. Its pretty fast too. It makes it through the first million elements of foos in just 31 seconds on my 3 GHz Pentium D with GHC 6.6. It shouldn't be too hard to change squares and foos to be calculated just with addition, which should give it a good speed boost.

Mike said...

Oops. FOrgot the code. Here it is. I tried out the addition optimization I mentioned, and indeed its gives a bit better than 2x improvement, now doing 71805 foos per second.


--squares = map (^2) [1..]
squares = scanl1 (+) [1,3..]
--foos = map (\x -> 1597*x*x + 1) [1..]
foos = tail $ scanl (+) 1 [1597*1, 1597*3..]

check_match (s:ss) (f:ff) = case s `compare` f of
LT -> check_match ss (f:ff)
EQ -> s
GT -> check_match (s:ss) ff

answer :: Integer
answer = check_match squares (take 2000000 foos)

main = print answer

Paddy3118 said...

If you use this Python library the code can be:

from stormer import Pell
print Pell(1597).next()

I cheated in that the Reddit discussion mentioned Pells equation and a quick google found the library.

- Paddy.