#! /usr/bin/env python3
#
def problem0 ( ):

#*****************************************************************************80
#
## problem0(): middle square random sequence
#
  import matplotlib.pyplot as plt

  print ( '' )
  print ( 'problem0():' )
  print ( '  Carry out the middle square random number sequence.' )
  print ( '' )

  d = 2

  for i in range ( 0, 11 ):

    if ( i == 0 ):

      s = 3647

    else:
#
#  Square.
#
      s = s * s
#
#  Drop last D digits.
#
      s = ( s // 10**d )
#
#  Drop first D digits.
#
      s = ( s % ( 10 ** ( 2 * d ) ) )

    print ( i, ': ', s )

  return

def problem1 ( ):

#*****************************************************************************80
#
## problem1(): Estimate pi in two lines.
# 
  import numpy as np

  print ( '' )
  print ( 'problem1():' )
  print ( '  Estimate pi in two lines of code.' )
  xy = np.random.random ( [ 100000, 2 ] )
  pi_est = 4 * sum ( xy[:,0]**2 + xy[:,1]**2 <= 1.0 ) / 100000
  print ( '  pi_estimate = ', pi_est )

  return
  
def problem2 ( ):

#*****************************************************************************80
#
## problem2(): calendar

  print ( '' )
  print ( 'problem2():' )
  print ( '  Given 0 <= d < 365, day of year' )
  print ( '  print month and day of month.' )

  d = 68
  print ( '  Use d=68 which should be March 10.' )

  month_name = [ 'January', 'February', 'March', 'April', 'May', 'June', \
            'July', 'August', 'September', 'October', 'November', 'December' ]

  month_length = [ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 ]

  for month in range ( 0, 12 ):

    if ( d < month_length[month] ):
      print ( '' )
      print ( '  The calendar date is', month_name[month], d + 1 )
      break

    d = d - month_length[month]

  return

def problem3 ( ):

#*****************************************************************************80
#
## problem3(): given card index, identify rank and suit.
#
  print ( '' )
  print ( 'problem3(): given card index c, identify rank and suit.' )

  c = 29
  print ( '  Card index c = ', c )

  suit = [ 'Hearts', 'Clubs', 'Diamonds', 'Spades' ]
  rank_index = ( c % 13 )
  rank = [ 'Ace', 'Two', 'Three', 'Four', 'Five', 'Six', 'Seven', \
           'Eight', 'Nine', 'Ten', 'Jack', 'Queen', 'King' ]
  k = c // 13

  print ( '  Card', c, 'is', rank[rank_index], 'of', suit[k] )
  
  return
  
def problem4 ( ):

#*****************************************************************************80
#
## problem4(): Newton-Pepys dice problem
#
  import numpy as np

  print ( '' )
  print ( 'problem4(): Newton-Pepys dice problem' )
  print ( '  Process 1: 6 dice are rolled, does at least 1 six show?' )
  print ( '  Process 2: 12 dice are rolled, does at least 2 sixes show?' )
  print ( '  Process 3: 18 dice are rolled, does at least 3 sixes show?' )

  n_test = 1000
  k1 = 0
  k2 = 0
  k3 = 0

  for i in range ( 0, n_test ):

    p1 = np.random.randint ( low = 1, high = 7, size = 6 )
    if ( 1 <= np.sum ( p1 == 6 ) ):
      k1 = k1 + 1
    p2 = np.random.randint ( low = 1, high = 7, size = 12 )
    if ( 2 <= np.sum ( p2 == 6 ) ):
      k2 = k2 + 1
    p3 = np.random.randint ( low = 1, high = 7, size = 18 )
    if ( 3 <= np.sum ( p3 == 6 ) ):
      k3 = k3 + 1

  print ( '' )
  print ( '  Process 1 probability estimate = ', k1 / n_test )
  print ( '  Process 2 probability estimate = ', k2 / n_test )
  print ( '  Process 3 probability estimate = ', k3 / n_test )

  return

def problem5 ( ):

#*****************************************************************************80
#
## problem5(): Probability of a straight in poker.
#
  import numpy as np

  print ( '' )
  print ( 'problem5():' )
  print ( '  Probability of a straight in poker.' )

  n_test = 1000000
  deck = np.arange ( 0, 52 )

  straight = 0

  for i in range ( 0, n_test ):
    hand = np.random.choice ( deck, 5, replace = False )
    hand = hand % 13
    hand.sort ()
    if ( is_straight ( hand ) ):
      straight = straight + 1
        
  print ( '' )
  print ( '  Estimated probability of a straight is ', straight / n_test )
  print ( '  Theoretical value is ', 10240 / 2598960 )
  return

def is_straight ( hand ):

  value = True
  for i in range ( 0, 4 ):
    if ( hand[i+1] != hand[i] + 1 ):
      value = False
      break

  return value
 
def problem6 ( ):

#*****************************************************************************80
#
## problem6(): Probability of a flush in poker.
#
  import numpy as np

  print ( '' )
  print ( 'problem6():' )
  print ( '  Probability of a flush in poker.' )

  n_test = 100000
  deck = np.arange ( 0, 52 )

  flush = 0

  for i in range ( 0, n_test ):
    c = np.random.choice ( deck, 5, replace = False )
    suit = c // 13
    d = np.unique ( suit )
    if ( len ( d ) == 1 ):
      flush = flush + 1

  print ( '' )
  print ( '  Estimated probability of a flush is ', flush / n_test )
  print ( '  Theoretical value is ', 5108 / 2598960 )

  return

def problem7 ( ):

#*****************************************************************************80
#
## problem7(): shared birthdays?
#
  import numpy as np

  print ( '' )
  print ( 'problem7(): shared birthdays among a group?' )
  n = 25
  b = np.random.choice ( range ( 0, 365 ), n, replace = True )
  u = np.unique ( b )
  print ( '' )
  print ( '  There were ', len ( u ), ' unique birthdays in group of ', n, 'people.' )

  return
  
def problem8 ( ):

#*****************************************************************************80
#
## problem8(): Stop children after first boy.
#
  import numpy as np

  print ( '' )
  print ( 'problem8(): Stop having children after first boy.' )

  n_test = 1000
  family_max = 0
  total_kids = 0
  total_boys = 0
  total_girls = 0

  for i in range ( 0, n_test ):

    family_girls = 0
    family_boys = 0

    while ( True ):
      baby = np.random.choice ( [ 'F', 'M' ], 1 )
      if ( baby == 'F' ):
        family_girls = family_girls + 1
      else:
        family_boys = family_boys + 1
        break

    family = family_girls + family_boys
    family_max = max ( family_max, family )
    total_boys = total_boys + family_boys
    total_girls = total_girls + family_girls
    total_kids = total_kids + family

  print ( '' )
  print ( '  Total kids born = ', total_kids )
  print ( '  Total girls born = ', total_girls )
  print ( '  Total boys born = ', total_boys )
  print ( '  Maximum family size was ', family_max )
  print ( '  Average family size was ', total_kids / n_test )

  return

def problem9 ( ):

#*****************************************************************************80
#
## problem9(): average length of monotone sequence using random values.
#
  import numpy as np

  print ( '' )
  print ( 'problem9():' )
  print ( '  How long does a random sequence stay monotone, on average?' )

  n_ave = 0
  n_test = 10000

  for test in range ( 0, n_test ):

    xmax = np.random.random ( )
    n = 1
    while ( True ):
      x = np.random.random ( )
      n = n + 1
      if ( xmax < x ):
        xmax = x
      else:
        n_ave = n_ave + n
        break

  n_ave = n_ave / n_test
  print ( '  Tested ', n_test, ' sequences' )
  print ( '  Average length of monotone sequence = ', n_ave )
  print ( '  Theoretical value is ', np.e )

  return

if ( __name__ == "__main__" ):
  problem0 ( )
  problem1 ( )
  problem2 ( )
  problem3 ( )
  problem4 ( )
  problem5 ( )
  problem6 ( )
  problem7 ( )
  problem8 ( )
  problem9 ( )