#! /usr/bin/env python
#
def pascal2 ( n ):

#*****************************************************************************80
#
## PASCAL2 returns the PASCAL2 matrix.
#
#  Formula:
#
#    If ( I = 1 or J = 1 )
#      A(I,J) = 1
#    else
#      A(I,J) = A(I-1,J) + A(I,J-1)
#
#  Example:
#
#    N = 5
#
#    1 1  1  1  1
#    1 2  3  4  5
#    1 3  6 10 15
#    1 4 10 20 35
#    1 5 15 35 70
#
#  Properties:
#
#    A is a "chunk" of the Pascal binomial combinatorial triangle.
#
#    A is positive definite.
#
#    A is symmetric: A' = A.
#
#    Because A is symmetric, it is normal.
#
#    Because A is normal, it is diagonalizable.
#
#    A is integral, therefore det ( A ) is integral, and 
#    det ( A ) * inverse ( A ) is integral.
#
#    A is nonsingular.
#
#    det ( A ) = 1.
#
#    A is unimodular.
#
#    Eigenvalues of A occur in reciprocal pairs.
#
#    The condition number of A is approximately 16^N / ( N*PI ).
#
#    The elements of the inverse of A are integers.
#
#    A(I,J) = (I+J-2)% / ( (I-1)! * (J-1)! )
#
#    The Cholesky factor of A is a lower triangular matrix R,
#    such that A = R * R'.  The matrix R is a Pascal
#    matrix of the type generated by subroutine PASCAL.  In other
#    words, PASCAL2 = PASCAL * PASCAL'.
#
#    If the (N,N) entry of A is decreased by 1, the matrix is singular.
#
#    Gregory and Karney consider a generalization of this matrix as
#    their test matrix 3.7, in which every element is multiplied by a
#    nonzero constant K.  They point out that if K is the reciprocal of
#    an integer, then the inverse matrix has all integer entries.
#
#    The family of matrices is nested as a function of N.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 February 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Robert Brawer, Magnus Pirovino,
#    The linear algebra of the Pascal matrix,
#    Linear Algebra and Applications,
#    Volume 174, 1992, pages 13-23.
#
#    Robert Gregory, David Karney,
#    Example 3.7,
#    A Collection of Matrices for Testing Computational Algorithms,
#    Wiley, 1969, page 32, 
#    LC: QA263.G68.
#
#    Nicholas Higham,
#    Accuracy and Stability of Numerical Algorithms,
#    Society for Industrial and Applied Mathematics,
#    Philadelphia, PA, USA, 1996; section 26.4.
#
#    Sam Karlin,
#    Total Positivity, Volume 1,
#    Stanford University Press, 1968.
#
#    Morris Newman, John Todd,
#    The evaluation of matrix inversion programs,
#    Journal of the Society for Industrial and Applied Mathematics,
#    Volume 6, Number 4, pages 466-476, 1958.
#
#    Heinz Rutishauser,
#    On test matrices,
#    Programmation en Mathematiques Numeriques,
#    Centre National de la Recherche Scientifique,
#    1966, pages 349-365.
#
#    John Todd,
#    Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
#    Academic Press, 1977, page 172.
#
#    HW Turnbull,
#    The Theory of Determinants, Matrices, and Invariants,
#    Blackie, 1929.
#
#  Parameters:
#
#    Input, integer N, the order of A.
#
#    Output, real A(N,N), the matrix.
#
  import numpy as np

  a = np.zeros ( ( n, n ) )

  for i in range ( 0, n ):
    for j in range ( 0, n ):

      if ( i == 0 ):
        a[i,j] = 1.0
      elif ( j == 0 ):
        a[i,j] = 1.0
      else:
        a[i,j] = a[i,j-1] + a[i-1,j]

  return a

def pascal2_determinant ( n ):

#*****************************************************************************80
#
## PASCAL2_DETERMINANT computes the determinant of the PASCAL2 matrix.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 February 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of the matrix.
#
#    Output, real VALUE, the determinant.
#
  value = 1.0

  return value

def pascal2_determinant_test ( ):

#*****************************************************************************80
#
## PASCAL2_DETERMINANT_TEST tests PASCAL2_DETERMINANT.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 February 2015
#
#  Author:
#
#    John Burkardt
#
  import platform
  from pascal2 import pascal2
  from r8mat_print import r8mat_print

  print ( '' )
  print ( 'PASCAL2_DETERMINANT_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  PASCAL2_DETERMINANT computes the PASCAL2 determinant.' )

  m = 5
  n = m
 
  a = pascal1 ( n )

  r8mat_print ( m, n, a, '  PASCAL2 matrix:' )

  value = pascal2_determinant ( n )

  print ( '  Value =  %g' % ( value ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'PASCAL2_DETERMINANT_TEST' )
  print ( '  Normal end of execution.' )
  return

def pascal2_inverse ( n ):

#*****************************************************************************80
#
## PASCAL2_INVERSE returns the inverse of the PASCAL2 matrix.
#
#  Formula:
#
#    A(I,J) = sum ( max(I,J) <= K <= N )
#      (-1)^(J+I) * COMB(K-1,I-1) * COMB(K-1,J-1)
#
#  Example:
#
#    N = 5
#
#      5 -10  10  -5   1
#    -10  30 -35  19  -4
#     10 -35  46 -27   6
#     -5  19 -27  17  -4
#      1  -4   6  -4   1
#
#  Properties:
#
#    A is symmetric: A' = A.
#
#    Because A is symmetric, it is normal.
#
#    Because A is normal, it is diagonalizable.
#
#    A is integral, therefore det ( A ) is integral, and 
#    det ( A ) * inverse ( A ) is integral.
#
#    The first row sums to 1, the others to 0.
#
#    The first column sums to 1, the others to 0.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    16 March 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of A.
#
#    Output, real A(N,N), the matrix.
#
  import numpy as np
  from r8_choose import r8_choose
  from r8_mop import r8_mop

  a = np.zeros ( ( n, n ) )

  for i in range ( 0, n ):
    for j in range ( 0, n ):
      klo = max ( i + 1, j + 1 )
      for k in range ( klo, n + 1 ):
        a[i,j] = a[i,j] + r8_mop ( i + j ) * r8_choose ( k - 1, i ) \
          * r8_choose ( k - 1, j )

  return a

def pascal2_llt ( n ):

#*****************************************************************************80
#
## PASCAL2_LLT returns the Cholesky factor of the PASCAL2 matrix.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    16 March 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of A.
#
#    Output, real A(N,N), the matrix.
#
  from pascal1 import pascal1

  a = pascal1 ( n )

  return a

def pascal2_plu ( n ):

#*****************************************************************************80
#
## PASCAL2_PLU returns the PLU factors of the PASCAL2 matrix.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    24 March 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of the matrix.
#
#    Output, real P(N,N), L(N,N), U(N,N), the PLU factors.
#
  import numpy as np
  from pascal1 import pascal1

  p = np.zeros ( ( n, n ) )

  for j in range ( 0, n ):
    p[j,j] = 1.0

  l = pascal1 ( n )

  u = np.transpose ( l )

  return p, l, u

def pascal2_test ( ):

#*****************************************************************************80
#
## PASCAL2_TEST tests PASCAL2.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 February 2015
#
#  Author:
#
#    John Burkardt
#
  import platform
  from r8mat_print import r8mat_print

  print ( '' )
  print ( 'PASCAL2_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  PASCAL2 computes the PASCAL2 matrix.' )

  m = 5
  n = m

  a = pascal2 ( n )
 
  r8mat_print ( m, n, a, '  PASCAL2 matrix:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'PASCAL2_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  pascal2_test ( )
  timestamp ( )