#! /usr/bin/env python3
#
def roots_rc ( n, x, fx, q ):

#*****************************************************************************80
#
## roots_rc() solves a system of nonlinear equations using reverse communication.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    28 November 2016
#
#  Author:
#
#    Original FORTRAN77 version by Gaston Gonnet.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Gaston Gonnet,
#    On the Structure of Zero Finders,
#    BIT Numerical Mathematics,
#    Volume 17, Number 2, June 1977, pages 170-183.
#
#  Input:
#
#    integer N, the number of equations.
#
#    real X(N).  Before the first call, the user should
#    set X to an initial guess or estimate for the root.  Thereafter, the input
#    value of X should be the output value of XNEW from the previous call.
#
#    real FX(N), the value of the function at XNEW.
#
#    real Q(2*N+2,N+2).  Before the first call
#    for a given problem, the user must set Q(2*N+1,1) to 0.0.
#
#  Output:
#
#    real XNEW(N), a new point at which a function
#    value is requested.
#
#    real FERR, the function error, that is, the sum of
#    the absolute values of the most recently computed function vector.
#
#    real Q(2*N+2,N+2).  
#
  import numpy as np

  ferr = np.sum ( abs ( fx[0:n] ) )
#
#  Initialization if Q(2*N+1,1) = 0.0.
#
  if ( q[2*n,0] == 0.0 ):

    for i in range ( 0, n ):
      for j in range ( 0, n + 1 ):
        q[i,j] = 0.0
        q[i+1,j] = 0.0
      q[i,i] = 100.0
      q[i+n,i] = 1.0

    q[2*n,0:n] = 1.0E+30
    q[2*n+1,0:n] = float ( n )

    for i in range ( 0, n ):
      q[i+n,n] = x[i]

    q[0:n,n] = fx[0:n]

    q[2*n,n] = ferr
    q[2*n+1,n] = 0.0
    damp = 0.99

  else:

    jsus = 0

    for i in range ( 1, n + 1 ):

      if ( 2 * n <= q[2*n+1,i] ):
        q[2*n,i] = 1.0E+30

      if ( q[2*n+1,jsus] <  ( ( n + 3 ) // 2 ) ):
        jsus = i
 
      if ( ( n + 3 ) // 2 <= q[2*n+1,i] and q[2*n,jsus] < q[2*n,i] ):
        jsus = i

    for i in range ( 0, n ):
      q[i+n,jsus] = x[i]
      q[i,jsus] = fx[i]

    q[2*n,jsus] = ferr
    q[2*n+1,jsus] = 0
    jsma = 1
    damp = 0.0

    for j in range ( 0, n + 1 ):
 
      if ( 1.0E+30 / 10.0 < q[2*n,j] ):
        damp = 0.99
 
      if ( q[2*n,j] < q[2*n,jsma] ):
        jsma = j

    if ( jsma != n ):
      for i in range ( 0, 2 * n + 2 ):
        t = q[i,jsma]
        q[i,jsma] = q[i,n]
        q[i,n] = t

  q[0:n,n+1] = q[0:n,n]
#
#  Call the linear equation solver, which should not destroy the matrix
#  in Q(1:N,1:N), and should overwrite the solution into Q(1:N,N+2).
#
# q[0:n,n+1] = np.linalg.solve ( q[0:n,0:n], q[0:n,n+1] )
  q[0:n,n+1], res, rank, s = np.linalg.lstsq ( q[0:n,0:n], q[0:n,n+1], rcond = None )

  sump = np.sum ( q[0:n,n+1] )

  if ( abs ( 1.0 - sump ) <= 1.0E-10 ):
    print ( '' )
    print ( 'ROOT - Fatal error!' )
    print ( '  SUMP almost exactly 1.' )
    print ( '  SUMP = %g' % ( sump ) )
    raise Exception ( 'ROOTS_RC - Fatal error!' )

  xnew = np.zeros ( n )

  for i in range ( 0, n ):
    xnew[i] = q[i+n,n]
    for j in range ( 0, n ):
      xnew[i] = xnew[i] - q[i+n,j] * q[j,n+1]
#
#  If system not complete, damp the solution.
#
    xnew[i] = xnew[i] / ( 1.0 - sump ) * ( 1.0 - damp ) + q[i+n,n] * damp

  for j in range ( 0, n + 1 ):
    q[2*n+1,j] = q[2*n+1,j] + 1.0

  return xnew, ferr, q

def roots_rc_test ( ):

#*****************************************************************************80
#
## roots_rc_test() tests roots_rc().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    28 November 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  n = 4
  it_max = 30
  x = np.zeros ( n )
  xnew = np.zeros ( n )
  fx = np.zeros ( n )

  print ( '' )
  print ( 'roots_rc_test():' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  roots_rc() seeks a solution of' )
  print ( '  the N-dimensional nonlinear system F(X) = 0.' )
  print ( '' )
  print ( '       FERR          X' )
  print ( '' )
#
#  Initialization.
#
  q = np.zeros ( [ 2 * n + 2, n + 2 ] )

  xnew[0] = 1.2
  xnew[1:n] = 1.0

  it = 0

  while ( True ):

    x = xnew.copy ( )

    fx[0] = 1.0 - x[0]
    for i in range ( 1, n ):
      fx[i] = 10.0 * ( x[i] - x[i-1] ** 2 )

    if ( it == 0 ):
      print ( '                ' ),
    else:
      print ( '  %14.6g' % ( ferr ) ),
 
    for i in range ( 0, n ):
      print ( '  %14.6g' % ( x[i] ) ),
    print ( '' )

    xnew, ferr, q = roots_rc ( n, x, fx, q )

    if ( ferr < 1.0E-07 ):
      print ( '' )
      print ( '  Sum of |f(x)| less than tolerance.' )
      break

    if ( it_max < it ):
      print ( '' )
      print ( '  Too many iterations!' )
      break
 
    it = it + 1
#
#  Terminate.
#
  print ( '' )
  print ( 'roots_rc_test():' )
  print ( '  Normal end of execution.' )
  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    06 April 2013
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return None

if ( __name__ == '__main__' ):
  timestamp( )
  roots_rc_test ( )
  timestamp ( )