! math.h - Extended Mathematical routines.
!   Version 1.0 (19-Sep-2001)
!
! by Matt Albrecht - groboclown@users.sourceforge.net
!
! (If you need to edit this file, note that indentations are 4 spaces
! and tabs are not used.)
!
! This has been donated to the Public Domain.  Use, abuse, and don't blame me.
!
! To prevent dirtying the global namespace, all members of this file begin with
! "math_", while members private to this file begin with "math__".

! Due to unsigned stuff, requires Z-Machine version 5 or older.

! Note that there is a bug in Nitfol v0.5 - if you unsign divide or mod $8000
! by anything less than $8000, it returns the wrong answer (the negative
! of the correct answer).


System_file;



! C-like header!
Ifndef MATH__INCLUDED;
Constant MATH__INCLUDED;

Message "Adding Extended Math library";

! Depends upon the longint.h library
Ifndef LONGINT__INCLUDED;
Constant LONGINT__INCLUDED;
Include "longint";
Endif;



! Used for interior operations
Array math__long1 -> 4;
Array math__long2 -> 4;
Array math__long3 -> 4;


!-------------------------------------------------------------------------------
! Print routines
!-------------------------------------------------------------------------------

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Print the 16-bit word as unsigned.
[ math_unsigned
    x;  ! parameter
    
    if (x < 0)
    {
        ! use long longint library to print the sucker
        math__set_unsigned_word_to_long( x, math__long3 );
        print (longunsign)math__long3;
    }
    else
    {
        print x;
    }
];


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Displays a word as a hex string.
! (Stolen from N.G.'s Inform help pages)
!
! Argument 1: the word to convert to hex.
[ math_hex
    x   ! parameter
    y;  ! local
    
	y = (x & $ff00) / $100;
	x = x & $ff;
	print (math_hdigit) y/$10, (math_hdigit) y,
          (math_hdigit) x/$10, (math_hdigit) x;
];

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Displays 4 bits as a hex digit.
! (Stolen from N.G.'s Inform help pages)
!
! Argument 1: the 4 bits to convert to hex.
[ math_hdigit
    x; ! parameter
    
	x = x & $000f;
	if (x < 10) print x; else print (char)'a'+x-10;
];


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Displays a word as a binary string.
! (Stolen from N.G.'s Inform help pages)
!
! Argument 1: the word to convert to hex.
[ math_binary
    x   ! parameter
    y z;    ! locals
! $8000 is considered negative, so we can't use y > 0
!print "[binary of ",x,":";
    if (x < 0)
    {
!print "... first is negative: ";
        print "1";
    }
    else
    {
        print "0";
    }
    for (y = $4000 : y > 0 : y = y/2)
    {
        ! perform a binary AND
        @and x y -> z;
!print "... ",(hex)x," & ",(hex)y," = ",(hex)z,": ";
        print (z/y);
    }
!print "]^";
];




!-------------------------------------------------------------------------------
! Unsigned Math routines
!-------------------------------------------------------------------------------

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Shift the first parameter by the given number of bits in the second parameter,
! shifting left if it is positive, or right if it is negative.  Note that
! since this is unsigned, shifting right will not do sign extension.
[ math_unsigned_shift
    x count ! parameters
    r;
    
    @log_shift x count -> r;
    return r;
];


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Unsigned division of 2 16-bit words (returns x / y).
!
[ math_unsigned_div
    x y     ! parameters
    ret;    ! locals
    
    math__unsigned_div_mod( x, y );

    ret = math__get_unsigned_long_to_word( math__long1 );
!print (math_unsigned)x," / ",(math_unsigned)y," = ",(math_unsigned)ret,"^";
    return ret;
];



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Unsigned modulo of 2 16-bit words (returns x % y).
!
[ math_unsigned_mod
    x y     ! parameters
    ret;    ! locals
    
    math__unsigned_div_mod( x, y );

    ret = math__get_unsigned_long_to_word( math__long2 );
!print (math_unsigned)x," % ",(math_unsigned)y," = ",(math_unsigned)ret,"^";
    return ret;
];


!-------------------------------------------------------------------------------
! Signed Math routines
!-------------------------------------------------------------------------------

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Shift the first parameter by the given number of bits in the second parameter,
! shifting left if it is positive, or right if it is negative.  Note that
! since this is signed, shifting right will perform a sign extension.
[ math_signed_shift
    x count ! parameters
    r;
    
    @art_shift x count -> r;
    return r;
];


!---------------------------------------------------
[ math_square
        x;
    return x * x;
];


!---------------------------------------------------
[ math_squareRoot
        c       ! parameter
        prev diff last orig; ! locals
     ! Based on Newton-Raphson method: f(x) = x*x - c, f'(x) = 2*x.
     ! So, a[n+1] = a[n] - (f(a[n])/f'(a[n]))
     ! Therefore,
     !     a[n+1] = a[n] - (a[n]*a[n] - c)/(2*a[n])
     ! In order to avoid overflows of precision:
     !     a[n+1] = a[n] - (a[n]*a[n]/(2*a[n])) - (c/(2*a[n]))
     !            = a[n] - (a[n]/2) - (c/2/a[n])


     ! the square root of c is <= c/2 by definition of primes.
     orig = c;
     c = c / 2;
     prev = c;
     if (prev < 0) return -1; ! bad input value

     if (prev == 0) return 0;
     for (::) {
         last = prev;
         diff = prev/2 - c/prev;
         if (diff == 0) break;
         
         prev = prev - diff;
         
         ! There are situations where we can toggle between two values,
         ! never reaching a good result.
         if (last == prev+1) break;
     };
     return prev;
];


!-------------------------------------------------------------------------------
! Private Math routines
!-------------------------------------------------------------------------------

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
[ math__set_unsigned_word_to_long
    x longArray; ! parameters
    
    LongSet( longArray,
        0, 0,
        math_unsigned_shift( x, -8) & $ff, x & $ff );
!print "[",x," to unsigned long ",(longunsign)longArray,"]^";
];


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Puts a word value into the given longint array
[ math__set_signed_word_to_long
    x longArray    ! parameters
    top;           ! locals
    
    top = 0;
    if (x < 0)
    {
        top = $ff;
    }
    LongSet( longArray,
        top, top,
        math_unsigned_shift( x, -8) & $ff, x & $ff );
];


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Puts a word value into the given longint array
[ math__get_unsigned_long_to_word
    longArray      ! parameters
    x y ret;       ! locals
    
    x = longArray->2;
    y = longArray->3;
    ret = math_unsigned_shift( x, 8 ) + y;
!print "[converted unsigned long ",(longunsign)longArray," to ",
!(math_unsigned)ret,". ([0] = ",longArray->0,", [1] = ",longArray->1,", [2] = ",
!x,", [3] = ",y,")]^";
    return ret;
];




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Sets the result of unsigned (x / y) in math__long1
! and (x % y) in math__long2
!
[ math__unsigned_div_mod
    x y     ! parameters
    t;
    
    ! special cases for optimization
        ! include ~= $8000 for covering Nitfol 0.5 bug
    if (x > 0 && y > 0 && x ~= $8000 && y ~= $8000)
    {
!        print "Unsigned division!^";
        math__set_unsigned_word_to_long( x / y, math__long1 );
        math__set_unsigned_word_to_long( x % y, math__long2 );
    }
! Uncomment these if you feel the performance drop from the if-statements
! is worth the reduced computation for these few special cases.
!    else if (x == y)
!    {
!        math__set_unsigned_word_to_long( 1, math__long1 );
!        math__set_unsigned_word_to_long( 0, math__long2 );
!    }
!    else if (x == 0)
!    {
!        math__set_unsigned_word_to_long( 0, math__long1 );
!        math__set_unsigned_word_to_long( 0, math__long2 );
!    }
!    else if (y == 1)
!    {
!        math__set_unsigned_word_to_long( x, math__long1 );
!        math__set_unsigned_word_to_long( 0, math__long2 );
!    }
    else
    {
        ! long1 = x
        math__set_unsigned_word_to_long( x, math__long1 );
        ! long2 = y
        math__set_unsigned_word_to_long( y, math__long2 );
!print " - ",(math_unsigned)x," (",(longunsign)math__long1,") / ",
!(math_unsigned)y," (",(longunsign)math__long2,"): ";
    
        ! long1 = long1 / long2, long2 = long1 % long2
        LongUnsignDivMod( math__long1, math__long2,
            math__long1, math__long2 );
    }
    
    ! There is a bug in Nitfol v0.5 - see header for details
!    if (x == $8000 && y ~= $8000 && y > 0)
!    {
!        ! need to turn the quotient and mod positive
!        math__set_unsigned_word_to_long(
!            -math__get_unsigned_long_to_word( math__long1 ), math__long1 );
!        math__set_unsigned_word_to_long(
!            -math__get_unsigned_long_to_word( math__long2 ), math__long2 );
!    }
    
!print "quotient = ",(longunsign)math__long1,", remainder = ",
!    (longunsign)math__long2,".^";
];

Endif;