Pascal to Sierpiński with help from Lucas

Tags: / math /

Finding Sierpiński triangle in Pascal's triangle with help of Lucas's theorem.

I am one of the TAs for the lab sessions of a UG-level course that our teacher is teaching this semester and during the first lab session, the teacher gave us a bonus task, since a lot of the students had finished the programs to be done early.

This blog post is a description of this task and a solution to it.


The problem

Print the Pascal's triangle upto a given height, but with the numbers in it made modulo 2 of its actual value.

ie, if the Pascal triangle is like this:

          1
        1   1
      1   2   1
    1   3   3   1
  1   4   6   4   1
1   5  10  10   5   1

the output should be the same triangle but with their values modulo 2. Like:

        1
      1   1
    1   0   1
  1   1   1   1
1   0   0   0   1

where the even numbers get replaced with 0 and odd numbers with 1.

Pascal's triangle

(I had mostly forgotten what Pascal's triangle was. So a refresher.)

As was shown above, Pascal's triangle looks like:


                01
              01  01
            01  02  01
          01  03  03  01
        01  04  06  04  01
      01  05  10  10  05  01
    01  06  15  20  15  06  01
  01  07  21  35  35  21  07  01
01  08  28  56  70  56  28  08  01

where the value of a position in a row is the sum of the nearest two values in the row directly above that position.

Several people from different parts of the world had studied this triangle long before the time of Blaise Pascal (1623-1662 ¹⁹).

But Pascal's Traité du triangle arithmétique (published posthumously in 1664), which dealt with the properties of this triangle quite comprehensively, was sort of popular enough for Pascal's name to get attached to the triangle.

Since it was known well before Pascal's time it is known by many names. Among them are:

Was a bit surprised to know that Omar Khayyam, to whom The Rubáiyát of Omar Khayyám is attributed to, was more famous as a mathematician and had studied (what later became known as) the Pascal's triangle among other things.

(I once tried to read the Rubáiyát, but couldn't really understand much. 😬)

Pascal's triangle attracts the attention of mathematicians because it has several interesting properties.

One of them is that its values correspond to coefficients of expansion of Binomial theorem.

            01
          01  01
        01  02  01
      01  03  03  01
    01  04  06  04  01
  01  05  10  10  05  01
01  06  15  20  15  06  01
           ...
           ...

Recall that binomial theorem is like:

            n  
(x + y)ⁿ =  Σ  nCi(xⁿ⁻ⁱ)(yⁱ)
           i=0

The elements of Pascal's triangle are same as the coefficients from binomial theorem.

Row i of the Pascal triangle consists of the coefficients of the expansion of (x+y)ⁿ.
(where i starts from 0)

                  0C0                      | Coeffs of (x + y)⁰
               1C0   1C1                   | Coeffs of (x + y)¹
            2C0   2C1   2C2                | Coeffs of (x + y)²
         3C0   3C1   3C2   3C3             | Coeffs of (x + y)³
      4C0   4C1   4C2   4C3   4C4          | Coeffs of (x + y)⁴
   5C0   5C1   5C2   5C3   5C4   5C5       | Coeffs of (x + y)⁵
6C0   6C1   6C2   6C3   6C4   6C5   6C6    | Coeffs of (x + y)⁶
                 ...
                 ...

where nCk is 'n choose k' (ie, ⁿCₖ).

A solution

For our problem, we can use this property of the Pascal triangle elements being coefficients of binomial expansion.

All we got to do is to find the appropriate aCb values, divide them by 2 and take the remainder.

           a!
aCb = ─────────── 
      (a-b)! . b!

ie, find these coefficients and do modulo 2.

Still.. that sounds like some unnecessary calculations are involved, right?

I mean, we figure out the full value of the aCb-s just to see if it's odd or even. We calculate the whole thing and just truncate it by doing a modulo 2.

If we write a program to print the first n rows of Pascal's triangle by actually calculating the whole aCb values and then doing modulo 2, it's going to be a bit slow.

Would've been nice if we could just cut to the case and get to the mod 2 value 'more directly'.

Turns out there are such way(s). One way is by taking advantage of Lucas's theorem.

Lucas's theorem

This theorem is named after the French mathematician Édouard Lucas, who came up with it around 1878.

Lucas's theorem gives us some information about the remainder resulting from the division of a aCb by a prime number p when a and b are expressed in base-p expansion.

a and b in base-p expansion meaning something like:

a = aₖpᵏ + aₖ₋₁pᵏ⁻¹ + .... + a₀
b = bₖpᵏ + bₖ₋₁pᵏ⁻¹ + .... + b₀

then, Lucas's theorem tells us that the following is relation holds:


        ⎡   k        ⎤ 
        ⎢  ┬─┬       ⎥
aCb  ≡  ⎢  │ │ aᵢCbᵢ ⎥  (mod p)
        ⎣  i=0       ⎦  

where aCb is 'a choose b'.

This means that the ∏ aᵢCbᵢ would be the reminder if we divide aCb with p.

ie,


             k                   
            ┬─┬                 
aCb % p  =  │ │ aᵢCbᵢ 
            i=0                   

((mod p) and mod p are not the same. See this.)

A faster solution

If we use binary representations of a and b, we can use the equation from Lucas's theorem with p = 2 to get what we want.

ie, figure out whether a given aCb value is odd or even.


             k                   
            ┬─┬                 
aCb % 2  =  │ │ aᵢCbᵢ 
            i=0                   

This means significantly lesser number of operations to our earlier solution as we knock off a few big factorials in favour a few small multiplications.

For example, consider 10C4, which is 210. An even number.

Let's try using Lucas theorem to figure that out.

p = 2
a = 10₁₀ = 1010₂
b = 4₁₀ = 0100₂


By Lucas theorem, we got


                  3        
                 ┬─┬       
    10C4 % 2  =  │ │ aᵢCbᵢ
                 i=0        


=>  10C4 % 2 =  a₀Cb₀ *  a₁Cb₁ * a₂Cb₂ * a₃Cb₃


We know values for aᵢ and bᵢ as:

|    | 0 | 1 | 2 | 3 |
|----+---+---+---+---|
| aᵢ | 0 | 1 | 0 | 1 |
| bᵢ | 0 | 0 | 1 | 0 |

Substituting,


=>  10C4 % 2  =  0C0 *  1C0 * 0C1 * 1C0


We know the following values:
 - 0C0 = 1
 - 0C1 = 0
 - 1C0 = 1
 - 1C1 = 1


=>  10C4 % 2  =  1 *  1 * 0 * 1
=>  10C4 % 2  =  0      

ie, 10C4 divided by 2 would give remainder 0.
That means 10C4 results in an even value.

Which is correct, since 210 is even.

(A couple more examples is included in the addendum section of this blog post.)

In fact, we could git rid of a few more calculations since the Pascal's triangle is symmetrical. ie, values appearing in the right half are the same values as appearing in the left half, but in the reverse order.

An implementation

I made a Haskell program that prints the Pascal's triangle modulo 2 using facts that we can tell from Lucas's theorem:

module Main where

-- Row indexing starts from zero.
-- Rows on top are of smaller index.
--
-- Row 0  |        1
-- Row 1  |       1 1
-- Row 2  |      1 0 1
-- Row 3  |     1 1 1 1
-- Row 4  |    1 0 0 0 1


-- | Find binary representation of an integer
bin
  :: Int    -- ^ number n
  -> [Int]  -- ^ digits of n to the base n. LSB first.
bin n
  | n < 2 = [n]
  | otherwise = rem:bin nxt
     where
      rem = mod n 2
      nxt = quot n 2


-- | Find value for aCb where a and b ∈ {0, 1}
aCb01
  :: (Int, Int)  -- ^ a and b values as tuple
  -> Int         -- ^ aCb value
aCb01 (0, 1) = 0
aCb01 _ = 1


-- | Calculate aCb value
aCb
  :: Int   -- ^ a
  -> Int   -- ^ b
  -> Int   -- ^ aCb
aCb a b = product $ map aCb01 $ zip apad bpad
 where
  abin = bin a
  bbin = bin b
  alen = length abin
  blen = length bbin
  maxlen = max alen blen
  apad = abin ++ replicate (maxlen-alen) 0
  bpad = bbin ++ replicate (maxlen-blen) 0


-- | Find values for one row
getRow
  :: Int    -- ^ row index
  -> [Int]  -- ^ values of the row
getRow idx = [aCb idx r | r <- [0..idx]]


-- | Find values of all rows
getRows
  :: Int      -- ^ index of largest row
  -> [[Int]]  -- ^ values of rows. Values of each row is in a separate list
getRows rowMaxIdx = [getRow i | i <- [0..rowMaxIdx ]]


-- | Find values of a row and format it as a string
fmtRow
  :: Int     -- ^ index of largest row
  -> Int     -- ^ row index
  -> String  -- ^ row data as formatted string
fmtRow maxIdx idx = begg ++ bodyhead ++ bodytail
 where
  begg = replicate (maxIdx-idx) ' '
  digitstr = map show $ getRow idx
  bodyhead = head digitstr
  bodytail = foldl (++) "" $ map ((++) " ") $ tail digitstr


-- | Find values of the entire Pascal triangle and format it as a string
pascal
  :: Int     -- ^ maximum row
  -> String  -- ^ formatted, read-to-print, string
pascal maxIdx = foldl (++) ""  $ map ((++) "\n") strLines
 where
  strLines = [fmtRow maxIdx idx | idx <- [0..maxIdx]]


main = do
  putStrLn "Enter number of rows:"
  rowCountStr <- getLine
  let rowCount = read rowCountStr :: Int
  let ll = getRows rowCount
  putStrLn $ pascal $ rowCount - 1

Sample output:

Enter number of rows:
10

         1
        1 1
       1 0 1
      1 1 1 1
     1 0 0 0 1
    1 1 0 0 1 1
   1 0 1 0 1 0 1
  1 1 1 1 1 1 1 1
 1 0 0 0 0 0 0 0 1
1 1 0 0 0 0 0 0 1 1

Okay, cool. Now that we've got an answer to our problem along with a program to try it out, let's take a look at another 'triangle' which has some connection to Pascal's triangle: Sierpiński triangle.

Sierpiński triangle

Sierpiński (pronounced like 'shee-yer-piyen-ski') triangle is a self-replicating pattern (a fractal).

Named after Wacław Sierpiński, a mathematician from Poland, it's also known as Sierpiński gasket or Sierpiński sieve.

It's construction goes like this:

  1. Start with an equilateral triangle.
  2. At each step do the following for each constituent non-empty triangle:
    1. Form four smaller congruent equilateral triangles by joining the midpoints of each side.
    2. Remove the central triangle.

Here's a representation of the construction of Sierpiński triangle:
(Parts of the triangle which are not 'filled' with 1-s are the 'removed'/empty portions.)


            +                             +                             +            
           /1\                           /1\                           /1\           
          /111\                         /111\                         /111\          
         /11111\                       /11111\                       +-----+         
        /1111111\                     /1111111\                     /1\   /1\        
       /111111111\                   /111111111\                   /111\ /111\       
      /11111111111\                 +-----------+                 +-----+-----+      
     /1111111111111\               /1\         /1\               /1\         /1\     
    /111111111111111\             /111\       /111\             /111\       /111\    
   /11111111111111111\           /11111\     /11111\           +-----+     +-----+   
  /1111111111111111111\         /1111111\   /1111111\         /1\   /1\   /1\   /1\  
 /111111111111111111111\       /111111111\ /111111111\       /111\ /111\ /111\ /111\ 
+-----------------------+     +-----------+-----------+     +-----+-----+-----+-----+

        Level 1                         Level 2                      Level 3

As you might have already noticed, the number of non-empty triangles in the pattern is where l is the level.

Okay, so what's the connection with Sierpiński triangle and Pascal's triangle?

Let's try running the haskell program to print upto a not-so-small number of rows. Say 32 rows.

Enter number of rows:
32

                                 1
                                1 1
                               1 0 1
                              1 1 1 1
                             1 0 0 0 1
                            1 1 0 0 1 1
                           1 0 1 0 1 0 1
                          1 1 1 1 1 1 1 1
                         1 0 0 0 0 0 0 0 1
                        1 1 0 0 0 0 0 0 1 1
                       1 0 1 0 0 0 0 0 1 0 1
                      1 1 1 1 0 0 0 0 1 1 1 1
                     1 0 0 0 1 0 0 0 1 0 0 0 1
                    1 1 0 0 1 1 0 0 1 1 0 0 1 1
                   1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
                  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
                 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
                1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1
               1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1
              1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1
             1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1
            1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 1
           1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1
          1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
         1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1
        1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1
       1 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 1
      1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1
     1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1
    1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1
   1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Notice anything? Yes? No?

(I couldn't spot it when I first saw it.)

Try filling in the 1-s and removing the 0-s.

                               █
                              ███
                             █  ██
                            ███████
                           █      ██
                          ███    ████
                         █  ██  ██  ██
                        ███████████████
                       █              ██
                      ███            ████
                     █  ██          ██  ██
                    ███████        ████████
                   █      ██      ██      ██
                  ███    ████    ████    ████
                 █  ██  ██  ██  ██  ██  ██  ██
                ███████████████████████████████
               █                              ██
              ███                            ████
             █  ██                          ██  ██
            ███████                        ████████
           █      ██                      ██      ██
          ███    ████                    ████    ████
         █  ██  ██  ██                  ██  ██  ██  ██
        ███████████████                ████████████████
       █              ██              ██              ██
      ███            ████            ████            ████
     █  ██          ██  ██          ██  ██          ██  ██
    ███████        ████████        ████████        ████████
   █      ██      ██      ██      ██      ██      ██      ██
  ███    ████    ████    ████    ████    ████    ████    ████
 █  ██  ██  ██  ██  ██  ██  ██  ██  ██  ██  ██  ██  ██  ██  ██
███████████████████████████████████████████████████████████████

And we have something that looks pretty much like the Sierpiński triangle. 😃

From Wikipedia:

If one takes Pascal's triangle with 2ⁿ rows and colours the even numbers white, and the odd numbers black, the result is an approximation to the Sierpinski triangle. More precisely, the limit as n approaches infinity of this parity-colored 2ⁿ-row Pascal's triangle is the Sierpinski triangle.

That's why I went with 32 rows. Because it's a power of 2.

Acknowledgements: Thanks to Piyush sir, the teacher who told us about this problem and the way of using Lucas's theorem to solve it.

References

Addendum

Couple more examples using Lucas theorem

Suppose we need to find if 25C13 is odd or even.

p = 2
a = 25₁₀ = 11001₂
b = 13₁₀ = 01101₂

Lucas theorem gives us:


            ⎡   4        ⎤ 
            ⎢  ┬─┬       ⎥
   25C13 ≡  ⎢  │ │ aᵢCbᵢ ⎥  (mod 2)
            ⎣  i=0       ⎦  


            ⎡                                        ⎤ 
=> 25C13 ≡  ⎢ a₀Cb₀ *  a₁Cb₁ * a₂Cb₂ * a₃Cb₃ * a₄Cb₄ ⎥  (mod 2)
            ⎣                                        ⎦  


            ⎡                              ⎤ 
=> 25C13 ≡  ⎢ 1C1 *  0C0 * 0C1 * 1C1 * 1C0 ⎥  (mod 2)
            ⎣                              ⎦  


            ⎡                    ⎤ 
=> 25C13 ≡  ⎢ 1 *  1 * 0 * 1 * 1 ⎥  (mod 2)
            ⎣                    ⎦  


=> 25C13 ≡  0 (mod 2)

So, 25C13 is even.

Which is correct, because its value is 5200300.

Just for the fun of it, let's try one more number: 123C56

p = 2
a = 123₁₀ = 1111011₂
b = 56₁₀ = 0111000₂

Lucas theorem gives us:


            ⎡   6        ⎤ 
            ⎢  ┬─┬       ⎥
   25C13 ≡  ⎢  │ │ aᵢCbᵢ ⎥  (mod 2)
            ⎣  i=0       ⎦  


            ⎡                                 ⎤ 
=> 25C13 ≡  ⎢ a₀Cb₀ * a₁Cb₁ * a₂Cb₂ * a₃Cb₃ * ⎥  (mod 2)
            ⎢ a₄Cb₄ * a₅Cb₅ * a₆Cb₆           ⎥
            ⎣                                 ⎦  


            ⎡                                         ⎤ 
=> 25C13 ≡  ⎢ 1C0 * 1C0 * 0C0 * 1C1 * 1C1 * 1C1 * 1C0 ⎥  (mod 2)
            ⎣                                         ⎦  


            ⎡                           ⎤ 
=> 25C13 ≡  ⎢ 1 * 1 * 1 * 1 * 1 * 1 * 1 ⎥  (mod 2)
            ⎣                           ⎦  


=> 25C13 ≡  1 (mod 2)

ie, 25C13 is an odd value.

That's so because 25C13 is 468410292423894497225305443312892167.