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

#*****************************************************************************80
#
## babylonian_test() tests babylonian().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    21 June 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'babylonian_test():' )
  print ( '  Python version: ' + platform.python_version ( ) )
  print ( '  babylonian() demonstrates some Babylonian arithmetic.' )

  babylonian_sqrt_test ( )
  digits_decimal_to_r8_test ( )
  digits_sexagesimal_print_test ( )
  digits_sexagesimal_to_r8_test ( )
  r8_to_digits_decimal_test ( )
  r8_to_digits_sexagesimal_test ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'babylonian_test():' )
  print ( '  Normal end of execution.' )

  return

def babylonian_sqrt ( x, tol = 0.000001 ):

#*****************************************************************************80
#
## babylonian_sqrt() estimates sqrt(x) using a Babylonian method.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    21 June 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real X: the number whose square root is desired.
#
#    real TOL: a tolerance for the error in |X-R^2|
#
#  Output:
#
#    real R: the estimated square root.
#
#    integer N: the number of iterations used.
#
  n = 0

  while ( True ):

    if ( n == 0 ):
      r = x
    else:
      r = ( r + x / r ) / 2.0

    e = abs ( x - r * r )
    print ( n, r, x / r, e )
    if ( e < tol ):
      break
    n = n + 1

  return r, n
  
def babylonian_sqrt_test ( ):

#*****************************************************************************80
#
## babylonian_sqrt_test() tests babylonian_sqrt().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    21 June 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
 
  print ( '' )
  print ( 'babylonian_sqrt_test():' )
  print ( '  babylonian_sqrt() estimates square roots using a' )
  print ( '  method known to the Babylonians.' )

  x = 2.0
  tol = 0.000001
  r, n = babylonian_sqrt ( x, tol )
  print ( '  Seek sqrt ( ', x, ' )' )
  print ( '  Estimate =  ', r )
  print ( '  sqrt(r) =   ', np.sqrt ( x ) )
  print ( '  Err =       ', r - np.sqrt ( x ) )
  print ( '  Tolerance = ', tol )
  print ( '  Iterations = ', n )

  return

def digits_decimal_to_r8 ( d, p ):

#*****************************************************************************80
#
## digits_decimal_to_r8() evalulates an N digit decimal digit number.
#
#  Discussion:
#
#    R = d(0) * 10^p + d(1) * 10^(p-1) + ... + d(n-1) * 10^(p+1-n).
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 June 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer D(N): the decimal digits.
#
#    integer P: the power of 10 multiplying the first digit.
#
#  Output:
#
#    real R: the value represented.
#
  power = 10.0 ** p
  n = len ( d )
  r = 0.0
  for i in range ( 0, n ):
    r = r + d[i] * power
    power = power / 10.0

  return r

def digits_decimal_to_r8_test ( ):

#*****************************************************************************80
#
## digits_decimal_to_r8_test() tests digits_decimal_to_r8().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 June 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'digits_decimal_to_r8_test():' )
  print ( '  digits_decimal_to_r8() converts base 10 digits of a' )
  print ( '  number into a real value;' )

  for test in range ( 0, 3 ):

    if ( test == 0 ):
      d = np.array ( [ 1, 2, 0, 3, 4 ] )
      p = -2
    elif ( test == 1 ):
      d = np.array ( [ 9, 8, 7, 6, 5, 4, 3 ] )
      p = 3
    elif ( test == 2 ):
      d = np.array ( [ 3, 1, 4, 1, 5, 9, 2, 6, 5 ] )
      p = 0

    r8 = digits_decimal_to_r8 ( d, p )

    print ( '' )
    print ( '  Given decimal digits:' )
    print ( d )
    print ( '  and "decimal place" p = ', p )
    print ( '  corresponding value is r8 = ', r8 )

  return

def digits_sexagesimal_print ( d, p ):

#*****************************************************************************80
#
## digits_sexagesimal_print() prints a sexagesimal value.
#
#  Discussion:
#
#    R = d(0) * 60^p + d(1) * 60^(p-1) + ... + d(n-1) * 60^(p+1-n)
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    21 June 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer D(N): the digits
#
#    integer P: the power of 60 multiplying the first digit.
#
  n = len ( d )
  p1 = max ( p, 0 )
  p2 = min ( p + 1 - n, 0 )

  for p3 in range ( p1, p2 - 1, -1 ):

    if ( p + 1 - n <= p3 and p3 <= p ):
      i = p - p3 
      if ( p3 == 0 ):
        print ( '%02d.' % ( d[i] ), end = '' )
      elif ( p3 == p2 ):
        print ( '%02d' % ( d[i] ), end = '' )
      else:
        print ( '%02d,' % ( d[i] ), end = '' )
    else:
      if ( p3 == 0 ):
        print ( '00.', end = '' )
      else:
        print ( '00,', end = '' )

  print ( '' )

  return

def digits_sexagesimal_print_test ( ):

#*****************************************************************************80
#
## digits_sexagesimal_print_test() tests digits_sexagesimal_print().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    21 June 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'digits_sexagesimal_print_test():' )
  print ( '  digits_sexagesimal_print() prints a real number ' )
  print ( '  as a sexagesimal quantity;' )

  d = np.array ( [ 11, 12, 0, 14, 15 ] )
  print ( '  Sexagesimal digit vector:' )
  print ( d );

  for p in range ( -5, 6 ):

    print ( '  p = %2d:  ' % ( p ), end = '' )
    digits_sexagesimal_print ( d, p )

  return

def digits_sexagesimal_to_r8 ( d, p ):

#*****************************************************************************80
#
## digits_sexagesimal_to_r8() evalulates an N digit sexagesimal digit number.
#
#  Discussion:
#
#    R = d(0) * 60^p + d(1) * 60^(p-1) + ... + d(n-1) * 60^(p+1-n).
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    21 June 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer D(N): the sexagesimal digits.
#
#    integer P: the power of 60 multiplying the first digit.
#
#  Output:
#
#    real R: the value represented.
#
  power = 60.0 ** p
  r = 0.0
  for c in d:
    r = r + c * power
    power = power / 60.0

  return r

def digits_sexagesimal_to_r8_test ( ):

#*****************************************************************************80
#
## digits_sexagesimal_to_r8_test() tests digits_sexagesimal_to_r8().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    21 June 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'digits_sexagesimal_to_r8_test():' )
  print ( '  digits_sexagesimal_to_r8() converts base 60 digits of a' )
  print ( '  number into a real value;' )

  for test in range ( 0, 3 ):

    if ( test == 0 ):
      d = np.array ( [ 4, 37, 46, 40 ] )
      p = 3
    elif ( test == 1 ):
      d = np.array ( [ 3, 8, 29, 44, 0, 47, 25 ] )
      p = 0
    elif ( test == 2 ):
      d = np.array ( [ 1, 24, 51, 10, 7, 46, 6, 4 ] )
      p = 0

    r8 = digits_sexagesimal_to_r8 ( d, p )

    print ( '' )
    print ( '  Given sexagesimal digits:' )
    print ( d )
    print ( '  and sexagesimal "decimal place" p = ', p )
    print ( '  corresponding value is r8 = ', r8 )

  return

def r8_to_digits_decimal ( r, n ):

#*****************************************************************************80
#
## r8_to_digits_decimal() computes N decimal digits representing R.
#
#  Discussion:
#
#    R = d(0) * 10^p + d(1) * 10^(p-1) + ... + d(n-1) * 10^(p+1-n) approximately
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 June 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real R, the number to be represented.
#
#    integer N: the number of digits to compute.
#
#  Output:
#
#    integer D(N): the decimal digits.
#
#    integer P: the power of 10 multiplying the first digit.
#
  import numpy as np

  d = np.zeros ( n, dtype = np.int )
  p = 0

  if ( r == 0.0 ):
    return d, p
#
#  Try to deal with negative input.
#
  if ( r < 0.0 ):
    sgn = -1
    r = - r
  else:
    sgn = +1
#
#  Determine P such that 10^P <= R < 10^(P+1)
#
  rcopy = r
  base = 10.0
  while ( base <= rcopy ):
    rcopy = rcopy / base
    p = p + 1

  while ( rcopy < 1.0 ):
    rcopy = rcopy * base
    p = p - 1
#
#  Now read off N digits.
#
  rcopy = r
  div = base ** p
  for i in range ( 0, n ):
    d[i] = int ( np.floor ( rcopy / div ) )
    rcopy = ( rcopy % div )
    div = div / base
#
#  Multiply by sgn.
#
  d = sgn * d

  return d, p

def r8_to_digits_decimal_test ( ):

#*****************************************************************************80
#
## r8_to_digits_decimal_test() tests r8_to_digits_decimal().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 June 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'r8_to_digits_decimal_test():' )
  print ( '  r8_to_digits_decimal() converts a real number' )
  print ( '  to its representation in base 10;' )

  hundredth = 1.0 / 100.0
  sqrt2 = np.sqrt ( 2.0 )
  million = 1000000.0

  for r8 in [ hundredth, sqrt2, np.pi, million ]:

    n = 8
    d, p = r8_to_digits_decimal ( r8, n )

    print ( '' )
    print ( '  ', n, ' decimal digits of ', r8 )
    for i in range ( 0, n ):
      print ( i, d[i], 10**(p-i) )

  return

def r8_to_digits_sexagesimal ( r, n ):

#*****************************************************************************80
#
## r8_to_digits_sexagesimal() computes N sexagesimal digits representing R.
#
#  Discussion:
#
#    R = d(0) * 60^p + d(1) * 60^(p-1) + ... + d(n-1) * 60^(p+1-n) approximately
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    21 June 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real R, the number to be represented.
#
#    integer N: the number of digits to compute.
#
#  Output:
#
#    integer D(N): the sexagesimal digits.
#
#    integer P: the power of 60 multiplying the first digit.
#
  import numpy as np

  d = np.zeros ( n, dtype = np.int )
  p = 0

  if ( r == 0.0 ):
    return d, p
#
#  Try to deal with negative input.
#
  if ( r < 0.0 ):
    sgn = -1
    r = - r
  else:
    sgn = +1
#
#  Determine P such that 60^P <= R < 60^(P+1)
#
  rcopy = r
  base = 60.0
  while ( base <= rcopy ):
    rcopy = rcopy / base
    p = p + 1

  while ( rcopy < 1.0 ):
    rcopy = rcopy * base
    p = p - 1
#
#  Now read off N digits.
#
  rcopy = r
  div = base ** p
  for i in range ( 0, n ):
    d[i] = int ( np.floor ( rcopy / div ) )
    rcopy = ( rcopy % div )
    div = div / base
#
#  Multiply by sgn.
#
  d = sgn * d

  return d, p

def r8_to_digits_sexagesimal_test ( ):

#*****************************************************************************80
#
## r8_to_digits_sexagesimal_test() tests r8_to_digits_sexagesimal().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    21 June 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'r8_to_digits_sexagesimal_test():' )
  print ( '  r8_to_digits_sexagesimal() converts a real number' )
  print ( '  to its representation in base 60;' )

  hundredth = 1.0 / 100.0
  sqrt2 = np.sqrt ( 2.0 )
  million = 1000000.0

  for r8 in [ hundredth, sqrt2, np.pi, million ]:

    n = 8
    d, p = r8_to_digits_sexagesimal ( r8, n )

    print ( '' )
    print ( '  ', n, ' sexagesimal digits of ', r8 )
    for i in range ( 0, n ):
      print ( i, d[i], 60**(p-i) )

  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

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

  return

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