(a+b)^2 = a || b

output.txt

-- pairs below 10^9
0,1
8,1
20,25
30,25
88,209
494,209
494,1729
6048,1729
744,1984
5288,1984
2450,2500
2550,2500
3008,14336
68320,14336
4938,17284
60494,17284
20408,122449
734694,122449
21948,126201
725650,126201
33058,148761
669420,148761
35010,152100
660790,152100
43470,165025
626480,165025
101558,217124
464194,217124
108878,221089
448944,221089
123448,227904
420744,227904
127194,229449
413908,229449
152344,237969
371718,237969
213018,248521
289940,248521
217930,248900
284270,248900
249500,250000
250500,250000
1975308,2469136
3086420,2469136
2428460,2499481
2572578,2499481
1860868,11780496
74578140,11780496
2752468,13838096
69571340,13838096
3714494,15558529
65168448,15558529
3909744,15863329
64363598,15863329
6126850,18625625
56621900,18625625
6376900,18875625
55871850,18875625
9538328,21345856
47769960,21345856
13223140,23140496
40495868,23140496
15119568,23764321
37351790,23764321
19681060,24682281
30954378,24682281
19825244,24700304
30774148,24700304
24995000,25000000
25005000,25000000
55022100,179545801
585886298,179545801
110888778,222110889
444889444,222110889
186997808,245434624
322132944,245434624

squarePairs.hs

import Math.NumberTheory.Powers.Squares   (exactSquareRoot)
import System.Environment                 (getArgs)
import Text.Read                          (readMaybe)

reduce :: Integral a => a -> a -> a
reduce b n = 25*10^(2*n) - (10^(n+1)-1)*b

roots :: Integral a => a -> [a]
roots n = case exactSquareRoot n of
            Nothing -> []
            Just t  -> [-t, t]

standardForm :: (Integral a, RealFrac b, Floating b) => a -> (a, b)
standardForm = properFraction . logBase 10 . fromIntegral

solve :: Integral a => a -> [(a, a)]
solve k = [(5*10^n - b + r,b) | b <- [1..k],
                                let (n,_) = standardForm b, --t < logBase 10 2.5,
                                r <- roots (reduce b n)]

display :: Show a => (a, a) -> String
display (x, y) = show x ++ "," ++ show y

main :: IO ()
main = getArgs >>= run

run :: [String] -> IO ()
run [outputFile, inputNum] = do
  upper <- maybe (putStrLn "defaulting to 1000000" >> return 1000000) return (readMaybe inputNum)
  let values = solve upper
  writeFile outputFile (unlines $ display <$> values)
  putStrLn $ "Found " ++ show (length values) ++ " pairs"
run _ = print "usage: ./squarePairs outputFile upperBound"