#!/usr/bin/python
# -*- coding: utf-8 -*-
#-----------------------------------------------------------------------------
# Copyright © 2009, by Phil D. Howard - all other rights reserved
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
#-----------------------------------------------------------------------------
# program       sqrt_big
#
# purpose       Calculate lots of digits for the square root of a given number
#               using a scaled integer adaptation of Newton's method.
#
# usage         sqrt_big [number [digits]]
#
# note          The number defaults to zero.  Since the square root of zero
#               is not defined, this program calculates the Golden Ratio as
#               a special case.
#
# note          The number of digits defaults to zero.  If it is zero, this
#               program keeps running until something stops it, like user
#               intervention or out of memory.
#-----------------------------------------------------------------------------
import sys
import math

stderr = sys.stderr
stdout = sys.stdout

def main():
    write_limit = 10000

    # Get 0 or 1 or 2 arguments (defaults: 0 and 0)
    this_arg = '0'
    digits_max = 0
    if len(sys.argv) > 1:
        this_arg = sys.argv[1]
        if len(sys.argv) > 2:
            digits_max = int(sys.argv[2])

    num_parts = this_arg.split('.')
    if len(num_parts) == 1:
        num_parts = num_parts + ['0']
    if len(num_parts) != 2:
        stderr.write( "Number needs 0 or 1 '.'\n" )
        return 1
    if len(num_parts[0]) < 1:
        num_parts[0] = '0'
    if len(num_parts[1]) < 8:
        num_parts[1] = (num_parts[1] + '00000000')[0:8]
    num_whole = int( num_parts[0] )
    num_fract = int( num_parts[1] )
    num_scale = 10 ** len( num_parts[1] )

    if num_whole < 0 or num_fract < 0:
        stderr.write( "Number is negative\n" )
        return 1

    if num_whole == 0 and num_fract == 0:
        goldratio =         1
        num       = 125000000
        num_scale = 100000000
    else:
        goldratio =         0
        num       = num_whole * num_scale + num_fract

    # Use floating point sqrt() function to make an estimate.
    num_float = float(num)
    num_float /= float(num_scale)
    num_float = math.sqrt(num_float)
    num_float *= float(num_scale)
    root = int(num_float)
    root_scale = num_scale
    root_str = str(root)[0:1]

    # Do 12 cycles with short precision growth.
    for n in range( 0, 12 ):
        root = ( ( num * root_scale * root_scale + num_scale * root * root ) * root_scale ) / ( root + root )
        root /= root_scale
        root_scale *= num_scale
        root_str = str(root)

    digits = 1 + len(str(root)) - len(str(root_scale))
    if goldratio:
        # Fake the addition of 0.5 for the Golden Ratio.
        stdout.write( '1.6' )
        digits_done = 2
    elif digits > 0:
        # Output just the whole digits to get the fraction point in there.
        stdout.write( root_str[0:digits] )
        stdout.write( '.' )
        digits_done = digits
    elif digits <= 0:
        # Output enough 0 digits to align the fraction.
        stdout.write( '0.' )
        stdout.write( '0'*(-digits) )
        digits_done = 0

    # Do remaining cycles with long precision growth.
    for n in range( 0, 60 ):
        root = ( ( num * root_scale * root_scale + num_scale * root * root ) * root_scale ) / ( root + root )
        root_scale = root_scale * root_scale * num_scale
        root_str = str(root)
        digits = len(root_str) * 7 / 8;

        # Output the digits, with a limit per write, up to the maximum.
        while digits_done < digits:
            count = digits - digits_done
            if count > write_limit:
                count = write_limit
            if digits_max > 0:
                remain = digits_max - digits_done
                if count > remain:
                    count = remain
            stdout.write( root_str[ digits_done : digits_done + count ] )
            digits_done += count
            if digits_max > 0 and digits_done >= digits_max:
                break
        if digits_max > 0 and digits_done >= digits_max:
            break
    stdout.write( '\n' )
    return 0

main()