nxu/fix16.tal

( fix16.tal                                               )
(                                                         )
( use a signed 16-bit short as a fixed point number.      )
(                                                         )
( numbers are interpreted as fractions with an implicit   )
( 256 denominator. the upper byte is signed and           )
( represents the "whole" part of the number, and the      )
( lower byte is unsigned and represents the               )
( "fractional" part of the number.                        )
(                                                         )
( 16-bit fixed point can represent fractional values      )
( in the range -128 <= x < 128. the smallest fraction it  )
( can represent is 1/256, which is about 0.0039.          )
(                                                         )
( SHORT     FRACTION     DECIMAL                          )
( #0000        0/256      0.0000                          )
( #0001        1/256      0.0039                          )
( #0002        2/256      0.0078                          )
( #0040       64/256      0.2500                          )
( #0080      128/256      0.5000                          )
( #0100      256/256      1.0000                          )
( #0700     1792/256      7.0000                          )
( #7f00    32512/256    127.0000                          )
( #7fff    32767/256    127.9961                          )
( #8000      invalid     invalid                          )
( #8001   -32767/256   -127.9961                          )
( #8100   -32767/256   -127.0000                          )
( #ff00     -256/256     -1.0000                          )
( #ffff       -1/256     -0.0039                          )
(                                                         )
( due to the limited range operations saturate at the     )
( terminal values (i.e. #7fff and #8001). #8000 should    )
( never be generated, and can be considered a Nan value   )
( although this library does not use it as such.          )
(                                                         )
( many 8.8 operations are equivalent to unsigned int16:   )
(   * addition                                            )
(   * subtraction                                         )
( or signed int16:                                        )
(   * comparisons/equality                                )
(                                                         )
( but due to 16-bit truncation multiplication differs...  )
(                                                         )
( x*y = x0*y0 + x0*y1/256 + x1*y0/256 + x1*y1/65536       )
(                                                         )
( since we only have 16-bits:                             )
(   1. we need to drop the 8 high bits from x0*y0         )
(   2. we need to drop the 8 low bits from x1*y1          )
(   3. we need to use all the bits from x0*y1 and x1*y0   )
(                                                         )
( that said, if either x or y is whole (i.e. ends in 00)  )
( then we can just shift that argument right by 8 and use )
( MUL2.                                                   )
(                                                         )
( similarly with division we have:                        )
(                                                         )
(     x = x'/256                                          )
(     y = y'/256                                          )
(   x/y = z = z'/256                                      )
(   (x'/256)/(y'/256) = (x'*256 / y)/256                  )
(    z' = (x' * 256 / y)/256                              )
(                                                         )
( trigonometry, etc:                                      )
(  - sin() and cos() are supported                        )
(  - tan() is mostly supported                            )
(  - tan(#0192) and tan(#04b6) throw an error             )
(  - a few tan() values are inaccurate due to range       )
(    + values "next to" pi/2 and 3pi/2 are affected       )
(    + tan(#0191) returns 127.996 not 227.785             )
(    + tan(#0193) returns -127.996 not -292.189           )
(    + etc.                                               )
(  - log() is supported                                   )
(  - log(0) throws an error                               )
(  - log values farther from zero may be a bit inaccurate )
(                                                         )
( rounding, truncation, conversion:                       )
(                                                         )
(  - x16-to-s16: converts the whole to signed 16-bit,     )
(                discarding the fractional part. the      )
(                magnitude of numbers is never increased. )
(                                                         )
(  - x16-to-s8:  same as x16-to-s16 but returns one byte. )
(                                                         )
(  - x16-ceil:   returns a x16 value with any fractional  )
(                part rounded up.                         )
(                                                         )
(  - x16-floor:  returns a x16 value with any fractional  )
(                part rounded down.                       )
(                                                         )
(  - x16-round:  returns a x16 value with any fractional  )
(                part rounded up or down depending on the )
(                fractional part. in the case of numbers  )
(                with a fractional part of 0.5, they will )
(                round towards the nearest even value.    )

( useful constants                       )
(                                        )
( to generate your own:                  )
(                                        )
(   1. take true value, e.g. 3.14159...  )
(   2. multiply by 256                   )
(   3. round to nearest whole number     )
(   4. emit hex output                   )
(                                        )
( in python: hex(round(x * 256))         )
%x16-zero      { #0000 } (    0.0        )
%x16-one       { #0100 } (    1.0        )
%x16-two       { #0200 } (    2.0        )
%x16-ten       { #0a00 } (   10.0        )
%x16-hundred   { #6400 } (  100.0        )
%x16-minus-one { #ff00 } (   -1.0        )
%x16-minus-two { #fe00 } (   -2.0        )
%x16-pi/2      { #0192 } (    1.57079... )
%x16-pi        { #0324 } (    3.14159... )
%x16-3pi/2     { #04b6 } (    4.71239... )
%x16-pi*2      { #0648 } (    6.28318... )
%x16-e         { #02b8 } (    2.71828... )
%x16-phi       { #019e } (    1.61803... )
%x16-sqrt-2    { #016a } (    1.41421... )
%x16-sqrt-3    { #01bb } (    1.73205... )
%x16-epsilon   { #0001 } (    0.00390... )
%x16-minimum   { #8000 } ( -128.0        )
%x16-maximum   { #7fff } (  127.99609... )
%x16-max-whole { #7f00 } (  127.0        )

( utils )
@x16-is-non-neg ( x* -> bool^ ) x16-minimum LTH2 JMP2r
@x16-is-neg ( x* -> bool^ ) x16-maximum GTH2 JMP2r
@x16-emit-dec-digit ( d^ -> ) #30 ADD #18 DEO JMP2r
@error [ #0000 DIV ]

@x16-emit ( x* -> )
    DUP2 #8000 EQU2 ?&is-min
    DUP2 #8000 GTH2 ?&is-neg
    SWP DUP #64 LTH ?&<100
        #64 DIVk DUP x16-emit-dec-digit MUL SUB !&>=10
    &is-min POP2
        [ LIT "- #18 DEO LIT "1 #18 DEO LIT "2 #18 DEO LIT "8 #18 DEO ]
        [ LIT ". #18 DEO LIT "0 #18 DEO LIT "0 #18 DEO LIT "0 #18 DEO ]
        JMP2r
    &is-neg
        [ LIT "- #18 DEO ] #ffff EOR2 INC2 !x16-emit
    &<100 DUP #0a LTH ?&<10
        &>=10 #0a DIVk DUP x16-emit-dec-digit MUL SUB
    &<10  x16-emit-dec-digit
        [ LIT ". #18 DEO ]
    ( emit fractional part )
    #00 SWP ( lo* )
    #000a MUL2 #0100 DIV2k DUP x16-emit-dec-digit MUL2 SUB2
    #000a MUL2 #0100 DIV2k DUP x16-emit-dec-digit MUL2 SUB2
    #000a MUL2 #0100 DIV2k DUP x16-emit-dec-digit MUL2 SUB2
    #000a MUL2 #0100 DIV2k STH2k MUL2 SUB2 #0080 LTH2 ?&no-round INC2r
    &no-round STH2r NIP !x16-emit-dec-digit

( comparison between x and y. )
(  - ff: x < y )
(  - 00: x = y )
(  - 01: x > y )
@x16-cmp ( x* y* -> c^ )
    STH2k x16-maximum GTH2 ?&yn      ( x* [y*] ; ? )
        DUP2 x16-minimum LTH2 ?&same ( x* [y*] ; y>=0 )
              POP2 POP2r #ff JMP2r   ( -1      ; x<0 y>=0 )
    &yn DUP2 x16-maximum GTH2 ?&same ( x* [y*] ; y<0 )
              POP2 POP2r #01 JMP2r   ( 1       ; x>=0  y<0 )
        &same STH2r ( fall-thru )    ( res     ; x<0  y<0 b )

( unsigned comparison between x and y. )
(  - ff: x < y )
(  - 00: x = y )
(  - 01: x > y )
@x16-ucmp ( x* y* -> c^ )
    LTH2k ?&lt GTH2 JMP2r
    &lt POP2 POP2 #ff JMP2r

@x16-eq   ( x* y* ->  x=y^ ) EQU2 JMP2r
@x16-ne   ( x* y* -> x!=0^ ) NEQ2 JMP2r
@x16-lt   ( x* y* ->  x<y^ ) x16-cmp #ff EQU JMP2r
@x16-lteq ( x* y* ->  x<y^ ) x16-cmp #01 NEQ JMP2r
@x16-gt   ( x* y* ->  x<y^ ) x16-cmp #01 EQU JMP2r
@x16-gteq ( x* y* ->  x<y^ ) x16-cmp [ #ff NEQ ] JMP2r

@unsigned-min ( x* y* -> min* )
    LTH2k JMP SWP2 POP2 JMP2r

@x16-min ( x* y* -> min[x, y]* )
    OVR2 OVR2 x16-lt JMP SWP2 POP2 JMP2r

@x16-max ( x* y* -> max[x, y]* )
    OVR2 OVR2 x16-gt JMP SWP2 POP2 JMP2r

@x16-is-whole ( x* -> bool^ )
    NIP #00 EQU JMP2r

@x16-add ( x* y* -> x+y* )
    ADD2 JMP2r

@x16-sub ( x* y* -> x-y* )
    SUB2 JMP2r

@x16-negate ( x* -> -x* )
    #0000 SWP2 SUB2 JMP2r

@x16-mul ( x* y* -- xy* )
    ;x16-mul-unsigned !x16-signed-op

@x16-mul8 ( x^ y^ -> xy* )
    #0000 SWP2 ROT SWP MUL2 JMP2r

@x16-mul-unsigned ( ab* cd* -> ac+ad+bc+bd* )
    OVR2 OVR2 STH2 STH2 ROT SWPr ROTr             ( a c d b [c b a d] )
    SWP2  x16-mul8 DUP2 #007f GTH2 ?&o1 SWP       ( d b acc* [c b a d] )
    SWP2  x16-mul8 #08 SFT2 ADD2                  ( acc* [c b a d] )
    STH2r x16-mul8 ADD2 OVR #7f GTH ?&o3          ( acc* [c d] )
    STH2r x16-mul8 ADD2 OVR #7f GTH ?&o4 JMP2r    ( acc* )
    &o1 POP2 POP2r &o3 POP2r &o4 POP2 #7fff JMP2r ( ; handle overflow )

@x16-div ( x* y* -- x/y* )
    ;x16-div-unsigned !x16-signed-op

( The idea here is a bit complicated. )
( )
( First we look at `x y DIV2`. If that is >255 then we are )
( going to overflow and we can return 7fff and be done. )
( )
( If not, we multiply that result by 256 and start looking )
( for smaller and smaller fractions of y to star chipping )
( away at the remainder, if any. )
( )
( However we want to avoid rounding errors caused by needlessly )
( truncating y. For that reason, if the remainder is less than )
( 0x8000 we will double it rather than diving y. Once it cannot )
( be safely multiplied we will start dividing y by 2. After 16 )
( rounds of multiplying the remainder or dividing the quotient )
( we get our final answer. )
@x16-div-unsigned ( x* y* -> x/y* )
    DIV2k DUP2 #007f GTH2 ?&o   ( x* y* x/y* )
    STH2k LITr 80 SFT2r         ( x* y* x/y* [div=(x/y)<<8*] )
    OVR2 STH2                   ( x* y* x/y* [div* y*] )
    MUL2 SUB2                   ( x%y* [div* y*] )
    STH2r LIT2r 0100            ( x%y* y* [div* 0100*] )

    ( We know rem < y, so start left-shifting rem. )
    &ploop                       ( rem* y* [div* s*] )
        STH2kr #0000 EQU2 ?&done ( rem* y* [div* s*] ; done when s=0 )
        OVR2 #7fff GTH2 ?&loop   ( rem* y* [div* s*] ; rem too big, start shifting y )
        SWP2 #10 SFT2 SWP2       ( rem<<1* y* [div* s*] )
        LITr 01 SFT2r            ( rem<<1* y* [div* s>>1*] )
        LTH2k ?&ploop            ( rem<<1* y* [div* s>>1*] ; rem too small )
        SWP2 OVR2 SUB2 SWP2      ( rem<<1-y* y* [div* s>>1*] )
        DUP2r ROT2r ADD2r SWP2r  ( rem<<1-y* y* [div+s>>1* s>>1*] )
        !&ploop                  ( rem<<1-y* y* [div+s>>1* s>>1*] )

    ( We know rem < y, so start right-shifting y. )
    &loop                        ( rem* y* [div* s*] )
        STH2kr #0000 EQU2 ?&done ( rem* y* [div* s*] ; done when s=0 )
        #01 SFT2 LITr 01 SFT2r   ( rem* y>>1* [div* s>>1*] )
        LTH2k ?&loop             ( rem* y>>1* [div* s>>1**] ; rem too small )
        SWP2 OVR2 SUB2 SWP2      ( rem-y>>1* y>>1* [div* s>>1*] )
        DUP2r ROT2r ADD2r SWP2r  ( rem-y>>1* y>>1* [div+s>>1* s>>1*] )
        !&loop                   ( rem-y>>1* y>>1* [div+s>>1* s>>1*] )

    &done
    POP2 POP2 POP2r STH2r        ( div* )
    DUP2 #7fff GTH2 ?&oo JMP2r   ( div* )
    &o POP2 POP2 &oo POP2 #7fff JMP2r ( 7fff ; saturate on overflow )

@x16-signed-op ( x* y* f* -> f(x,y)* )
    STH2 LIT2r 0001
    DUP2 #8000 LTH2 ?&ypos x16-negate SWPr
    &ypos SWP2 DUP2 #8000 LTH2 ?&xpos x16-negate SWPr
    &xpos SWP2 SWP2r STH2r JSR2
    STHr ?&xypos x16-negate &xypos POPr
    JMP2r

@x16-quotient ( x* y* -> x//y* )
    ;x16-quot-unsigned !x16-signed-op

@x16-quot-unsigned ( x* y* -> x//y* )
    DIV2 DUP2 #007f GTH2 ?{ #80 SFT2 JMP2r }
    POP2 #7fff JMP2r

@x16-remainder ( x* y* -> x%y* )
    DIV2k MUL2 SUB2 JMP2r

@x16-from-s8 ( n^ -> x* )
    #00 JMP2r

@x16-from-s16 ( n* -> x* )
    DUP2 #ff80 GTH2 ?&neg
    DUP2 #007f GTH2 ?error
    NIP #00 SWP JMP2r
    &neg NIP #ff SWP JMP2r

( 1.5 -> 1, 0.5 -> 0, -1.5 -> -1 )
@x16-to-s16 ( x* -> whole* )
    DUP2 #7fff GTH2 ?&neg ( x0 x1 )
        DUP EOR   ( x0 00 )
        SWP JMP2r ( 0 x0 )
    &neg #00ff STHk      ( x0 x1 00 ff [ff] )
        ADD2 POP         ( x0' [ff] )
        DUP #07 SFT STHr ( x0' s' ff )
        MUL SWP JMP2r    ( 00-or-ff x0' )

( 1.5 -> 1, 0.5 -> 0, -1.5 -> -1 )
@x16-to-s8 ( x* -> whole^ )
    DUP2 #0f SFT2 #00ff MUL2 ADD2 POP JMP2r

( e.g. #0812 -> #0800 )
@x16-floor ( x* -> floor[x]* )
    DUP EOR JMP2r

( e.g. #0812 -> #0900 )
@x16-ceil ( x* -> floor[x]* )
    #00ff ADD2 DUP EOR JMP2r

( round half-even, #0080 -> #0000, #0180 -> #0200 )
@x16-round ( x* -> round[x]* )
    OVR #01 AND ?&odd
    ( even ) #007f !&rest
    &odd #0080
    &rest ADD2 DUP EOR JMP2r

( use up to 256 iterations of heron's algorithm )
( )
( NOTE: there are some inaccuracies here, currently )
( the algorithm doesn't always converge perfectly. )
( it should be "good enough" for many use cases. )
@x16-sqrt ( x* -> sqrt[x]* )
    LIT2r ff00
    LIT2r 0200                  ( [c* 2*] )
    DUP2 STH2kr x16-div         ( x* s=x/2* [c* 2*] )
    &loop                       ( x* s0* [c* 2*] )
        OVR2 OVR2 x16-div       ( x* s0* x/s0* [c* 2*] )
        OVR2 x16-add            ( x* s0* (x/s0)+s0* [c* 2*] )
        STH2kr x16-div          ( x* s0* s1=((x/s0)+s0)/2* [c* 2*] )
        SWP2 OVR2 x16-ne        ( x* s1* go^ [c* 2*] )
        SWP2r INC2r ORAkr STHr  ( x* s1* go^ ok^ [2* c+1*] )
        SWP2r ?&continue        ( x* s1* go^ [c+1* 2*] )
        POP !&done              ( x* s1* [c+1* 2*] )
    &continue                   ( x* s1* go^ [c+1* 2*] )
        ?&loop                  ( x* s1* [c+1* 2*] )
    &done                       ( x* s1* [c* 2*] )
        POP2r POP2r NIP2 JMP2r  ( s1* )

@x16-unit-circle ( x* -> x'* )
    x16-pi*2 STH2            ( x       [2pi] )
    DUP2 STH2kr x16-quotient ( x x/2pi [2pi] )
    DUP2 #1400 DIV2 STH2 SWP2r ( x x/2pi [adj* 2pi*] )
    STH2r x16-mul STH2r ADD2 SUB2       ( x' ; 0 <= x' < 2pi )
    JMP2r

@x16-cos ( x* -> cos[x]* )
    x16-unit-circle x16-pi/2 ADD2 ( fall-through )

@x16-sin ( x* -> sin[x]* )
    DUP2 #8000 LTH2 ?&non-negative
    x16-negate x16-sin/non-negative !x16-negate
    &non-negative
        x16-unit-circle
        DUP2 x16-3pi/2 LTH2 ?&c1
    ( -sin(2pi - x) ) x16-pi*2 SWP2 SUB2 x16-sin-q !x16-negate
    &c1 DUP2 x16-pi LTH2 ?&c2
    ( -sin(x - pi) )  x16-pi         SUB2 x16-sin-q !x16-negate
    &c2 DUP2 x16-pi/2 LTH2 ?&c3
    (  sin(pi - x) )  x16-pi   SWP2 SUB2 !x16-sin-q
    &c3
    (  sin[x] )                          ( fall-thru )

( 0 <= x < 2pi )
@x16-sin-q ( x* -> sin[x] )
    DUP2 ADD2 ;x16-sin-table ADD2 LDA2 JMP2r

@x16-sin-table
    0000 0001 0002 0003 0004 0005 0006 0007 0008 0009 000a 000b 000c 000d 000e 000f
    0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 001a 001b 001c 001d 001e 001f
    0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 002a 002b 002c 002d 002e 002f
    0030 0031 0032 0033 0034 0035 0036 0037 0038 0039 003a 003a 003b 003c 003d 003e
    003f 0040 0041 0042 0043 0044 0045 0046 0047 0048 0049 004a 004b 004c 004d 004e
    004f 0050 0051 0052 0053 0053 0054 0055 0056 0057 0058 0059 005a 005b 005c 005d
    005e 005f 0060 0061 0061 0062 0063 0064 0065 0066 0067 0068 0069 006a 006b 006c
    006c 006d 006e 006f 0070 0071 0072 0073 0074 0075 0075 0076 0077 0078 0079 007a
    007b 007c 007c 007d 007e 007f 0080 0081 0082 0083 0083 0084 0085 0086 0087 0088
    0089 0089 008a 008b 008c 008d 008e 008e 008f 0090 0091 0092 0093 0093 0094 0095
    0096 0097 0097 0098 0099 009a 009b 009b 009c 009d 009e 009f 009f 00a0 00a1 00a2
    00a2 00a3 00a4 00a5 00a6 00a6 00a7 00a8 00a9 00a9 00aa 00ab 00ac 00ac 00ad 00ae
    00ae 00af 00b0 00b1 00b1 00b2 00b3 00b4 00b4 00b5 00b6 00b6 00b7 00b8 00b8 00b9
    00ba 00bb 00bb 00bc 00bd 00bd 00be 00bf 00bf 00c0 00c1 00c1 00c2 00c3 00c3 00c4
    00c4 00c5 00c6 00c6 00c7 00c8 00c8 00c9 00ca 00ca 00cb 00cb 00cc 00cd 00cd 00ce
    00ce 00cf 00d0 00d0 00d1 00d1 00d2 00d2 00d3 00d4 00d4 00d5 00d5 00d6 00d6 00d7
    00d7 00d8 00d8 00d9 00da 00da 00db 00db 00dc 00dc 00dd 00dd 00de 00de 00df 00df
    00e0 00e0 00e1 00e1 00e2 00e2 00e2 00e3 00e3 00e4 00e4 00e5 00e5 00e6 00e6 00e7
    00e7 00e7 00e8 00e8 00e9 00e9 00ea 00ea 00ea 00eb 00eb 00ec 00ec 00ec 00ed 00ed
    00ed 00ee 00ee 00ef 00ef 00ef 00f0 00f0 00f0 00f1 00f1 00f1 00f2 00f2 00f2 00f3
    00f3 00f3 00f4 00f4 00f4 00f4 00f5 00f5 00f5 00f6 00f6 00f6 00f6 00f7 00f7 00f7
    00f8 00f8 00f8 00f8 00f8 00f9 00f9 00f9 00f9 00fa 00fa 00fa 00fa 00fb 00fb 00fb
    00fb 00fb 00fb 00fc 00fc 00fc 00fc 00fc 00fd 00fd 00fd 00fd 00fd 00fd 00fd 00fe
    00fe 00fe 00fe 00fe 00fe 00fe 00fe 00ff 00ff 00ff 00ff 00ff 00ff 00ff 00ff 00ff
    00ff 00ff 00ff 0100 0100 0100 0100 0100 0100 0100 0100 0100 0100 0100 0100 0100
    0100 0100 0100

@x16-tan ( x* -> tan[x]* )
    x16-pi*2 STH2            ( x       [2pi] )
    DUP2 STH2kr x16-quotient ( x x/2pi [2pi] )
    DUP2 #1400 DIV2 STH2 SWP2r ( x x/2pi [adj* 2pi*] )
    STH2r x16-mul STH2r ADD2 SUB2       ( x' ; 0 <= x' < 2pi )

    ( tan(pi/2) = tan(3pi/2) = error )
    DUP2 x16-3pi/2 EQU2 ?error
    DUP2 x16-pi/2 EQU2 ?error

        DUP2 x16-3pi/2 LTH2 ?&c1
    ( -tan(2pi - x) ) x16-pi*2 SWP2 SUB2 x16-tan-q !x16-negate
    &c1 DUP2 x16-pi LTH2 ?&c2
    (  tan(x - pi) )  x16-pi        SUB2 !x16-tan-q
    &c2 DUP2 x16-pi/2 LTH2 ?&c3
    ( -tan(pi - x) )  x16-pi   SWP2 SUB2 x16-tan-q !x16-negate
    &c3
    (  tan[x] )                          ( fall-thru )

( 0 <= x < 2pi )
@x16-tan-q ( x* -> sin[x] )
    DUP2 ADD2 ;x16-tan-table ADD2 LDA2 JMP2r

@x16-tan-table
    0000 0001 0002 0003 0004 0005 0006 0007 0008 0009 000a 000b 000c 000d 000e 000f
    0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 001a 001b 001c 001d 001e 001f
    0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 002a 002b 002c 002d 002f 0030
    0031 0032 0033 0034 0035 0036 0037 0038 0039 003a 003b 003c 003d 003e 003f 0040
    0041 0042 0044 0045 0046 0047 0048 0049 004a 004b 004c 004d 004e 004f 0051 0052
    0053 0054 0055 0056 0057 0058 0059 005b 005c 005d 005e 005f 0060 0061 0062 0064
    0065 0066 0067 0068 0069 006b 006c 006d 006e 006f 0071 0072 0073 0074 0075 0077
    0078 0079 007a 007b 007d 007e 007f 0080 0082 0083 0084 0085 0087 0088 0089 008b
    008c 008d 008e 0090 0091 0092 0094 0095 0096 0098 0099 009a 009c 009d 009f 00a0
    00a1 00a3 00a4 00a6 00a7 00a8 00aa 00ab 00ad 00ae 00b0 00b1 00b3 00b4 00b6 00b7
    00b9 00ba 00bc 00bd 00bf 00c0 00c2 00c4 00c5 00c7 00c8 00ca 00cc 00cd 00cf 00d1
    00d2 00d4 00d6 00d7 00d9 00db 00dc 00de 00e0 00e2 00e4 00e5 00e7 00e9 00eb 00ed
    00ee 00f0 00f2 00f4 00f6 00f8 00fa 00fc 00fe 0100 0102 0104 0106 0108 010a 010c
    010e 0110 0113 0115 0117 0119 011b 011e 0120 0122 0124 0127 0129 012b 012e 0130
    0133 0135 0137 013a 013c 013f 0142 0144 0147 0149 014c 014f 0152 0154 0157 015a
    015d 0160 0162 0165 0168 016b 016e 0171 0175 0178 017b 017e 0181 0185 0188 018b
    018f 0192 0196 0199 019d 01a0 01a4 01a8 01ac 01af 01b3 01b7 01bb 01bf 01c3 01c7
    01cc 01d0 01d4 01d8 01dd 01e1 01e6 01eb 01ef 01f4 01f9 01fe 0203 0208 020d 0212
    0218 021d 0223 0228 022e 0234 023a 0240 0246 024c 0252 0259 025f 0266 026d 0274
    027b 0282 0289 0291 0299 02a0 02a8 02b1 02b9 02c1 02ca 02d3 02dc 02e5 02ef 02f9
    0302 030d 0317 0322 032d 0338 0343 034f 035b 0368 0374 0382 038f 039d 03ab 03ba
    03c9 03d9 03e9 03f9 040a 041c 042e 0441 0454 0468 047d 0492 04a9 04c0 04d8 04f1
    050b 0526 0542 055f 057d 059d 05bf 05e1 0606 062c 0654 067e 06aa 06d9 070a 073e
    0775 07af 07ed 082f 0876 08c1 0911 0967 09c4 0a28 0a95 0b0a 0b8b 0c17 0cb2 0d5d
    0e1a 0eed 0fdb 10e8 121b 137d 1519 1700 1946 1c0c 1f80 23ed 29cc 31f5 3e13 51f2
    7888 7fff 7fff

@x16-log ( x* -> log[x]* )

    [ DUP2 #0000 GTH2 ] STH
    [ DUP2 #8000 LTH2 ] STHr AND ?&0<x<128
    ( error ) !error

    &0<x<128     DUP2 #0800 GTH2 ?&8<x<128
    ( 0<x<=8 )   DUP2 #0200 GTH2 ?&2<x<=8
    ( 0<x<=2 )                         !x16-log-q
    &2<x<=8      DUP2 #0400 GTH2 ?&4<x<=8
    ( 2<x<=4 )     #0200 x16-div x16-log-q #00b1 ADD2 JMP2r
    &4<x<=8        #0400 x16-div x16-log-q #0163 ADD2 JMP2r
    &8<x<128     DUP2 #2000 GTH2 ?&32<x<128
    ( 8<x<=32 )  DUP2 #1000 GTH2 ?&16<x<=32
    ( 8<x<=16 )    #0800 x16-div x16-log-q #0214 ADD2 JMP2r
    &16<x<=32      #1000 x16-div x16-log-q #02c6 ADD2 JMP2r
    &32<x<128    DUP2 #4000 GTH2 ?&64<x<128
    ( 32<x<=64 )   #2000 x16-div x16-log-q #0377 ADD2 JMP2r
    &64<x<128      #4000 x16-div x16-log-q #0429 ADD2 JMP2r

( 0 < x <= 2 )
@x16-log-q ( x* -> log[x]* )
    DUP2 ADD2 ;x16-log-table ADD2 LDA2 JMP2r

( the first entry, i.e. log[0], is invalid and should not be used. )
( the last entry is log[2]. )
@x16-log-table
    8000 fa74 fb26 fb8e fbd7 fc10 fc3f fc67 fc89 fca7 fcc2 fcda fcf1 fd05 fd18 fd2a
    fd3a fd4a fd58 fd66 fd73 fd80 fd8c fd97 fda2 fdac fdb7 fdc0 fdc9 fdd2 fddb fde4
    fdec fdf4 fdfb fe03 fe0a fe11 fe18 fe1e fe25 fe2b fe31 fe37 fe3d fe43 fe49 fe4e
    fe53 fe59 fe5e fe63 fe68 fe6d fe72 fe76 fe7b fe7f fe84 fe88 fe8d fe91 fe95 fe99
    fe9d fea1 fea5 fea9 fead feb0 feb4 feb8 febb febf fec2 fec6 fec9 fecc fed0 fed3
    fed6 fed9 fedd fee0 fee3 fee6 fee9 feec feef fef2 fef4 fef7 fefa fefd ff00 ff02
    ff05 ff08 ff0a ff0d ff0f ff12 ff14 ff17 ff19 ff1c ff1e ff21 ff23 ff25 ff28 ff2a
    ff2c ff2f ff31 ff33 ff35 ff38 ff3a ff3c ff3e ff40 ff42 ff44 ff46 ff48 ff4b ff4d
    ff4f ff51 ff53 ff54 ff56 ff58 ff5a ff5c ff5e ff60 ff62 ff64 ff65 ff67 ff69 ff6b
    ff6d ff6e ff70 ff72 ff74 ff75 ff77 ff79 ff7b ff7c ff7e ff80 ff81 ff83 ff84 ff86
    ff88 ff89 ff8b ff8c ff8e ff90 ff91 ff93 ff94 ff96 ff97 ff99 ff9a ff9c ff9d ff9f
    ffa0 ffa2 ffa3 ffa4 ffa6 ffa7 ffa9 ffaa ffab ffad ffae ffb0 ffb1 ffb2 ffb4 ffb5
    ffb6 ffb8 ffb9 ffba ffbc ffbd ffbe ffc0 ffc1 ffc2 ffc3 ffc5 ffc6 ffc7 ffc8 ffca
    ffcb ffcc ffcd ffcf ffd0 ffd1 ffd2 ffd3 ffd5 ffd6 ffd7 ffd8 ffd9 ffda ffdc ffdd
    ffde ffdf ffe0 ffe1 ffe2 ffe3 ffe5 ffe6 ffe7 ffe8 ffe9 ffea ffeb ffec ffed ffee
    ffef fff1 fff2 fff3 fff4 fff5 fff6 fff7 fff8 fff9 fffa fffb fffc fffd fffe ffff
    0000 0001 0002 0003 0004 0005 0006 0007 0008 0009 000a 000b 000c 000d 000e 000f
    0010 0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 001a 001b 001b 001c 001d
    001e 001f 0020 0021 0022 0023 0023 0024 0025 0026 0027 0028 0029 0029 002a 002b
    002c 002d 002e 002f 002f 0030 0031 0032 0033 0033 0034 0035 0036 0037 0038 0038
    0039 003a 003b 003c 003c 003d 003e 003f 003f 0040 0041 0042 0043 0043 0044 0045
    0046 0046 0047 0048 0049 0049 004a 004b 004c 004c 004d 004e 004f 004f 0050 0051
    0052 0052 0053 0054 0054 0055 0056 0057 0057 0058 0059 0059 005a 005b 005c 005c
    005d 005e 005e 005f 0060 0060 0061 0062 0062 0063 0064 0064 0065 0066 0066 0067
    0068 0068 0069 006a 006a 006b 006c 006c 006d 006e 006e 006f 0070 0070 0071 0072
    0072 0073 0074 0074 0075 0075 0076 0077 0077 0078 0079 0079 007a 007a 007b 007c
    007c 007d 007e 007e 007f 007f 0080 0081 0081 0082 0082 0083 0084 0084 0085 0085
    0086 0087 0087 0088 0088 0089 0089 008a 008b 008b 008c 008c 008d 008e 008e 008f
    008f 0090 0090 0091 0092 0092 0093 0093 0094 0094 0095 0095 0096 0097 0097 0098
    0098 0099 0099 009a 009a 009b 009c 009c 009d 009d 009e 009e 009f 009f 00a0 00a0
    00a1 00a1 00a2 00a3 00a3 00a4 00a4 00a5 00a5 00a6 00a6 00a7 00a7 00a8 00a8 00a9
    00a9 00aa 00aa 00ab 00ab 00ac 00ac 00ad 00ad 00ae 00ae 00af 00af 00b0 00b0 00b1
    00b1