bigsqrt - Python vs Pike

[Skaperen]

this is a commandline program i wrote many years ago in both Python (my first Python program) and Pike. both are given below.  the Pike code is substantially faster.  this program does its work using big ints, so that ends up being what is compared.  give it one argument, the number to find the square root of.  if no arguments are given, it will compute the golden ratio.  it just keeps outputting digits to get more and more precise (never ending).

My code

#!/usr/bin/python
# -*- coding: ISO-8859-1 -*-
#-----------------------------------------------------------------------------
# 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
#
# 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()

#!/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;
}

[wavic]

What if you import gmpy2 instead of math module?

[nilamo]

Once you hit line 57, your program is done. The python version never computes anything, regardless of arguments, as you unconditionally return out of the main function.

[micseydel]

Once that return issue is resolved, I'd be curious about actual benchmarks, including using PyPy instead of CPython (which I imagine you're using).

[Skaperen]

line 57 should be indented.  it is in my copy of the code which works for me. how can i insert plain code, now?

Output:
lt1/forums /home/forums 3> head -n57 bigsqrt.py|tail -n3
    if len(num_parts) != 2:
        stderr.write( "Number needs 0 or 1 '.'\n" )
        return 1
lt1/forums /home/forums 4> 

when i c&p the code to this post editing window (from a terminal window in a different desktop), the "return" line was unindent by 4 spaces so that it was just as far as the "if".  i typed in 4 spaces to fix it (looks right in preview now).  maybe that is why the program code was inserted wrong.

---

the new code display refuses to scroll when doing a highlight select ... we need a "select all" or better yet, a "download" link or button.  the old way to display code did have a "select all".  a "download" would be better, IMHO.

---

i put the 2 bigsqrt sources on the website as links to avoid the indentation issues:

http://stratusrelay.com/phil/bigsqrt.py (4f5b52aa9d8043a2d3d9d090f04dec7e)

http://stratusrelay.com/phil/bigsqrt.pike (3928561f6465d036852701a7a2421b21)

[snippsat]

i put the 2 bigsqrt sources on the website as links to avoid the indentation issues:

It looks more like a problem in the way you store code on web.
I save bigsqrt.py and then try to open it PyCharm,Atom,NotePad++ and indentation is wrong in all.
The same happen with copy-paste from website.

[nilamo]

Since we're still talking about line 57, it looks like the "if" clause line is indented with spaces, while the body of that same if block is indented with tabs. Mixing the two causes... strange things to happen :p

[Skaperen]

Since we're still talking about line 57, it looks like the "if" clause line is indented with spaces, while the body of that same if block is indented with tabs.  Mixing the two causes... strange things to happen :p

yes, that could be causing it ... anyone got a good untabify script?

[nilamo]

Most text editors can convert back and forth.  Otherwise, it should be as simple as:

with open('python-file.py', 'r') as source:
    with open('python-file-fixed.py', 'w') as out:
        for line in source:
            if line.startswith('\t'):
                line = ' '*4 + line[1:]
            print(line, end='', file=out)

[Skaperen]

emacs is the only editor i know (besides tools like sed) and it has don untabify wrong before so i want to avoid that.

i tried that code:

Output:
lt1/forums /home/forums 6> py3 untab.py bigsqrt.py fixed.py
Traceback (most recent call last):
  File "untab.py", line 13, in <module>
    main()
  File "untab.py", line 7, in main
    for line in source:
  File "/usr/lib/python3.5/encodings/ascii.py", line 26, in decode
    return codecs.ascii_decode(input, self.errors)[0]
UnicodeDecodeError: 'ascii' codec can't decode byte 0xa9 in position 138: ordinal not in range(128)
lt1/forums /home/forums 7> cat untab.py
def main():
    from sys import argv
    fn = argv[1]
    tf = argv[2]
    with open(fn, 'r') as source:
        with open(tf, 'w') as out:
            for line in source:
                if line.startswith('\t'):
                    line = ' '*4 + line[1:]
                print(line, end='', file=out)
    return

main()

google tells me about the GNU expand command which i found is already in Ubuntu 16.04.1 LTS :)  it seems to work so i will attach the fixed file.  the attach system here does not like .pike so an untabify of that is not attached.