[Haskell-beginners] Informal proof of bijective mapping between Naturals and Natural pairs

[Brent Yorgey]

Hi Lyndon,

On Tue, Feb 01, 2011 at 12:43:14AM +0800, Lyndon Maydwell wrote:
> Hi Beginners.
>
> I was musing over weather Naturals were isomorphic to Natural pairs
> and wrote the following "proof".

Indeed they are.  This is a classic result first due to Cantor, who
came at it from the point of view of discovering that the natural
numbers have the same cardinality (size) as the positive rational
numbers (which can be thought of as pairs of naturals).

Your proof is great -- it has the advantage of working on rather
elementary principles.

However, there are many such bijections.  Here's another one, shorter
than yours although somewhat harder to understand just by looking at
the code:

  f :: (Integer, Integer) -> Integer
  f (x,y) = (r*r + r) `div` 2 + x
    where r = x + y

  f' :: Integer -> (Integer, Integer)
  f' n = (n - t, r - (n - t))
    where r = floor $ (sqrt(1 + 8*(fromIntegral n)) - 1) / 2
          t = (r*r + r) `div` 2

  *Main Test.QuickCheck> quickCheck (\(NonNegative x) (NonNegative y) -> f' (f (x,y)) == (x,y))
  +++ OK, passed 100 tests.
  *Main Test.QuickCheck> quickCheck (\(NonNegative x) -> f (f' x) == x)
  +++ OK, passed 100 tests.

This one looks a bit complicated but at heart it's rather simple: we
think of the pairs of naturals as residing in a grid, like so:

  .       .    	.
  :       :	:
  (0,2) (1,2) (2,2)  ...
  (0,1) (1,1) (2,1)  ...
  (0,0) (1,0) (2,0)  ...

and we imagine listing them by diagonals, like

  (0,0) (0,1) (1,0) (0,2) (1,1) (2,0) (0,3) ...

Now, the math to go back and forth between a pair and its index in
this listing by diagonals isn't quite as nice-looking as one might
hope (essentially we have to use the formula for triangular numbers
and solve a quadratic equation) but it's not too bad.  (Actually, what
I wrote above won't work once you get to integers big enough that the
sqrt starts losing precision, but you could write an accurate integer
square root operation to get around that.)

Some comments on your code style are interspersed below.

> import Control.Monad.Logic
> import Data.Maybe
> import Data.List
> import Data.Numbers.Primes
> import Test.QuickCheck
>
> -- This program demonstrates a mapping between the pairs of natural
> numbers and a subset of co-primes.
> -- The property should hold for all sized lists, not just pairs.
> Ordering is preserved.
> -- Question: Are we able to compress the range to create a bijection?
> -- Answer: Yes! We can use the breadth-wise indces of the products,
> rather than the products themselves.
>
> -- Pair to Number
>
> x_p = (evens !!)
> y_p = (odds  !!)
>
> xy_p x y = x_p x * y_p y
>
> xy_c x = ns_nc . xy_p x

The definition of xy_c strikes me as a gratuitous use of points-free
style just for its own sake.  It's hard to read because of the
asymmetric treatment of x and the (implicit) second argument.  I would
find it much clearer to just write

  xy_c x y = ns_nc (xy_p x y)

If you really want to do it points-free you could use a 'double
composition' operator,

  oo = (.) . (.)
  xy_c = ns_nc `oo` xy_p

but I don't recommend it.

Also, most of the difficulty I had reading your code was due to your
choice of function names, which were too telegraphic and didn't give
me any clues as to what they did.  Of course naming is rather an
idiomatic thing, but if you want others to read your code I suggest
trying to use names that are a bit more descriptive/suggestive.

>
> -- Number to Pair
>
> p_xy n = (x,y)
>   where
>     x = fromJust $ findIndex (`divides` n) evens
>     y = fromJust $ findIndex (`divides` n) odds
>
> c_xy = p_xy . nc_ns
>
> -- Sparse Number to Compact Number
>
> numbers = odds >>- \x -> evens >>- \y -> return (x*y) -- I don't
> really understand LogicT... :-(
>
> ns_nc n = fromJust $ findIndex (==n) numbers
>
> -- Compact Numner to Sparse Number
>
> nc_ns = (numbers !!)
>
> -- Helpers
>
> evens = map (primes !!) [0,2..]
> odds  = map (primes !!) [1,3..]

This strikes me as a rather inefficient way to define evens and odds.
It would be better to do something like

  deinterlace :: [a] -> ([a],[a])
  deinterlace (x1:x2:xs) = (x1:l1, x2:l2)
    where (l1,l2) = deinterlace xs

  (evens, odds) = deinterlace primes

>
> x `divides` y = y `mod` x == 0
>
> -- Id referencing properties
>
> prop_p1 = forAll (elements [(x,y) | x <- [1..6], y <- [1..7]]) f
>   where
>     f (x,y) = (x,y) == c_xy (xy_c x y)
>
> prop_p2 = forAll (elements [1..100]) f
>   where
>     f c = c == uncurry xy_c (c_xy c)