\ arithmetic on big signed-magnitude numbers
\ requires dmath.fth

\ based on code by leonard francis zettel, jr.
\ https://www.taygeta.com/fsl/library/big.fth
\ released to the forth scientific library.

\ this is a forth implementation of the "classical algorithms".
\ see knuth, donald the art of computer programming vol 2 p 250.

\ the internal representation of the big numbers is "little-endian
\ signed magnitude".  the cell at the address addr contains n, the
\ size of the number in digits.  n is positive for positive numbers,
\ negative for negative.  each succeeding cell contains a digit of
\ the number in base 2**(cell size-1), least significant digit first.

base @ decimal \ housekeeping

\ helper words
: leave unloop exit ;

32767 constant biggest            \ largest representable signed number
biggest 1+ constant bigbase       \ big number base as an unsigned number
biggest s>d 1 m+ 2constant bigbd  \ big number base as a double number

: >body ( xt -- a ) ;  \ Do nothing?
: cell- ( a1 -- a2 ) 1 cells - ;
: digit ( c n1 -- n2 f ) \ returns valid
    over [char] 0 <
    if   2drop false   \ characters below the zero character can't be digits
    else over [char] : <
         if   drop [char] 0 - true
         else over [char] a <
              if   2drop false
              else swap [char] 7 -        \ convert to numeric value
                   dup rot <
                   if   true              \ valid digit
                   else drop false
                   then
              then
         then
    then ;

: carry ( digit -- carry digit)  
  \ check for a carry, remove it,  it under the result.
  biggest over u< if bigbase - else 0 then swap ;

: d>carry ( low high -- carry digit) \ convert a double number to a
  \ low-order digit and a carry.
  bigbase um/mod swap ;

: overflow? ( borrow uj -- uj new_borrow) \ if uj is negative (indicating
  \ a result out of range on the previous subtraction), bring it in
  \ range and increment the borrow that will be necessary on the next
  \ digit.
  dup 0< if  bigbase +  1  else 0  then rot + ;

: big_digit_pointer ( -- ) ( n -- a )  \ create a word <name>.
  \ when <name> is executed return the address of the nth cell after
  \ the address in <name>'s data field.
  create 1 cells allot           \ create the word & allot the data space
  does> @                        \ put the address in <name>'s data field
                                 \ on the stack
        swap cells + ;           \ increment the address by index cells


: to_pointer ( <name> -- a1 ) \ compiling: a1 is <name>'s data field.
             ( a2 --) \ execution: a2 is placed in <name>'s
                         \ data field.
   ' >body postpone literal  postpone ! ; immediate
.s cr

\ TODO figure out what the heck is happening here and fix it 
\ I think the issue is the word represntation is different here
big_digit_pointer clippee
.s
\ broken 
: clip ( addr --)      \ remove leading zeroes from the number at addr
  to_pointer clippee
  0                    \ default - no non-zero digits
  1  0 clippee @ abs   \ loop from present number of digits to one.
  do       i clippee @               \ next big digit
           0<> if drop i leave then  \ index of first non-zero digit
                                     \ on stack
  -1 +loop
  ?dup if   0 clippee @              \ original sign & size
            0< if negate then        \ minus sign on new size
       else 1                        \ number is exactly zero, keep one
                                     \ of the zeros, show plus number
       then
  0 clippee ! ;                 \ store new size
\ something added to stack

: big_digit ( addr n1 -- n2k) \ return digit n1 of the big number at addr.
  \ if n1 is greater than the number of digits, return a leading zero.
  over @ abs         \ number of digits
  over <
  if   2drop 0       \ return leading zero
  else cells + @     \ return digit
  then ;

: bignegate ( addr --) \ change the sign of the big number at addr
  dup @ dup                          \ number of digits
  1 = if   over cell+ @              \ check for zero
           if   negate swap !        \ non-zero, negate
           else 2drop                \ zero, do nothing
           then
      else negate swap !
      then ;

: bigabs ( addr --) \ give the big number at addr its absolute value,
  dup @ abs swap ! ;

: big>here ( addr --)    \ "big to here" append the big number at addr
                         \ to the end of data space.
  here                   \ address to move to
  over @ abs 1+ cells    \ number of address units in the number
  dup allot              \ allot space for the number
  cmove ;


: adjust_sign ( a1 a2 a3 -- a3)  \ adjust the sign of the big
  \ number at a3 according to the rules for forming the algebraic
  \ product from the operands at a1 and a2
  rot @  rot @  xor 0< if dup bignegate then ;

\ move the number at a1 to a2 and free any data space beyond it.
: reposition ( a1 a2 -- )
  swap 2dup @                    \ (a2 a1 a2 size)
  abs 1+ cells                   \ (a2 a1 a2 bytes)
  dup >r  cmove                   \ (a2) (bytes)
  r> + here - allot ;

big_digit_pointer |big|1  big_digit_pointer |big|2

: |big|= ( a1 a2 -- flag) \ true if the big number at a1 has the
  \ same absolute value as the big number at a2.  false otherwise
  over @ abs over @ abs =     \  are the numbers the same size?
  if   to_pointer |big|1
       to_pointer |big|2
       true                   \ default initial flag.
       1  0 |big|1 @  abs
       do i |big|1 @
          i |big|2 @  <>
          if drop false leave then
          -1
       +loop
  else 2drop false
  then ;

: |big|< ( a1 a2 -- flag) \ true if the absolute value of the
  \ big number at a1 is less than the absolute value of the big number
  \ at a2.  false otherwise.

  to_pointer |big|2
  to_pointer |big|1
  0 |big|1 @ abs
  0 |big|2 @ abs
  2dup <
  if   2drop true
  else = false                   \ default flag if equal, result if <>.
       swap
       if   1 0 |big|1 @ abs     \ from the high order digit to the first
            do                   \ digit
               i |big|1 @
               i |big|2 @
               2dup
               <> if < nip leave  then
               2drop
            -1 +loop
       then
  then ;


: big0= ( addr -- flag) \ return true if the big number at addr is zero.
  dup @ 1 = if cell+ @ 0= else drop false then ;

: big0<> ( addr -- flag) \ return true if the big number at addr is not zero.
  big0= 0= ;

: big0< ( addr -- flag) @ 0< ;

: big< ( a1 a2 -- flag) \ true if the operand at a1 is less than
                            \ the operand at a2.  false otherwise.
  over @ over @ <           \ look at operand sign & number of digits
  if   2drop true
  else                      \ ( a1 a2)
       over @ over @ >
       if   2drop false
       else                 \ to get here the operands must be the
                            \ same sign & be of equal length
            dup @ 0<
            if swap then    \ if the numbers are negative, the one with
                            \ the larger absolute value is the lesser.
            dup @ abs >r    \ park number of digits
            r@ cells
            dup rot +       \ high order cell of operand 2.
            rot rot +       \ high order cell of operand 1.
            swap
            false           \ dummy initial flag
            r> 0
            do                      \ ( a1 a2 flag)
               drop                 \ flag from previous cycle
               over @ over @        \ ( a1 a2 digit1 digit2)
               < dup if leave then
               rot cell- rot cell-
               rot
            loop
            nip nip
       then
 then
 ;
: big= ( a1 a2 -- flag) \ true if the big number at a1 has the
  \ same absolute value as the big number at a2.  false otherwise
  over @  over @  =     \  are the numbers the same size?
  if   to_pointer |big|1
       to_pointer |big|2
       true                   \ default initial flag.
       1  0 |big|1 @  abs
       do i |big|1 @
          i |big|2 @  <>
          if drop false leave then
          -1
       +loop
  else 2drop false
  then ;


\ words doing mixed single-precision and big number arithmetic

big_digit_pointer big_addend

: big+s ( addr n --) \ add n to the number at addr.  n must be non-negative
  \ the number at addr must be non-negative and end at here.
  swap to_pointer big_addend  \ ( n)
  0 big_addend @  abs 1+      \ loop limit
  1                           \ loop start
  do                          \ ( n)
       i big_addend @  +      \ ( ui+n)
       carry
       i big_addend !         \ store new ui
       dup 0= if leave then   \ no carry, we are done
  loop
                              \ carry in high-order digit?
  if
      1 ,                     \ append carry to the number
      1  0 big_addend +!      \ increment number size
  then ;

big_digit_pointer big_multiplicand
: big*s ( addr n -- )         \ multiply the number at addr by n.
                              \ n must be positive
                              \ the number at addr must end at "here"
  swap
  to_pointer big_multiplicand
  0                           \ ( n carry)
  0 big_multiplicand @ abs 1+
  1
  do                          \ ( n carry)
       over
       i big_multiplicand @
       m*                     \ ( carry n low[ui*n] high[ui*n])
       rot  m+                \ ( n low[ui*n+carry] high[ui*n+carry])
       d>carry                \ ( n carry ui*n)
       i big_multiplicand !   \ store digit i back in u ( n carry)
  loop
  nip
  ?dup if  0 big_multiplicand  @  \ ( carry n)
            dup 0< if   1-
                   else 1+
                   then           \ ( carry n)
            0 big_multiplicand  ! ,
       then ;


big_digit_pointer big_dividend

: big/mods ( addr n1 -- n2) \ "big slash-mod s".  divide the big number at
  \ addr by n1, leaving the quotient at addr.  n2 is the remainder.
  swap
  to_pointer big_dividend
  0                         \ ( divisor remainder)
  1  0 big_dividend @ abs
  do                        \ ( divisor remainder)
      bigbase um*           \ ( divisor lowr highr)
      i big_dividend @      \ ( divisor divisor lowr highr uj)
      m+
      2 pick
      um/mod                \ ( divisor r wj )
      i big_dividend !
      -1
  +loop
  nip  0 big_dividend clip ;


\ words for going from characters to big numbers

: >big_number ( a1 a2 u1 -- a1 a3 u2) \ "to big number"
  \ extend the big number at a1 by the number represented by the
  \ string of u1 characters at a2.  a3 is the address of the first
  \ unconverted character and u2 is the number of unconverted characters
  2dup + >r                    \ address just beyond end of string on
                               \ return stack
  0 do                         \ ( a1 a2)
        2dup c@                \ ( a1 a2 a1 char)
        base @ digit           \ ( a1 a2 a1 n flag)
        if   over base @ big*s \ ( a1 a2 a1 n)
             big+s             \ ( a1 a2)
        else                   \ ( a1 a2 a1 char)
             drop leave \ ( a1 a2)
        then
        1+
    loop
    r> over - ;


: make_big_number ( a1 u -- a2)  \ convert the u characters at a1
                                       \ to a big number at a2
  \ if the first character is "-" (ascii 45) the result will be negative.
  \ embedded commas are ignored 
  \ (usa representation convention for large numbers)
  \ conversion stops at the first non-convertible character.
             
  over c@                          \ get the first character
  [char] - =                       \ is it a minus sign?
  dup >r
  if   swap 1+ swap 1- then     \ adjust to next character
                                   \ ( a1 u)
  here 1 , 0 ,                     \ create big number = 0
                                   \ ( a1 u a2)
  rot rot
  begin                            \ ( a2 a1 u)
        >big_number                \ ( a2 a1 u)
        over c@ [char] , =         \ ( a2 a1 u flag)
        over and
  while
        swap 1+ swap 1-
  repeat
  2drop
  r> if dup bignegate then ;


\ words for big number output
\ the words <big# through big#s and big. are adapted from
\ descriptions of their pictured numeric output string counterparts
\ in "all about forth" 2nd ed by glen haydon.  used with permission.

create big_string 256 allot

variable bighld

: <big# ( --) \ "less big number sign"
              \ initialize the big number pictured numeric output area
  big_string 256 + bighld ! ;  \ haydon p 67

: bighold ( c -- ) \ add c to the beginning of the big pictured numeric
                   \ output string
  -1 bighld +!  bighld @ c! ; \ haydon p 170.

: #big> ( a1 -- a2 +n) \ "number sign big less".  end big number
                             \ pictured output conversion
  drop  bighld @             \ start of string
  big_string 256 +     \ one past end of string
  over - 1 / ;         \ length of string

: big# ( addr -- addr)  \ "big number sign"
  \ generate the next ascii character from the big number at addr.
  \ afterward the big number at addr will hold the quotient obtained
  \ by dividing its previous value by the value in base.
  \ this result can then be used for further processing.
  \ haydon p 18

  dup
  base @ big/mods       \ next digit
  9 over <              \ is it bigger than a decimal digit?
  if 7 + then           \ add seven to its character representation,
                        \ thus skipping the ascii codes between 9 and a.
  48 +                  \ convert from number to ascii character code.
  bighold ;             \ add the character to the front of the output
                        \ string

: big#s ( addr -- addr) \ "big number sign s"  convert all digits of the
                        \ big number at addr to big numeric output, leaving
                        \ zero at addr
  begin big# dup @ 1 =      \ down to length 1
        over cell+ @ 0=     \ remaining cell is zero
        and
  until ;               \ haydon p 21.

: bigsign ( n --) \ put a minus sign in the big pictured numeric
                   \ character output string if n is negative
  0< if 45 bighold then ; \ haydon p 222.

: bigstring ( a1 sign --  ) \ display the big number at a1 with
  \ the sign of the number in sign.

  <big# big#s swap bigsign #big> type ;

\ words doing arithmetic on two big numbers

big_digit_pointer long_addend  big_digit_pointer short_addend

: sum ( a1 a2 - a3) \ a3 has the result of adding the absolute
  \ value of the big number at a1 to the absolute value of the big
  \ number at a2.

  over @ abs  over @ abs <     \ compare the size of the addends
  if swap then
  to_pointer short_addend
  to_pointer long_addend

  here                         \ address of result
  0 ,                          \ dummy placeholder for the count of the
                               \ result
  0                            \ initialize carry

  0 short_addend @ abs 1+      \ for each digit in the short addend
  1                            \ starting at the first
  do
       i short_addend @  +     \ add digit to carry
       i long_addend  @  +     \ add digit to previous sum
       carry ,                 \ new carry, append digit to result
  loop

  0 long_addend @ abs          \ number of digits in long operand
  1+                           \ jog to make do end on last digit
  0 short_addend @ abs         \ number of digits in short operand
  1+                           \ jog to start do on first digit
                               \ not yet used
  ?do  i long_addend @  +      \ append any remaining digits to the
       carry ,                 \ result, rippling the carry as
  loop                         \ necessary

  0 long_addend @ abs          \ result size so far
  swap
  if 1 , 1+ then               \ if final carry, append to result,
                               \ bump size
  over ! ;                     \ store result size.

big_digit_pointer minuend  big_digit_pointer subtrahend

: difference ( a1 a2 -- a3 )
  \ a3 is the address of the difference of the absolute values of
  \ the big number at a1 and the big number at a2.

  here >r                       \ park address of result
  2dup |big|=
  if   2drop  1 , 0 ,           \ equal absolute values, result is zero
  else
       2dup |big|<
       if swap then
       to_pointer subtrahend
       to_pointer minuend
       0 minuend @ abs  ,       \ count of the result
       0                        \ initialize borrow

       0 minuend @ abs 1+       \ for each minuend digit
       1                        \ starting with the first
       do                       \ ( borrow)
          0                     \ next borrow
          i minuend @           \ get the ith minuend digit
          rot -                 \ subtract previous borrow
          overflow?
          swap                  \ ( borrow result)
          0 subtrahend
          i big_digit  -        \ subtract the ith subtrahend digit
          overflow?             \ ( result borrow)
          swap ,                \ append result
       loop

       drop                 \ get rid of final borrow (it will be zero)
       r@ clip              \ remove leading zeroes
  then
  r> ;                      \ address of result on stack