//-----------------------------------------------------------------------------
// Copyright © 2003 - 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.
//-----------------------------------------------------------------------------
// package	libh/string
// homepage	http://libh.slashusr.org/
//-----------------------------------------------------------------------------
// author	Philip Howard
// email	libh at ipal dot org
// homepage	http://phil.ipal.org/
//-----------------------------------------------------------------------------
// This file is best viewed using a fixed spaced font such as Courier
// and in a display at least 120 columns wide.
//-----------------------------------------------------------------------------

#include <ctype.h>
#include <errno.h>
#include <math.h>
#include <stdio.h>

#include "string_lib.h"

__PROTO_BEGIN__
//-----------------------------------------------------------------------------
// function	str_to_d
//
// purpose	Convert an external representation of a double number to
//		an internal representation.
//
//		This implementation focuses on portability across internal
//		representations and accuracy, rather than speed.  If speed
//		is desired, it should come from an implementation custom
//		made for the specific platform and float format, probably
//		in assembly language.
//
// arguments	1 (const char *) string to decode
//		2 (char * *) where to store end pointer
//
// returns	(double) 0.0 and errno set : error
//		(double) converted value
//
// issues
//
// If the base digits given, disregarding the fraction point, represents a
// a number which overflows (not just exceeds precision of) the internal
// format then an undetected conversion error is likely.  This will generally
// require quite a number of digits, although as few as 38 in minimal but
// strictly conforming implementations.
//
// Suffixes of 'f', 'l', 'F', or 'L' are permitted, although the number is
// not actually rounded to the apparently expressed type precision under
// the premise that an external representation does not imply specific C
// implementation limits, and that double is preferred by the caller
// of this function.  The effect of the suffix is to move the end pointer
// on to the next character.
//-----------------------------------------------------------------------------
double
str_to_d (
    const char *	arg_string
    ,
    char * *		arg_endptr
    )
__PROTO_END__
{
    static double		num_d		[16] = {
	0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0,
	8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0
    };

    static double		bits		[4] = {
	4194304.0, 8388608.0, 16777216.0, 33554432.0
    };

    const char *		str_ptr		;
    double			p2		;
    double			p5		;
    double			tmp		;
    double			num_hi		;
    double			num_lo		;
    int				accum		;
    int				base		;
    int				ch		;
    int				digit		;
    int				exp_ch		;
    int				exponent	;
    int				neg		;
    int				neg_exp		;
    int				point		;
    int				scaling		;


    //----------------------------
    // Make sure we have a string.
    //----------------------------
    if ( ! arg_string ) {
	errno = EINVAL;
	return 0.0;
    }
    str_ptr = arg_string;

    //----------------------------------
    // Skip over any leading whitespace.
    //----------------------------------
    while ( isspace( * str_ptr ) ) {
	++ str_ptr;
    }

    //-----------------------------------------------------------
    // Base value sign, if present, must be first non-whitespace.
    //-----------------------------------------------------------
    neg = 0;
    if ( * str_ptr == '-' ) {
	neg = 1;
	++ str_ptr;
    } else if ( * str_ptr == '+' ) {
	++ str_ptr;
    }

    //---------------------------------------
    // Skip over insignificant leading zeros.
    // They do not mean octal here.
    //---------------------------------------
    digit = 0;
    if ( * str_ptr == '0' ) {
	digit = 1;
	while ( * str_ptr == '0' ) {
	    ++ str_ptr;
	}
    }

#ifdef ALLOW_NAN
#ifdef NAN
    //-----------------
    // Check for "NaN".
    //-----------------
    if ( ( str_ptr[0] == 'N' || str_ptr[0] == 'n' ) &&
	 ( str_ptr[1] == 'a' || str_ptr[1] == 'A' ) &&
	 ( str_ptr[2] == 'N' || str_ptr[2] == 'n' ) ) {
	if ( arg_endptr ) {
	    * (const char * *) arg_endptr = str_ptr + 3;
	}
	return neg ? -(NAN) : (NAN);
    }
#endif
#endif

    //--------------------------------------
    // Validate first significant character.
    //--------------------------------------
    ch = * str_ptr;
    if ( '1' <= ch && ch <= '9' ) {
	base = 10;
    } else if ( ch == '.' ) {
	base = 10;
    } else if ( ch == 'x' || ch == 'X' ) {
	base = 16;
	ch = * ++ str_ptr;
    } else {
	if ( arg_endptr ) {
	    * (const char * *) arg_endptr = str_ptr;
	}
	if ( ! digit ) {
	    errno = EINVAL;
	}
	return 0.0;
    }

    //-------------------------------------------------------------------------
    // Now use different code for different bases.
    //-------------------------------------------------------------------------
    num_hi = 0.0;
    num_lo = 0.0;
    scaling = 0;
    point = '.';
    accum = 1;

    //-------------------------------------------------------------------------
    // Decimal.
    //-------------------------------------------------------------------------
    if ( base == 10 ) {

	for (;;) {
	    if ( ch == point ) {
		point = -1;
	    }
	    else if ( '0' <= ch && ch <= '9' ) {
		if ( accum ) {
		    tmp = num_hi;
		    num_hi *= 10.0;
		    if ( num_hi / 10.0 != tmp ) {
			num_hi = tmp;
			accum = 0;
		    } else {
			num_lo *= 10.0;
		    }
		}
		if ( accum ) {
		    num_lo += num_d[ ch - '0' ];

		    //-- Move some bits over to num_hi before they get too high in num_lo.
		    if ( num_lo >= bits[3] ) { num_lo -= bits[3]; num_hi += bits[3]; }
		    if ( num_lo >= bits[2] ) { num_lo -= bits[2]; num_hi += bits[2]; }
		    if ( num_lo >= bits[1] ) { num_lo -= bits[1]; num_hi += bits[1]; }
		    if ( num_lo >= bits[0] ) { num_lo -= bits[0]; num_hi += bits[0]; }

		    if ( point != '.' ) ++ scaling;
		} else {
		    if ( point == '.' ) -- scaling;
		}
	    }
	    else {
		break;
	    }
	    ch = * ++ str_ptr;
	}

	//-- Check for an exponent specification.
	exponent = 0;
	if ( ch == 'e' || ch == 'E' ) {
	    neg_exp = 0;
	    ch = * ++ str_ptr;

	    //-- Check for exponent sign.
	    if ( * str_ptr == '-' ) {
		neg_exp = 1;
		ch = * ++ str_ptr;
	    } else if ( * str_ptr == '+' ) {
		ch = * ++ str_ptr;
	    }

	    //-- Scan exponent digits in decimal.
	    for (;;) {
		if ( '0' <= ch && ch <= '9' ) {
		    exponent *= 10;
		    exponent += ch - '0';
		} else {
		    break;
		}
		ch = * ++ str_ptr;
	    }
	    if ( neg_exp ) exponent = - exponent;
	}

	//-- Compensate exponent for point scaling.
	exponent -= scaling;

	//-- Rescale converted base value for exponent.
	p2 = 1.0;
	p5 = 1.0;
	if ( exponent < 0 ) {
	    while ( exponent <= -32 ) { p2 *= 4294967296.0; p5 *= 23283064365386962890625.0; exponent += 32; }
	    if ( exponent <= -16 ) { p2 *= 65536.0; p5 *= 152587890625.0; exponent += 16; }
	    if ( exponent <= -8 ) { p2 *= 256.0; p5 *= 390625.0; exponent += 8; }
	    if ( exponent <= -4 ) { p2 *= 16.0; p5 *= 625.0; exponent += 4; }
	    if ( exponent <= -2 ) { p2 *= 4.0; p5 *= 25.0; exponent += 2; }
	    if ( exponent <= -1 ) { p2 *= 2.0; p5 *= 5.0; exponent += 1; }
	    num_hi += num_lo;
	    num_hi /= p5;
	    num_hi /= p2;
	} else {
	    while ( exponent >= 32 ) { p2 *= 4294967296.0; p5 *= 23283064365386962890625.0; exponent -= 32; }
	    if ( exponent >= 16 ) { p2 *= 65536.0; p5 *= 152587890625.0; exponent -= 16; }
	    if ( exponent >= 8 ) { p2 *= 256.0; p5 *= 390625.0; exponent -= 8; }
	    if ( exponent >= 4 ) { p2 *= 16.0; p5 *= 625.0; exponent -= 4; }
	    if ( exponent >= 2 ) { p2 *= 4.0; p5 *= 25.0; exponent -= 2; }
	    if ( exponent >= 1 ) { p2 *= 2.0; p5 *= 5.0; exponent -= 1; }
	    num_hi *= p5;
	    num_lo *= p5;
	    num_hi += num_lo;
	    num_hi *= p2;
	}
    }

    //-------------------------------------------------------------------------
    // Cetal (hexadecimal).
    //-------------------------------------------------------------------------
    else if ( base == 16 ) {

	//-- Outer loop scans base digits before the fraction point.
	for (;;) {
	    if ( ch == '.' ) {
		//-- Inner loop scans base digits after the fraction point.
		for (;;) {
		    ch = * ++ str_ptr;
		    if ( ( '0' <= ch && ch <= '9' && ( ( digit = ch - '0' ), 1 ) ) ||
			 ( 'A' <= ch && ch <= 'F' && ( ( digit = 10 + ch - 'A' ), 1 ) ) ||
			 ( 'a' <= ch && ch <= 'f' && ( ( digit = 10 + ch - 'a' ), 1 ) ) ) {
			num_hi *= 16.0;
			num_hi += num_d[ digit ];
			++ scaling;
		    } else {
			break;
		    }
		}
		break;
	    }
	    if ( ( '0' <= ch && ch <= '9' && ( ( digit = ch - '0' ), 1 ) ) ||
		 ( 'A' <= ch && ch <= 'F' && ( ( digit = 10 + ch - 'A' ), 1 ) ) ||
		 ( 'a' <= ch && ch <= 'f' && ( ( digit = 10 + ch - 'a' ), 1 ) ) ) {
		num_hi *= 16.0;
		num_hi += num_d[ digit ];
	    } else {
		break;
	    }
	    ch = * ++ str_ptr;
	}

	//-- Check for an exponent specification.
	//-- 'p' or 'P' : radix 2 (C99)
	//-- 'q' or 'Q' : radix 16 (non-standard)
	exponent = 0;
	if ( ch == 'p' || ch == 'P' || ch == 'q' || ch == 'Q' ) {
	    exp_ch = ch;
	    neg_exp = 0;
	    ch = * ++ str_ptr;

	    //-- Check for exponent sign.
	    if ( * str_ptr == '-' ) {
		neg_exp = 1;
		ch = * ++ str_ptr;
	    } else if ( * str_ptr == '+' ) {
		ch = * ++ str_ptr;
	    }

	    //-- Scan exponent digits in decimal.
	    for (;;) {
		if ( '0' <= ch && ch <= '9' ) {
		    exponent *= 10;
		    exponent += ch - '0';
		} else {
		    break;
		}
		ch = * ++ str_ptr;
	    }
	    if ( neg_exp ) exponent = - exponent;

	    //-- If exponent radix is 'q' or 'Q' the exponent is radix 16 instead of 2.
	    if ( exp_ch == 'q' || exp_ch == 'Q' ) {
		exponent *= 4;
	    }
	}

	//-- Compensate exponent for radix 16 point scaling.
	exponent -= scaling * 4;

	//-- Rescale converted base value for exponent.
	p2 = 1.0;
	if ( exponent < 0 ) {
	    while ( exponent <= -64 ) { p2 *= 18446744073709551616.0; exponent += 64; }
	    if ( exponent <= -32 ) { p2 *= 4294967296.0; exponent += 32; }
	    if ( exponent <= -16 ) { p2 *= 65536.0; exponent += 16; }
	    if ( exponent <= -8 ) { p2 *= 256.0; exponent += 8; }
	    if ( exponent <= -4 ) { p2 *= 16.0; exponent += 4; }
	    if ( exponent <= -2 ) { p2 *= 4.0; exponent += 2; }
	    if ( exponent <= -1 ) { p2 *= 2.0; exponent += 1; }
	    num_hi /= p2;
	} else {
	    while ( exponent >= 64 ) { p2 *= 18446744073709551616.0; exponent -= 64; }
	    if ( exponent >= 32 ) { p2 *= 4294967296.0; exponent -= 32; }
	    if ( exponent >= 16 ) { p2 *= 65536.0; exponent -= 16; }
	    if ( exponent >= 8 ) { p2 *= 256.0; exponent -= 8; }
	    if ( exponent >= 4 ) { p2 *= 16.0; exponent -= 4; }
	    if ( exponent >= 2 ) { p2 *= 4.0; exponent -= 2; }
	    if ( exponent >= 1 ) { p2 *= 2.0; exponent -= 1; }
	    num_hi *= p2;
	}

    }

    //-------------------------------------------------------------------------
    // The base cannot be unknown unless there is some logic error.
    // EDOM is not a perfect choice here.  Any better ideas?  abort()?
    //-------------------------------------------------------------------------
    else {
	errno = EDOM;
	return 0.0;
    }

    //-----------------------------------------------------------
    // If the number is supposed to be negative, then make it so.
    //-----------------------------------------------------------
    if ( neg ) {
	num_hi = - num_hi;
    }

    //---------------------------------------------------------
    // Store the ending pointer and return the converted value.
    //---------------------------------------------------------
    if ( arg_endptr ) {
	* (const char * *) arg_endptr = str_ptr;
    }
    return num_hi;
}

