#!/usr/bin/env pike
//-----------------------------------------------------------------------------
// Copyright © 2009 - Philip Howard - All 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.pike
//
// purpose	Calculate lots of digits for the square root of a given number
//		using a scaled integer adaptation of Newton's method.
//
// usage	[pike] sqrt_big.pike [number [digits]]
//
// note		The number defaults to zero.  If it is zero, 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 it is out of memory.
//-----------------------------------------------------------------------------
int main( int argc, array(string) arga )
{
    array	num_parts;
    string	this_arg,root_str;
    float	num_float;
    int		num,num_whole,num_fract,num_scale;
    int		root,root_scale;
    int		digits,digits_done,digits_max;
    int		goldratio;
    int		write_limit;


    write_limit = 10000;

    // Get 0 or 1 or 2 arguments (defaults: 0 and 0)
    this_arg = argc <= 1 ? "0" : arga[1];
    digits_max = argc <= 2 ? 0 : (int) arga[2];

    // Break apart the argument number.
    num_parts = this_arg / ".";
    if ( sizeof( num_parts ) == 1 ) num_parts = num_parts + ({ "0" });
    if ( sizeof( num_parts ) != 2 ) return 1;
    if ( strlen( num_parts[0] ) < 1 ) num_parts[0] = "0";
    if ( strlen( num_parts[1] ) < 8 ) num_parts[1] = ( num_parts[1] + "00000000" )[0..7];
    num_whole = (int) num_parts[0];
    num_fract = (int) num_parts[1];
    num_scale = 10->pow( strlen( num_parts[1] ) );

    // Do not calculate square root of negative number.
    if ( num_whole < 0 ) return 1;
    if ( num_fract < 0 ) return 1;

    // If zero is specified, calculate the golden ratio instead.
    if ( num_whole == 0 && num_fract == 0 ) {
	// To find golden ratio, find square root of 1.25 and add 0.5 to result.
	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 = sqrt( num_float );
    num_float *= (float) num_scale;
    root = (int) num_float;
    root_scale = num_scale;
    root_str = ( (string) root )[0..0];

    // Do 13 cycles with short precision growth.
    for ( int n = 0; n < 13; ++ n ) {
	root = ( ( num * root_scale * root_scale + num_scale * root * root ) * root_scale ) / ( root + root );
	root /= root_scale;
	root_scale *= num_scale;
	root_str = (string) root;
    }

    // Display the whole part first.
    digits = 1 + strlen( (string) root ) - strlen( (string) root_scale );
    if ( goldratio ) {
	// Fake the addition of 0.5 for the Golden Ratio.
	write( "1.6" );
	digits_done = 2;
    } else if ( digits > 0 ) {
	// Output just the whole digits to get the fraction point in there.
	write( "%s.", root_str[ 0 .. digits - 1 ] );
	digits_done = digits;
    } else if ( digits <= 0 ) {
	// Output enough 0 digits to align the fraction.
	write( "0.%s", "0" * (-digits) );
	digits_done = 0;
    }

    // Do remaining cycles with long precision growth.
    for ( int n = 0; n < 60; ++ n ) {
	string root_str;
	root = ( ( num * root_scale * root_scale + num_scale * root * root ) * root_scale ) / ( root + root );
	root_scale = root_scale * root_scale * num_scale;

	root_str = (string) root;
	digits = ( strlen( root_str ) * 14 ) / 15;

	// Output the digits, with a limit per write, up to the maximum.
	while ( digits_done < digits ) {
	    int count;

	    count = digits - digits_done;
	    if ( count > write_limit ) count = write_limit;
	    if ( digits_max > 0 ) {
		int remain;
		remain = digits_max - digits_done;
		if ( count > remain ) count = remain;
	    }
	    write( "%s", root_str[ digits_done .. digits_done + count - 1 ] );
	    digits_done += count;
	    if ( digits_max > 0 && digits_done >= digits_max ) break;
	}
	if ( digits_max > 0 && digits_done >= digits_max ) break;
    }

    write( "\n" );
    return 0;
}