# ray.sed
# Mike Wilson, cmikewilsonAgmailDcom
# April 25, 2005 - June 3, 2005

# [FIRST LINE]
s/^#.*//
s/$/,/
G
s/\n//
/^,/d
s/,$/,@/

:cycle
t cycle    # clear sed substitution flag
#p
# print top function: cat | sort -r | uniq -c
#h
#s/^\([^, ]*\).*/\1/
#w c1
#:x
#tx
#g
/^A/{
    /^ADD /b add
    /^ABS /b abs
    /^ASSERT /b assert
    /^ASSERT2 /b assert2
    b nomatch
}
/^C/{
    /^COMPARE /b compare
    /^CROSS /b cross
    /^CONS /b cons
    b nomatch
}
/^D/{
    /^DI/{
        /^DIVREM2 /b divrem2
        /^DIVREM3 /b divrem3
        /^DIVREM /b divrem
        /^DIV /b div
        b nomatch
    }
    /^DECR /b decr
    /^DEFAULT /b default
    /^DTOB /b dtob
    /^DELETE /b delete
    b nomatch
}
/^F/{
    /^FMULT /b fmult
    /^FADD /b fadd
    /^FSUB /b fsub
    /^FAPPROX /b fapprox
    /^FCMP2 /b fcmp2
    /^FCMP /b fcmp
    /^FD/{
        /^FDIV4 /b fdiv4
        /^FDIV2 /b fdiv2
        /^FDIV3 /b fdiv3

        /^FDIV /b fdiv
        /^FDTOB /b fdtob
        /^FDTOB2 /b fdtob2
        b nomatch
    }
    /^FBTOD /b fbtod
    /^FBTOD2 /b fbtod2
    /^FTRIM /b ftrim
    b nomatch
}
/^I/{
    /^IADD /b iadd
    /^IMULT /b imult
    /^ISUB /b isub
    /^IN/{
        /^INTERSECTD /b intersectd
        /^INTERSECTC /b intersectc
        /^INTERSECTB /b intersectb
        /^INTERSECT /b intersect
        /^INTERSECTA /b intersecta
        /^INCR /b incr
        b nomatch
    }
    /^ICMP /b icmp
    /^IB/{
        /^IBTOD3 /b ibtod3
        /^IBTOD4 /b ibtod4
        /^IBTOD2 /b ibtod2
        /^IBTOD5 /b ibtod5
        /^IBTOD /b ibtod
        b nomatch
    }
    /^ICMP2 /b icmp2
    /^IDTOB /b idtob
    /^IDTOB2 /b idtob2
    /^IDIV /b idiv
    b nomatch
}
/^O/{
    /^OUTPUTLINE /b outputline
    /^OUTPUTPRE /b outputpre
    /^OUTPUTPOST /b outputpost
    /^OUTPUTPS1,/b outputps1
    /^OUTPUTPS2 /b outputps2
    /^OUTPUTPS3,/b outputps3
    /^OUTPUTPS4,/b outputps4
    b nomatch
}
/^P/{
    /^POW10 /b pow10
    /^PRINTOBJ /b printobj
    /^PRINT /b print
    /^PRINTSTACK,/b printstack
    b nomatch
}
/^R/{
    /^RE/{
        /^REALIZE1 /b realize1
        /^REALIZE2 /b realize2
        /^RETURN /b return
        /^REF /b ref
        /^REM /b rem
        /^READCMD,/b readcmd
        b nomatch
    }
    /^RA/{
        /^RAYTRACEHIT /b raytracehit
        /^RAYTRACE2 /b raytrace2
        /^RAYTRACEPI /b raytracepi
        /^RAYTRACEADDSPHERE /b raytraceaddsphere
        /^RAYTRACEADDLIGHT /b raytraceaddlight
        /^RAYTRACELOOKFROMPOINT /b raytracelookfrompoint
        /^RAYTRACELOOKATPOINT /b raytracelookatpoint
        /^RAYTRACEWINDOWWIDTH /b raytracewindowwidth
        /^RAYTRACEWINDOWHEIGHT /b raytracewindowheight
        /^RAYTRACEWINDOWWIDTHRESOLUTION /b raytracewindowwidthresolution
        /^RAYTRACEWINDOWHEIGHTRESOLUTION /b raytracewindowheightresolution
        /^RAYTRACEWINDOWDISTANCE /b raytracewindowdistance
        /^RAYTRACEOUTPUTTYPE /b raytraceoutputtype
        /^RAYTRACEUP /b raytraceup
        /^RAYTRACE,/b raytrace
        b nomatch
    }
    b nomatch
}
/^S/{
    /^SUBST /b subst
    /^SQUARE /b square
    /^SELECT /b select
    /^SCALARMULT /b scalarmult
    /^SQRT2 /b sqrt2
    /^SET /b set
    /^SQRT /b sqrt
    /^SUB /b sub
    b nomatch
}
/^UNPAD /b unpad
/^NEG /b neg
/^VADD /b vadd
/^BTOD /b btod
/^VSUB /b vsub
/^LREF /b lref
/^MULT /b mult
/^UNIT /b unit
/^MAG /b mag
/^LCHANGE /b lchange
/^FLUID /b fluid
/^ERROR /b error
/^EXIT[ ,]/b exit


:nomatch
# [AFTER JUMPS]
s/^@/READCMD,@/
t readcmd
s/^\([^,]\+\),.*/ERROR Unrecognized Command: \1,/
t error
b exit


##################################
## RAYTRACE INPUT FILE COMMANDS ##
##################################


:raytracelookfrompoint
s/^RAYTRACELOOKFROMPOINT \(-\?[[:digit:]]\+\) \(-\?[[:digit:]]\+\) \(-\?[[:digit:]]\+\),/DTOB \1,SUBST X _,DTOB \2,SUBST Y _,DTOB \3,SET LookFromPoint V#{X}#{Y}#_#,!,!,/
t dtob
b loop-error
:raytracelookatpoint
s/^RAYTRACELOOKATPOINT \(-\?[[:digit:]]\+\) \(-\?[[:digit:]]\+\) \(-\?[[:digit:]]\+\),/DTOB \1,SUBST X _,DTOB \2,SUBST Y _,DTOB \3,SET LookAtPoint V#{X}#{Y}#_#,!,!,/
t dtob
b loop-error
:raytraceup
# window up, a unit vector pointing up.
s/^RAYTRACEUP \(-\?[01]\) \(-\?[01]\) \(-\?[01]\),/SET Up V#B\1#B\2#B\3#,/
t set
b loop-error
:raytracewindowwidth
# raytrace window width
s/^RAYTRACEWINDOWWIDTH \(-\?[[:digit:]]\+\),/DTOB \1,SET WindowWidth _,/
t dtob
b loop-error
:raytracewindowheight
# raytrace window height
s/^RAYTRACEWINDOWHEIGHT \(-\?[[:digit:]]\+\),/DTOB \1,SET WindowHeight _,/
t dtob
b loop-error
:raytracewindowwidthresolution
s/^RAYTRACEWINDOWWIDTHRESOLUTION \(-\?[[:digit:]]\+\),/DTOB \1,SET WindowWidthResolution _,/
t dtob
b loop-error
:raytracewindowheightresolution
s/^RAYTRACEWINDOWHEIGHTRESOLUTION \(-\?[[:digit:]]\+\),/DTOB \1,SET WindowHeightResolution _,/
t dtob
b loop-error
:raytracewindowdistance
s/^RAYTRACEWINDOWDISTANCE \(-\?[[:digit:]]\+\),/DTOB \1,SET WindowDistance _,/
t dtob
b loop-error
:raytraceaddsphere
s/^RAYTRACEADDSPHERE \(-\?[[:digit:]]\+\) \(-\?[[:digit:]]\+\) \(-\?[[:digit:]]\+\) \(-\?[[:digit:]]\+\) \([0-9A-F]\{6\}\),/DTOB \1,SUBST X _,DTOB \2,SUBST Y _,DTOB \3,SUBST Z _,DTOB \4,SUBST R _,DEFAULT Spheres <>,REF Spheres,CONS V#{X}#{Y}#{Z}#{R}#\5# _,!,!,!,!,SET Spheres _,/
t dtob
b loop-error
:raytraceaddlight
s/^RAYTRACEADDLIGHT \(-\?[[:digit:]]\+\) \(-\?[[:digit:]]\+\) \(-\?[[:digit:]]\+\),/DTOB \1,SUBST X _,DTOB \2,SUBST Y _,DTOB \3,SUBST Z _,DEFAULT Lights <>,REF Lights,CONS V#{X}#{Y}#{Z}# _,!,!,!,SET Lights _,/
t dtob
b loop-error
:raytraceoutputtype
s/^RAYTRACEOUTPUTTYPE \(PS\|SCREEN\),/SET OutputType \1,/
t set
b loop-error


##############
## RAYTRACE ##
##############


:raytrace    # [RAYTRACE Algorithm]
# A = LAP-LFP
# Au = unit(A)
# WP = LFP + wd * Au
# B = A X UP
# RT = unit(B)
# pdw = WindowWidth / WindowWidthResolution
# DW = pdw * RT
# pdh = WindowHeight / WindowHeightResolution
# DH = pdh * UP
# Pi = WP + w*DW + h*DH, w = -WWR/2 .. WWR/2, h = -WHR/2 .. WHR/2
s/^RAYTRACE,/DEFAULT LookFromPoint V#B0#B0#B100000#,DEFAULT LookAtPoint V#B0#B0#B0#,DEFAULT Up V#B0#B1#B0#,DEFAULT WindowDistance B10000,DEFAULT WindowWidth B1000,DEFAULT WindowHeight B1000,DEFAULT WindowWidthResolution 20,DEFAULT WindowHeightResolution 20,DEFAULT Spheres <V#B0#B0#B0#B100#000000#:>,DEFAULT Lights <V#B0#B1000#B0#:>,DEFAULT OutputType SCREEN,REF LookFromPoint,PRINTOBJ LFP _,SUBST LFP _,REF LookAtPoint,PRINTOBJ LAP _,VSUB _ {LFP},PRINTOBJ A _,SUBST A _,UNIT {A},PRINTOBJ unitTAunsimp _,PRINTOBJ Au _,SUBST Au _,REF WindowDistance,PRINTOBJ wd _,SCALARMULT _ {Au},VADD {LFP} _,!,PRINTOBJ WP _,SUBST WP _,REF Up,PRINTOBJ UP _,SUBST UP _,CROSS {A} {UP},UNIT _,PRINTOBJ RT _,SUBST RT _,REF WindowWidth,PRINTOBJ ww _,SUBST ww _,REF WindowWidthResolution,PRINTOBJ wwr _,SUBST wwr _,FDIV {ww}.00000000 {wwr}.00000000,SCALARMULT _ {RT},PRINTOBJ DW _,SUBST DW _,REF WindowHeight,PRINTOBJ wh _,SUBST wh _,REF WindowHeightResolution,PRINTOBJ whr _,SUBST whr _,FDIV {wh}.00000000 {whr}.00000000,SCALARMULT _ {UP},PRINTOBJ DH _,SUBST DH _,DIV {wwr} B10,NEG _,PRINTOBJ wmin _,SUBST wmin _,ADD {wmin} {wwr},PRINTOBJ wmax _,SUBST wmax _,DIV {whr} B10,NEG _,PRINTOBJ hmin _,SUBST hmin _,ADD {hmin} {whr},PRINTOBJ hmax _,SUBST hmax _,REF OutputType,SUBST ot _,OUTPUTPRE {ot} {wwr} {whr},RAYTRACE2 {WP} {DW} {DH} {wmin} {wmin} {wmax} {hmin} {hmax} {hmax} <>,OUTPUTPOST {ot},!,!,!,!,!,!,!,!,!,!,!,!,!,!,!,!,/
t default
b loop-error
:raytrace2
# RAYTRACE2 WP DW DH wmin w wmax hmin h hmax <color-string>
# Traverse each pixel in the window, and fire a ray through
# that pixel and receive a color value.
# RAYTRACEPI calculates the pixel location (Pi) from integer w,h values,
# and RAYTRACEHIT fires the ray and returns the color value.
# Want to start at wmin,hmax and work right and down
# Pi = WP + w*DW + h*DH, w+=1, h+=1
# while h>hmin
#   while w<wmax
#     for H in Hs
#       if intersect(Po Pi H): r=color
#     r=black
#     append r to color-string
#     w=w+1
#   h=h-1
#   w=wmax
# return <color-string>
#
# RAYTRACE2 WP DW DH wmin w wmax hmin h hmax <color-string>
# if h=hmin, finished
s/^RAYTRACE2 V#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+# V#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+# V#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+# B-\?[01]\+ B-\?[01]\+ B-\?[01]\+ B\(-\?[01]\+\) B\1 B-\?[01]\+ <\([^>]\+\)>,//
t cycle
# if w=wmax, output line, w=wmin, h=h-1
s/^RAYTRACE2 V\(#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#\) V\(#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#\) V\(#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#\) B\(-\?[01]\+\) B\(-\?[01]\+\) B\5 B\(-\?[01]\+\) B\(-\?[01]\+\) B\(-\?[01]\+\) <\([^>]*\)>,/OUTPUTLINE <\9>,DECR B\7,RAYTRACE2 V\1 V\2 V\3 B\4 B\4 B\5 B\6 _ B\8 <>,/
t outputline
# otherwise, w=w+1, calc
s/^RAYTRACE2 \(V#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+# V#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+# V#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#\) B\(-\?[01]\+\) B\(-\?[01]\+\) B\(-\?[01]\+\) B\(-\?[01]\+\) B\(-\?[01]\+\) B\(-\?[01]\+\) <\([^>]*\)>,/RAYTRACEPI \1 B\3 B\6,SUBST Pi _,REF LookFromPoint,SUBST LFP _,REF Spheres,SUBST Hs _,INCR B\3,SUBST Wprime _,RAYTRACEHIT {LFP} {Pi} B\3 B\4 {Hs},RAYTRACE2 \1 B\2 {Wprime} B\4 B\5 B\6 B\7 <\8_>,!,!,!,!,/
t raytracepi
b loop-error

:raytracepi
# RAYTRACEPI WP DW DH w h => Pi
# Pi = WP + w*DW + h*DH
s/^RAYTRACEPI V\(#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#\) V\(#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#\) V\(#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#\) B\(-\?[01]\+\) B\(-\?[01]\+\),/SCALARMULT B\4 V\2,SUBST T _,SCALARMULT B\5 V\3,VADD {T} _,VADD V\1 _,RETURN _,!,/
t scalarmult
b loop-error


# We really need to calculate an intersection with each sphere
# and receive a distance, sphere pair.  Then find the sphere with
# the smallest distance.  That's the one we hit.
# Then we need to calculate a normal at the point of impact,
# calculate a reflected ray, and recurse.

:raytracehit
# RAYTRACEHIT LFP Pi w wmax Hs => color
# Hs==<>, must not have hit anything, return background color
s/^RAYTRACEHIT V#B-\?[01]\+#B-\?[01]\+#B-\?[01]\+# V#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+# B-\?[01]\+ B-\?[01]\+ <>,/RETURN 000000,/
t return
# if intersect with sphere, append color and recurse
# INTERSECT V1 V2 H1
s/^RAYTRACEHIT V\(#B-\?[01]\+#B-\?[01]\+#B-\?[01]\+#\) V\(#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#\) B\(-\?[01]\+\) B\(-\?[01]\+\) <V\(#B-\?[01]\+#B-\?[01]\+#B-\?[01]\+#B[01]\+#\([0-9A-F]\{6\}\)#\):\([^>]*\)>,/INTERSECT V\1 V\2 V\5,SELECT _ (SNONE RAYTRACEHIT V\1 V\2 B\3 B\4 <\7>) (STANGENT RETURN \6) (STWICE RETURN \6),_,/
t intersect
b loop-error


########################
## RAYTRACE INTERSECT ##
########################


# INTERSECT P1 P2 H3 => SNONE,STWICE,STANGENT [INTERSECT algorithm]
:intersect
# does the line P1->P2 intersect the sphere H3
# Points look like: P'x'y'z with x,y,z in binary form (B0100, ...)
# returns SNONE, STANGENT, STWICE
# (Alg. from http://astronomy.swin.edu.au/~pbourke/geometry/sphereline/)
# a = (x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2
# b = 2[(x2-x1)(x1-x3)+(y2-y1)(y1-y3)+(z2-z1)(z1-z3)]
# c = x3^2 + y3^2 + z3^2 + x1^2 + y1^2 + z1^2 - 2[x3x1 + y3y1 + z3z1] - r^2
# d = b^2 - 4ac {
#     d<0 no intersection; d=0 tangent at u=-b/2a; d>0 intersect at 2 pts
# P = P1 + u(P2-P1)
s/^INTERSECT V\(#B-\?[01]\+#B-\?[01]\+#B-\?[01]\+#\) V\(#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#B-\?[01]\+\.[01]\+#\) V\(#B-\?[01]\+#B-\?[01]\+#B-\?[01]\+#\)B\(-\?[01]\+\)#[0-9A-F]\{6\}#,/INTERSECTA V\1 V\2,SUBST A _,INTERSECTB V\1 V\2 V\3,SUBST B _,INTERSECTC V\1 V\3 B\4,INTERSECTD {A} {B} _,!,!,FCMP _ B0.00000000 (SNONE) (STWICE) (STANGENT),/
t intersecta
b loop-error
:intersecta
# INTERSECTA P1 P2
# a = (x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2
s/^INTERSECTA V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)# V#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#,/FSUB B\4 B\1.00000000,SQUARE _,SUBST T1 _,FSUB B\5 B\2.00000000,SQUARE _,FADD {T1} _,!,SUBST T2 _,FSUB B\6 B\3.00000000,SQUARE _,FADD {T2} _,!,/
t fsub
b loop-error
:intersectb
# INTERSECTB P1 P2 P3
# b = 2[(x2-x1)(x1-x3)+(y2-y1)(y1-y3)+(z2-z1)(z1-z3)]
s/^INTERSECTB V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)# V#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)# V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)#,/FSUB B\4 B\1.00000000,SUBST T1 _,ISUB B\1 B\7,FMULT {T1} _.00000000,!,SUBST T2 _,FSUB B\5 B\2.00000000,SUBST T3 _,ISUB B\2 B\8,FMULT {T3} _.00000000,FADD {T2} _,!,!,SUBST T4 _,FSUB B\6 B\3.00000000,SUBST T5 _,ISUB B\3 B\9,FMULT {T5} _.00000000,FADD {T4} _,!,!,FMULT B10.00000000 _,/
t fsub
b loop-error
:intersectc
# INTERSECTC P1 P3 R
# c = x3^2 + y3^2 + z3^2 + x1^2 + y1^2 + z1^2 - 2[x3x1 + y3y1 + z3z1] - r^2
s/INTERSECTC V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)# V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)# B\(-\?[01]\+\),/IMULT B\4 B\4,SUBST T1 _,IMULT B\5 B\5,IADD {T1} _,!,SUBST T2 _,IMULT B\6 B\6,IADD {T2} _,!,SUBST T3 _,IMULT B\1 B\1,IADD {T3} _,!,SUBST T4 _,IMULT B\2 B\2,IADD {T4} _,!,SUBST T5 _,IMULT B\3 B\3,IADD {T5} _,!,SUBST C1 _,IMULT B\4 B\1,SUBST T6 _,IMULT B\5 B\2,IADD {T6} _,SUBST T7 _,IMULT B\6 B\3,IADD {T7} _,IMULT B10 _,ISUB {C1} _,!,!,!,SUBST T8 _,IMULT B\7 B\7,ISUB {T8} _,!,/
t imult
b loop-error
:intersectd
# INTERSECTD A B C
# d = b^2 - 4ac
s/INTERSECTD B\(-\?[01]\+\.[01]\+\) B\(-\?[01]\+\.[01]\+\) B\(-\?[01]\+\),/FMULT B\2 B\2,SUBST T1 _,FMULT B100.00000000 B\1,FMULT B\3.00000000 _,FSUB {T1} _,!,/
t fmult
b loop-error


###################################################################
## RAYTRACE OUTPUT: OUTPUTPRE, OUTPUTLINE, OUTPUTPOST, OUTPUTPS* ##
###################################################################


:outputpre
# OUTPUTPRE PS/SCREEN wwr whr
s/^OUTPUTPRE SCREEN B[01]\+ B[01]\+,//
t cycle
s/^OUTPUTPRE PS B\([01]\+\) B\([01]\+\),/OUTPUTPS1,BTOD B\1,SUBST wwrd _,BTOD B\2,SUBST whrd _,OUTPUTPS2 {wwrd} {whrd},OUTPUTPS3,!,!,/
t outputps1
b loop-error


:outputline
s/^OUTPUTLINE <\([^>]*\)>,/PRINT \1,/
t print
b loop-error


:outputpost
s/^OUTPUTPOST SCREEN,//
t cycle
s/^OUTPUTPOST PS,/OUTPUTPS4,/
t outputps4
b loop-error


:outputps1
s/^OUTPUTPS1,//
i\
%!PS-Adobe-3.0\
%%Title: (ray.sed output)\
%%Creator: (Mike Wilson)\
%%Copyright: (Anyone may copy and/or modify this document without restriction.)\
%%LanguageLevel: 2\
%%Orientation: Portrait\
%%DocumentMedia: Letter 612 792 0 () ()\
%%Pages: 1\
%%Version: 1.0 1\
%%EndComments\
%%BeginDefaults\
%%PageOrientation: Portrait\
%%EndDefaults\
%%BeginSetup\
/inch { 72 mul } bind def\
/PageWidth 8.5 inch def\
/PageHeight 11 inch def
t cycle
b loop-error

:outputps2
# OUTPUTPS2 IMAGEX IMAGEY
h
s/^OUTPUTPS2 \([[:digit:]]\+\) \([[:digit:]]\+\),.*/\/ImageX \1 def\n\/ImageY \2 def/
p
g
s/^OUTPUTPS2 [[:digit:]]\+ [[:digit:]]\+,//
t cycle
b loop-error

:outputps3
s/^OUTPUTPS3,//
i\
/ImageWidth 6 inch def\
/ImageHeight ImageWidth ImageX div ImageY mul def\
%%EndSetup\
%%Page: (1) 1\
%%BeginPageSetup\
    /pgsave save def\
%%EndPageSetup\
    /picstr ImageX 3 mul string def\
    /DeviceRGB setcolorspace\
    PageWidth ImageWidth sub 2 div\
    PageHeight ImageHeight sub 2 div\
    translate\
    ImageWidth dup scale\
    ImageX ImageY 8\
    [ImageX 0 0 ImageY neg 0 ImageY]    % Sx - - Sy Tx Ty\
    { currentfile picstr readhexstring pop }\
    false 3\
    colorimage
t cycle
b loop-error

:outputps4
s/^OUTPUTPS4,//
i\
    pgsave restore\
    showpage\
%%PageTrailer\
%%Trailer\
%%EOF
t cycle
b loop-error


################################


:compare
# COMPARE A B <R1:R2:R3:>
# If A<B, return R1; if A==B, return R2; if A>B, return R3
# Is all this complexity really quicker than SUB A B, and sign check?
#   Yes, it's much faster -- like 10x.
# negative numbers?
/^COMPARE B\(-\|[01]\+ B-\)/{
    s/^COMPARE B-[01]\+ B[01]\+ <\([^:]\+\):[^:]\+:[^:]\+:>,/RETURN \1,/
    t return
    s/^COMPARE B[01]\+ B-[01]\+ <[^:]\+:[^:]\+:\([^:]\+\):>,/RETURN \1,/
    t return
    s/^COMPARE B-\([01]\+\) B-\([01]\+\) /COMPARE B\2 B\1 /
    t compare+
    b loop-error
}
:compare+
#
# Are both numbers equal?
s/^COMPARE B\([01]\+\) B\1 <[^:]\+:\([^:]\+\):[^:]\+:>,/RETURN \2,/
t return
#
# Is one number shorter than the other?
# count digits
s/^COMPARE B\([01]\+\) B\([01]\+\) <\([^:]\+\):[^:]\+:\([^:]\+\):>,/COMPC B\1 B\2 B\1 B\2 <\3:\4:>,/
t compc
# COMPARE A B <R1:R2:R3:>  =>  COMPC A B A B <R1:R3:>
/^COMPC /{
    :compc
    s/^COMPC B\([01]*\)[01] B\([01]*\)[01] B\([01]\+\) B\([01]\+\) /COMPC B\1 B\2 B\3 B\4 /
    t compc
    s/^COMPC B B B\([01]\+\) B\([01]\+\) /COMPSL B\1 B\2 /
    t compsl
    /^COMPSL /{
        # COMPSL A B <R1:R3:>
        # Both numbers are the same length; compare digits L to R
        :compsl
        s/^COMPSL B0[01]* B1[01]* <\([^:]\+\):[^:]\+:>,/RETURN \1,/
        t return
        s/^COMPSL B1[01]* B0[01]* <[^:]\+:\([^:]\+\):>,/RETURN \1,/
        t return
        s/^COMPSL B[01]\([01]\+\) B[01]\([01]\+\) /COMPSL B\1 B\2 /
        t compsl
        b loop-error
    }
    s/^COMPC B B[01]\+ B[01]\+ B[01]\+ <\([^:]\+\):[^:]\+:>,/RETURN \1,/
    t return
    s/^COMPC B[01]\+ B B[01]\+ B[01]\+ <[^:]\+:\([^:]\+\):>,/RETURN \1,/
    t return
}
b loop-error


######################################################
## GENERAL UTILITIES: RETURN, SUBST, SELECT, UNPAD  ##
######################################################


:return
# We have a value, substitute into the caller (where the _ is)
# Value may be anything, but it may not contain a comma.
# RETURN is a hotspot
# Check for a fluid frame in the way and get rid of it.
#s/^RETURN \([^,]\+\),FLUID [^,]\+,/RETURN \1,/
#t return
# If we return at the end of the stack, print it.
##s/^RETURN \([^,]\+\),@/PRINT \1,@/
##t print
# Replace the _.
s/^RETURN \([^,]\+\),\([^,_]*\)_/\2\1/
t cycle
# If the recipient doesn't contain a _, throw the value away.
s/^RETURN [^,]\+,//
t cycle
b loop-error


:subst
# SUBST id value,...,!,
# substitute a value for {id} in the stack
# only affects up to the first ! frame
# SUBST is a hotspot
s/^SUBST \([[:alnum:]]\+\) \([^,]\+\),\([^!]*\){\1}\([^!]*\),!,/SUBST \1 \2,\3\2\4,!,/
t subst
s/^SUBST [[:alnum:]]\+ \([^,]\+\),\([^!]\+\),!,/\2,/
t cycle
s/^SUBST [[:alnum:]]\+ \([^,]\+\),!,//
t cycle
b loop-error


:select
# SELECT A (X value) (Y value) (A value),
# return value associated with A; default to last
# SELECT is a hotspot
s/^SELECT \([^, ]\+\) \(([^,)]*) \)*(\1[ ]\+\([^,)]\+\))\( ([^,)]*)\)*,/RETURN \3,/
t return
# no match, use last as default; might or might not have label
s/^SELECT [^, ]\+ \(([^,)]*) \)*(\([^, ]\+[ ]\+\)\?\([^,)]*\)),/RETURN \3,/
t return
b loop-error


:unpad
# remove leading zeroes
# might be base 2 or 10; might have fractional part
# UNPAD is a hotspot
# binary numbers
s/^UNPAD B\(-\?\)0*\([01]\+\),/RETURN B\1\2,/
t return
# binary fractions
s/^UNPAD B\(-\?\)0*\([01]\+\.[01]\{8\}\),/RETURN B\1\2,/
t return
# base 10 numbers
s/^UNPAD \(-\?\)0*\([[:digit:]]\+\),/RETURN \1\2,/
t return
# base 10 fractions have 8 digits?
s/^UNPAD \(-\?\)0*\([[:digit:]]\+\.[[:digit:]]\+\),/RETURN \1\2,/
t return
b loop-error


###########################################
## GENERIC ABS, ADD, SUB, MULT, DIV, NEG ##
###########################################

:abs
s/^ABS B-\([01]\+\(\.[01]\{8\}\)\?\),/RETURN B\1,/
t return
s/^ABS B\([01]\+\(\.[01]\{8\}\)\?\),/RETURN B\1,/
t return
b loop-error


:add
s/^ADD B\(-\?[01]\+\) B\(-\?[01]\+\),/IADD B\1 B\2,/
t iadd
b loop-error


:sub
s/^SUB B\(-\?[01]\+\) B\(-\?[01]\+\),/ISUB B\1 B\2,/
t isub
b loop-error


:mult
s/^MULT B\(-\?[01]\+\) B\(-\?[01]\+\),/IMULT B\1 B\2,/
t imult
b loop-error

:div
s/^DIV B\(-\?[01]\+\) B\(-\?[01]\+\),/IDIV B\1 B\2,/
t idiv
b loop-error


:neg    # negate
s/^NEG B0,/RETURN B0,/
t return
s/^NEG B\([01]\+\(\.[01]\+\)\?\),/RETURN B-\1,/
t return
s/^NEG B-\([01]\+\(\.[01]\+\)\?\),/RETURN B\1,/
t return
# base 10
s/^NEG 0,/RETURN 0,/
s/^NEG -\([0-9]\+\),/RETURN \1,/
s/^NEG \([0-9]\+\),/RETURN -\1,/
t return
b loop-error


:square
s/^SQUARE B\(-\?[01]\+\),/IMULT B\1 B\1,/
t imult
s/^SQUARE B\(-\?[01]\+.[01]\+\),/FMULT B\1 B\1,/
t fmult
b loop-error


#######################################################################
## INTEGER INCR, DECR, IABS, ICMP, ADD, SUB, MULT, DIVREM, IDIV, REM ##
#######################################################################


:incr    # increment
s/^INCR B\([01]*\)0,/RETURN B\11,/
t return
s/^INCR B1,/RETURN B10,/
t return
s/^INCR B\([01]\+\)1,/INCR B\1,RETURN _0,/
t incr
s/^INCR B-\([01]\+\),/DECR B\1,NEG _,/
t decr
b loop-error


:decr    # decrement
s/^DECR B0,/RETURN B-1,/
t return
s/^DECR B\([01]*\)1,/RETURN B\10,/
t return
s/^DECR B\([01]\+\)0,/DECR B\1,UNPAD _1,/
t decr
s/^DECR B-\([01]\+\),/INCR B\1,NEG _,/
t incr
b loop-error


:icmp
# ICMP A B (A<B) (A==B) (A>B)
s/^ICMP B\(-\?[01]\+\) B\(-\?[01]\+\) (\([^),]*\)) (\([^),]*\)) (\([^),]*\)),/ISUB B\1 B\2,ICMP2 _ (\3) (\4) (\5),/
t isub
# Or with one arg
s/^ICMP B\(-\?[01]\+\) (\([^),]*\)) (\([^),]*\)) (\([^),]*\)),/ICMP2 B\1 (\2) (\3) (\4),/
t icmp2
b loop-error
:icmp2
s/^ICMP2 B-[01]\+ (\([^),]*\)) ([^),]*) ([^),]*),/RETURN \1,/
t return
s/^ICMP2 B0 ([^),]*) (\([^),]*\)) ([^),]*),/RETURN \1,/
t return
s/^ICMP2 B[01]\+ ([^),]*) ([^),]*) (\([^),]*\)),/RETURN \1,/
t return
b loop-error


:iadd    # [ADD algorithm]
# IADD is a hotspot
# negative numbers or zero
/^IADD B\([-0]\|[01]\+ B[-0]\)/{
    s/^IADD B0 B\(-\?[01]\+\),/RETURN B\1,/
    t return
    s/^IADD B\(-\?[01]\+\) B0,/RETURN B\1,/
    t return
    s/^IADD B-\([01]\+\) B\([01]\+\),/ISUB B\2 B\1,/
    t isub+
    s/^IADD B-\([01]\+\) B-\([01]\+\),/IADD B\1 B\2,NEG _,/
    t iadd+
    s/^IADD B\([01]\+\) B-\([01]\+\),/ISUB B\1 B\2,/
    t isub+
    b loop-error
}
:iadd+
# Reverse A and B and interleave the digits, padding with zero if one runs out
# Then add a carry digit to the front, and lookup the first three digits
# to get the rightmost digit of the result.
# IADD
#     A B => X=blank, IADD2 A B X
# IADD2
#     Aa Bb X => IADD2 A B Xab
#      x Bb X => IADD3 B X0b
#     Aa  x X => IADD4 A Xa0
#      x  x X => R=blank, IADD5 0X R
# IADD3
#     Bb X => IADD3 B X0b
#      x X => R=blank, IADD5 0X R
# IADD4
#     Aa X => IADD4 A Xa0
#      x X => R=blank, IADD5 0X R
# IADD5
#     cabX R => cr=(lookup cab), IADD5 cX rR
#     c    R => c=0, return R; c=1, return 1R
s/^IADD B\([01]\+\) B\([01]\+\),/IADD2 B\1 B\2 ,/
t iadd2
b loop-error
:iadd2
s/^IADD2 B\([01]*\)\([01]\) B\([01]*\)\([01]\) \([01]*\),/IADD2 B\1 B\3 \5\2\4,/
t iadd2
s/^IADD2 B B\([01]*\)\([01]\) \([01]*\),/IADD3 B\1 \30\2,/
t iadd3
s/^IADD2 B\([01]*\)\([01]\) B \([01]*\),/IADD4 B\1 \3\20,/
t iadd4
s/^IADD2 B B \([01]\+\),/IADD5 0\1 B,/
t iadd5
b loop-error
:iadd3
s/^IADD3 B\([01]*\)\([01]\) \([01]*\),/IADD3 B\1 \30\2,/
t iadd3
s/^IADD3 B \([01]\+\),/IADD5 0\1 B,/
t iadd5
b loop-error
:iadd4
s/^IADD4 B\([01]*\)\([01]\) \([01]*\),/IADD4 B\1 \3\20,/
t iadd4
s/^IADD4 B \([01]\+\),/IADD5 0\1 B,/
t iadd5
b loop-error
:iadd5
s/^IADD5 000\([01]*\) B\([01]*\),/IADD5 0\1 B0\2,/
t iadd5
s/^IADD5 001\([01]*\) B\([01]*\),/IADD5 0\1 B1\2,/
t iadd5
s/^IADD5 010\([01]*\) B\([01]*\),/IADD5 0\1 B1\2,/
t iadd5
s/^IADD5 011\([01]*\) B\([01]*\),/IADD5 1\1 B0\2,/
t iadd5
s/^IADD5 100\([01]*\) B\([01]*\),/IADD5 0\1 B1\2,/
t iadd5
s/^IADD5 101\([01]*\) B\([01]*\),/IADD5 1\1 B0\2,/
t iadd5
s/^IADD5 110\([01]*\) B\([01]*\),/IADD5 1\1 B0\2,/
t iadd5
s/^IADD5 111\([01]*\) B\([01]*\),/IADD5 1\1 B1\2,/
t iadd5
s/^IADD5 0 B\([01]\+\),/RETURN B\1,/
t return
s/^IADD5 1 B\([01]\+\),/RETURN B1\1,/
t return
b loop-error


:isub
# find negative number or zero
/^ISUB B\([-0]\|[01]\+ B[-0]\)/{
    s/^ISUB B0 B\(-\?[01]\+\),/NEG B\1,/
    t neg
    s/^ISUB B\(-\?[01]\+\) B0,/RETURN B\1,/
    t return
    s/^ISUB B-\([01]\+\) B\([01]\+\),/IADD B\1 B\2,NEG _,/
    t iadd+
    s/^ISUB B-\([01]\+\) B-\([01]\+\),/ISUB B\2 B\1,/
    t isub+
    s/^ISUB B\([01]\+\) B-\([01]\+\),/IADD B\1 B\2,/
    t iadd+
    b loop-error
}
# Reverse and interleave the digits of A and B, padding with zero if necessary
# Add a borrow bit to the front and lookup the first three bits
# and add to the result.
# A might turn out less than B, so we call ISUB B A,NEG _.
# Rather than keep up with the original A and B the whole time,
# push TMP A B on the stack and either use it or remove it.
# ISUB
#     A B => X=blank, ISUB2 A B X,TMP A B
# ISUB2
#     Aa Bb X => ISUB2 A B Xab
#     x  Bb X => (A<B), call ISUB B A,NEG _
#     Aa x  X => ISUB3 A Xa0
#     x  x  X => R=blank, ISUB4 0X R
# ISUB3
#     Aa X => ISUB3 A Xa0
#     x  X => R=blank, ISUB4 0X R
# ISUB4 (c = borrow)
#     cabX R => rc=(lookup cab), ISUB4 cX rR
#     0x   R => remove ISUB B A stack frame, return R
#     1x   R => (A<B), call ISUB B A,NEG _
:isub+
s/^ISUB B\([01]\+\) B\([01]\+\),/ISUB2 B\1 B\2 ,TMP B\1 B\2,/
t isub2
b loop-error
:isub2
s/^ISUB2 B\([01]*\)\([01]\) B\([01]*\)\([01]\) \([01]*\),/ISUB2 B\1 B\3 \5\2\4,/
t isub2
s/^ISUB2 B B[01]\+ [01]*,TMP B\([01]\+\) B\([01]\+\),/ISUB B\2 B\1,NEG _,/
t isub+
s/^ISUB2 B\([01]*\)\([01]\) B \([01]*\),/ISUB3 B\1 \3\20,/
t isub3
s/^ISUB2 B B \([01]\+\),/ISUB4 0\1 B,/
t isub4
b loop-error
:isub3
s/^ISUB3 B\([01]*\)\([01]\) \([01]*\),/ISUB3 B\1 \3\20,/
t isub3
s/^ISUB3 B \([01]\+\),/ISUB4 0\1 B,/
t isub4
b loop-error
:isub4
s/^ISUB4 000\([01]*\) B\([01]*\),/ISUB4 0\1 B0\2,/
t isub4
s/^ISUB4 001\([01]*\) B\([01]*\),/ISUB4 1\1 B1\2,/
t isub4
s/^ISUB4 010\([01]*\) B\([01]*\),/ISUB4 0\1 B1\2,/
t isub4
s/^ISUB4 011\([01]*\) B\([01]*\),/ISUB4 0\1 B0\2,/
t isub4
s/^ISUB4 100\([01]*\) B\([01]*\),/ISUB4 1\1 B1\2,/
t isub4
s/^ISUB4 101\([01]*\) B\([01]*\),/ISUB4 1\1 B0\2,/
t isub4
s/^ISUB4 110\([01]*\) B\([01]*\),/ISUB4 0\1 B0\2,/
t isub4
s/^ISUB4 111\([01]*\) B\([01]*\),/ISUB4 1\1 B1\2,/
t isub4
s/^ISUB4 0 B0*\([01]\+\),TMP B[01]\+ B[01]\+,/RETURN B\1,/
t return
s/^ISUB4 1 B[01]\+,TMP B\([01]\+\) B\([01]\+\),/ISUB B\2 B\1,NEG _,/
t isub+


:imult    # multiply; [MULT algorithm]
# IMULT is a hotspot
# get rid of signs
/^IMULT B\(-\|[01]\+ B-\)/{
	s/^IMULT B-\([01]\+\) B-\([01]\+\),/IMULT B\1 B\2,/
	t imult+
	s/^IMULT B-\?\([01]\+\) B-\?\([01]\+\),/IMULT B\1 B\2,NEG _,/
	t imult+
	b loop-error
}
# IMULT
#     A  0 => 0
#     A1 B => IMULT A B0,IADD _ B
#     A0 B => IMULT A B0
#     0  B => return 0
#     1  B => B
:imult+
# shortcut IMULT Bn B0; IMULT B0 Bn will be short anyway.
s/^IMULT B[01]\+ B0,/RETURN B0,/
t return
:imult2
# A0=1, A'=shr A; B'=shl B; (IMULT A' B')+B
s/^IMULT B\([01]\+\)1 B\([01]\+\),/IMULT B\1 B\20,IADD _ B\2,/
t imult2
# A0=0, A'=shr A; B'=shl B; IMULT A' B'
s/^IMULT B\([01]\+\)0 B\([01]\+\),/IMULT B\1 B\20,/
t imult2
# A=0, RETURN B0
s/^IMULT B0 B[01]\+,/RETURN B0,/
t return
# A=1, RETURN B
s/^IMULT B1 B\([01]\+\),/RETURN B\1,/
t return
b loop-error


:divrem
# remove signs
/^DIVREM B\(-\|[01]\+ B-\)/{
    s/^DIVREM B-\([01]\+\) B-\([01]\+\),/DIVREM B\1 B\2,/
    t divrem+
    s/^DIVREM B-\?\([01]\+\) B-\?\([01]\+\),/DIVREM B\1 B\2,LCHANGE B0 _,NEG _,/
    t divrem+
    b loop-error
}
:divrem+
# Check for divison by zero
/^DIVREM B[01]\+ B0,/{
    p
    s/^/ERROR DIVREM Divison by zero,/
    b error
}
s/^DIVREM B\([01]\+\) B\([01]\+\),/DIVREM2 B\1 B\2 B B B,/
t divrem2
b loop-error
:divrem2
# DIVREM2 D F C E R =>
#     if D=blank, and F=blank, then len(A)==len(B),
#                 if E<=C, R=R1., r=(C-E)
#                 else,    R=R0., r=C
#                 return R,r
#     if D=blank, then there are no more A digits, but some of B remains,
#                 R=R0, r=C, return R,r (r=EF?)
#     if F=blank, then all of B is copied (E==B) but some of A remains,
#                 if E<=C, R=R1, r=(C-E), C=rd, D=rightbits(D)
#                 else,    R=R0,  C=Cd, D=rightbits(D)
#     else, C=Cd, D=rightbits(D), E=Ef, F=rightbits(F)
# DIVREM2 D F C E R
# if D=blank and F=blank
s/^DIVREM2 B B B\([01]\+\) B\([01]\+\) B\([01]*\),/ISUB B\1 B\2,SUBST REMAINDER _,ICMP {REMAINDER} (DIVREM3 B\30 B\1) (DIVREM3 B\31 B0) (DIVREM3 B\31 {REMAINDER}),!,_,/
t isub
# if D=blank
s/^DIVREM2 B B[01]\+ B\([01]\+\) B[01]\+ B\([01]*\),/DIVREM3 B\20 B\1,/
t divrem3
# if F=blank
s/^DIVREM2 B\([01]\)\([01]*\) B B\([01]\+\) B\([01]\+\) B\([01]*\),/ISUB B\3 B\4,SUBST REMAINDER _,UNPAD B\3\1,SUBST C1 _,UNPAD B\50,SUBST R0 _,UNPAD B\51,SUBST R1 _,ICMP {REMAINDER} (DIVREM2 B\2 B {C1} B\4 {R0}) (DIVREM2 B\2 B B\1 B\4 {R1}) (DIVREM2 B\2 B {REMAINDER}\1 B\4 {R1}),!,!,!,!,_,/
t isub+
# else
s/^DIVREM2 B\([01]\)\([01]*\) B\([01]\)\([01]*\) B\([01]*\) B\([01]*\) B\([01]*\),/UNPAD B\5\1,SUBST C _,UNPAD B\6\3,SUBST E _,DIVREM2 B\2 B\4 {C} {E} B\7,!,!,/
t unpad
b loop-error
:divrem3
s/^DIVREM3 B0*\([01]\+\) B0*\([01]\+\),/RETURN <B\1:B\2:>,/
t return


:idiv
# divide and get quotient
s/^IDIV B\(-\?[01]\+\) B\(-\?[01]\+\),/DIVREM B\1 B\2,LREF B0 _,/
t divrem
b loop-error


:rem    # divide and get remainder
# This isn't the mathematical definition.
s/^REM B\(-\?[01]\+\) B\(-\?[01]\+\),/DIVREM B\1 B\2,LREF B1 _,/
t divrem


########################################################################
## FRACTIONAL REALIZE[12], FTRIM, FADD, FSUB, FMULT, FDIV, SQRT, FCMP ##
########################################################################


:realize1
# REALIZE1 Babbbbbbbb => Ba.bbbbbbbb
# REALIZE1 Bbbbbbbbb => B0.bbbbbbbb
# REALIZE1 B0 => B0.00000000
# REALIZE1 Bbbbb => B0.0000bbbb
# Add a decimal to a bunch of bits (make real = realize)
# fractional portion is always eight bits
s/^REALIZE1 B\(-\?\)0*\([01]\+\)\([01]\{8\}\),/RETURN B\1\2.\3,/
t return
# fewer than eight bits given
s/^REALIZE1 B\(-\?\)\([01]\+\),/REALIZE1 B\100000000\2,/
t realize1
b loop-error

:realize2
# REALIZE2 Babbbbbbbbcccccccc => Ba.bbbbbbbb
# REALIZE2 Bbbbbbbbbcccccccc => B0.bbbbbbbb
# REALIZE2 B0 => B0.00000000
# REALIZE2 Bbbbb => B0.00000000
# for example, A.B*C.D=>AB*CD=EFG=>E.F, drop the last 8 bits
s/^REALIZE2 B\(-\?\)0*\([01]\+\)\([01]\{8\}\)[01]\{8\},/RETURN B\1\2.\3,/
t return
# fewer than 16 bits given
s/^REALIZE2 B\(-\?\)\([01]\+\),/REALIZE2 B\10000000000000000\2,/
t realize2
b loop-error


:ftrim
# Take a (possibly left zero padded) fraction that might have too many
# or too few fractional parts.  Fix it.
s/^FTRIM B0*\([01]\+\.[01]\{8\}\)[01]*,/RETURN B\1,/
t return
s/^FTRIM B\([01]\+\.[01]*\),/FTRIM B\100000000,/
t ftrim


:fapprox
# FAPPROX A.B C.D
# Is A.B == C.D allowing the last decimal to vary
s/^FAPPROX B\(-\?[01]\+\)\.\([01]\{7\}\)[01] B\1\.\2[01],/RETURN STRUE,/
t return
s/^FAPPROX B0\.00000000 B-0\.00000001,/RETURN STRUE,/
t return
s/^FAPPROX B-0\.00000001 B0\.00000000,/RETURN STRUE,/
t return
s/^FAPPROX B-\?[01]\+\.[01]\+ B-\?[01]\+\.[01]\+,/RETURN SFALSE,/
t return
b loop-error


:fadd
# A.B + C.D => AB + CD = EF => E.F
s/^FADD B\(-\?[01]\+\)\.\([01]\{8\}\) B\(-\?[01]\+\)\.\([01]\{8\}\),/UNPAD B\1\2,SUBST T1 _,UNPAD B\3\4,IADD {T1} _,!,REALIZE1 _,/
t unpad
b loop-error


:fsub
# A.B + C.D => AB - CD = EF => E.F
s/^FSUB B\(-\?[01]\+\)\.\([01]\{8\}\) B\(-\?[01]\+\)\.\([01]\{8\}\),/UNPAD B\1\2,SUBST T1 _,UNPAD B\3\4,ISUB {T1} _,!,REALIZE1 _,/
t unpad
b loop-error


:fmult
# A.B * C.D => AB * CD = EFG => E.F
s/^FMULT B\(-\?[01]\+\)\.\([01]\{8\}\) B\(-\?[01]\+\)\.\([01]\{8\}\),/UNPAD B\1\2,SUBST T1 _,UNPAD B\3\4,IMULT {T1} _,!,REALIZE2 _,/
t unpad
b loop-error


:fdiv
# remove signs
/^FDIV B\(-\|[01]\+\.[01]\+ B-\)/{
    s/^FDIV B-\([01]\+\.[01]\+\) B-\([01]\+\.[01]\+\),/FDIV B\1 B\2,/
    t fdiv+
    s/^FDIV B-\?\([01]\+\.[01]\+\) B-\?\([01]\+\.[01]\+\),/FDIV B\1 B\2,NEG _,/
    t fdiv+
    b loop-error
}
:fdiv+
# Check for divison by zero
/^FDIV B[01]\+\.[01]\+ B0.00000000,/{
    p
    s/^/ERROR Divison by zero in FDIV,/
    b error
}
# FDIV W.X Y.Z => A=WX, B=YZ, D=A, F=B, C=E=blank, R=B0, FDIV2 D F C E R
s/^FDIV B\([01]\+\)\.\([01]\{8\}\) B\([01]\+\)\.\([01]\{8\}\),/UNPAD B\1\2,SUBST D _,UNPAD B\3\4,FDIV2 {D} _ B B B,!,/
t unpad
b loop-error

# The basic algorithm is like long division.
# assume b is the leftmost digit of B
# R=blank, E=blank
# A/B => C=a,
#     if B<C, R=R1, r=(C-B), C=ra, A=rightbits(A)
#     else: R=R0, C=Ca, A=rightbits(A)
#     if A=blank, R=R., a=0
# If len(B) < len(A), we have unnecessary subtractions.
# We drop these at the beginning by copying digits of A and B
# until one is empty.
# FDIV A B = A/B => C=blank, D=A, E=blank, F=B, R=blank, FDIV2 C D E F R
# FDIV2 D F C E R =>
#     if D=blank, and F=blank, then len(A)==len(B),
#                 if E<=C, R=R1., r=(C-E), C=r0
#                 else,    R=R0., C=C0
#                 FDIV4 C E R
#     if D=blank, then there are no more A digits, but some of B remains,
#                 E=Ed, F=rightbits(F), R=R0., C=C0
#                 FDIV3 F C E R
#     if F=blank, then all of B is copied (E==B) but some of A remains,
#                 if E<=C, R=R1, r=(C-E), C=rd, D=rightbits(D)
#                 else,    R=R0, C=Cd, D=rightbits(D)
#     else, C=Cd, D=rightbits(D), E=Ef, F=rightbits(F)
# FDIV3 F C E R =>
#     if F=blank, then len(A)==len(B), FDIV4 C E R
#     else,       E=Ef, F=rightbits(F), C=C0, R=R0, FDIV3 F C E R
# FDIV4 C B R = C/B =>
#     if R=n.nnnnnnnn, return (ftrim R)
#     if B<=C, R=R1, r=(C-B), C=r0
#     else,    R=R0, C=C0

:fdiv2
# FDIV2 D F C E R
# if D=blank and F=blank
s/^FDIV2 B B B\([01]\+\) B\([01]\+\) B\([01]*\),/ISUB B\1 B\2,SUBST REMAINDER _,ICMP {REMAINDER} (FDIV4 B\10 B\2 B\30.) (FDIV4 B0 B\2 B\31.) (FDIV4 {REMAINDER}0 B\2 B\31.),!,_,/
t isub+
# if D=blank
s/^FDIV2 B B\([01]\)\([01]*\) B\([01]\+\) B\([01]\+\) B\([01]*\),/FDIV3 B\2 B\30 B\4\1 B\50.,/
t fdiv3
# if F=blank
s/^FDIV2 B\([01]\)\([01]*\) B B\([01]\+\) B\([01]\+\) B\([01]*\),/ISUB B\3 B\4,SUBST REMAINDER _,ICMP {REMAINDER} (FDIV2 B\2 B B\3\1 B\4 B\50) (FDIV2 B\2 B B\1 B\4 B\51) (FDIV2 B\2 B {REMAINDER}\1 B\4 B\51),!,_,/
t isub+
# else
s/^FDIV2 B\([01]\)\([01]*\) B\([01]\)\([01]*\) B\([01]*\) B\([01]*\) B\([01]*\),/FDIV2 B\2 B\4 B\5\1 B\6\3 B\7,/
t fdiv2
b loop-error
:fdiv3
# FDIV3 F C E R
s/^FDIV3 B B\([01]\+\) B\([01]\+\) B\([01.]\+\),/FDIV4 B\1 B\2 B\3,/
t fdiv4
s/^FDIV3 B\([01]\)\([01]*\) B\([01]\+\) B\([01]\+\) B\([01.]\+\),/UNPAD B\30,FDIV3 B\2 _ B\4\1 B\50,/
t unpad
b loop-error
:fdiv4
# FDIV4 C B R = C/B =>
#     if R=n.nnnnnnnn, return (ftrim R)
#     if B<=C, R=R1, r=(C-B), C=r0
#     else,    R=R0, C=C0
s/^FDIV4 B[01]\+ B[01]\+ B\([01]\+\.[01]\{8\}\)[01]*,/FTRIM B\1,/
t ftrim
s/^FDIV4 B0*\([01]\+\) B\([01]\+\) B\([01.]\+\),/ISUB B\1 B\2,SUBST REMAINDER _,ICMP {REMAINDER} (FDIV4 B\10 B\2 B\30) (FDIV4 B0 B\2 B\31) (FDIV4 {REMAINDER}0 B\2 B\31),!,_,/
t isub+
b loop-error


:sqrt    # [SQRT algorithm]
# Babylonian Method (http://en.wikipedia.org/wiki/Square_root)
# sqrt(r) = ?
# x = guess
# repeat x = avg(x,r/x), until x = x2
s/^SQRT B\([01]\+\.[01]\+\),/SQRT2 B\1 B0.00000000 B101.00000000,/
t sqrt2
b loop-error
:sqrt2
# SQRT2 r xprev x; for little numbers (r<100)
# x' <- (r/x + x) * 1/2; if abs(xprev-x) < eps, return x, else xprev<-x; x<-x'
s/^SQRT2 B\([01]\+\.[01]\+\) B[01]\+\.[01]\+ B\([01]\+\.[01]\+\),/FDIV B\1 B\2,FADD _ B\2,FMULT _ B0.10000000,SUBST Xp _,FAPPROX {Xp} B\2,SELECT _ (STRUE RETURN {Xp}) (SFALSE SQRT2 B\1 B\2 {Xp}) (FAIL),!,_,/
t fdiv
b loop-error


:fcmp
# FCMP A B (R1) (R2) (R3)
# t <- A-B; if t<0, return R1; if t==0, return R2; if t>0, return R3
s/^FCMP B\(-\?[01]\+\.[01]\+\) B\(-\?[01]\+\.[01]\+\) \(([^),]*) ([^),]*) ([^),]*)\),/FSUB B\1 B\2,SUBST T _,FAPPROX {T} B0.00000000,FCMP2 _ {T} \3,!,/
t fsub
b loop-error
:fcmp2
# =0?
s/^FCMP2 STRUE B-\?[01+\.[01]\+ ([^),]*) (\([^),]\+\)) ([^),]*),/RETURN \1,/
t return
s/^FCMP2 SFALSE B-[01]\+\.[01]\+ (\([^),]\+\)) ([^),]*) ([^),]*),/RETURN \1,/
t return
s/^FCMP2 SFALSE B[01]\+\.[01]\+ ([^),]*) ([^),]*) (\([^),]\+\)),/RETURN \1,/
t return
b loop-error


#################################################
## VARIABLES: REF, SET, DEFAULT, DELETE, FLUID ##
#################################################


:ref
# reference a global variable
s/^REF \([^,]\+\),\([^@]*\)@\(.*,VAR \1 \([^,]\+\).*\)/RETURN \4,\2@\3/
t return
b loop-error


:set
# set a global variable
# if preexists, replace old value.
s/^SET \([^ ]\+\) \([^,]\+\),\([^@]*\)@\(.*\),VAR \1 [^,]\+\(.*\)/\3@\4,VAR \1 \2\5/
t cycle
# if doesn't preexist, add it
s/^SET \([^ ]\+\) \([^,]\+\),\([^@]*\)@/\3@,VAR \1 \2/
t cycle
b loop-error


:default
# set only if undefined.
# defined? ignore
s/^DEFAULT \([^ ]\+\) [^,]\+,\([^@]*@.*,VAR \1 [^,]\+\)/\2/
t cycle
s/^DEFAULT \([^ ]\+\) \([^,]\+\),/SET \1 \2,/
t set
b loop-error


:delete
# delete a global variable
s/^DELETE \([^,]\+\),\([^@]*\)@\(.*\),VAR \1 [^,]\+/\2@\3/
t cycle
b loop-error


:fluid
# Fluid frames hold stack allocated variables.
# If it's on top, just pop it.
s/^FLUID [^ ]\+ [^,]\+,//
t cycle
b loop-error


#########################################
## LIST UTILITIES: LREF, CONS, LCHANGE ##
#########################################


:lref
# a list is something like these:  <> <B1:> <B0:B1:B10:B11:>
# We can't have *syntactically* recursive structures,
# but we can have lists like this:  <V6:> where V6 is <B0:B1:>
# They have to be handled specially, though.
s/^LREF B[01]\+ <>,/ERROR LREF: attempted to index into empty list,/
t error
# one digit
/^LREF B[01] /{
    s/^LREF B0 <\([^:]\+\):[^>]*>,/RETURN \1,/
    t return
    s/^LREF B1 <[^:]\+:\([^:]\+\):[^>]*>,/RETURN \1,/
    t return
    b loop-error
}
# It's not the first two positions, so we can subtract two
# from the index and remove two objects from the list
s/^LREF B\([01]\+\) <[^:]\+:[^:]\+:\([^>]*\)>,/DECR B\1,DECR _,LREF _ <\2>,/
t decr
b loop-error


:cons
s/^CONS \([^ ]\+\) <\([^>]*\)>,/RETURN <\1:\2>,/
t return
b loop-error


:lchange
# LCHANGE <index-n> <a0,a1,a2,...,an,...>,f _,
# => f an,RETURN <a0,a1,a2,...,_,...>,
# This seems kinda dumb.
s/^LCHANGE B0 <\([^:]\+\):\([^>]*\)>,\([^,]\+\),/RETURN \1,\3,RETURN <_:\2>,/
t return
s/^LCHANGE B1 <\([^:]\+\):\([^:]\+\):\([^>]*\)>,\([^,]\+\),/RETURN \2,\4,RETURN <\1:_:\3>,/
t return
/^LCHANGE B[01]\+ <[^>]*>,/{
    p
    s/^/ERROR LCHANGE Unimplemented with indices other than B0 and B1.,/
    t error
}
b loop-error


######################################################################
## VECTOR/TUPLE FUNCTIONS: VADD, VSUB, SCALARMULT, CROSS, UNIT, MAG ##
######################################################################


:vadd
# VECTOR ADDITION: (x1 y1 z1) + (x2 y2 z2) = (x1+x2 y1+y2 z1+z2)
# integers
s/^VADD V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)# V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)#,/IADD B\1 B\4,SUBST X _,IADD B\2 B\5,SUBST Y _,IADD B\3 B\6,RETURN V#{X}#{Y}#_#,!,!,/
t add
# both fractions
s/^VADD V#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)# V#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#,/FADD B\1 B\4,SUBST X _,FADD B\2 B\5,SUBST Y _,FADD B\3 B\6,RETURN V#{X}#{Y}#_#,!,!,/
t fadd
# first integer vector, second fractional
s/^VADD V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)# V#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#,/FADD B\1.00000000 B\4,SUBST X _,FADD B\2.00000000 B\5,SUBST Y _,FADD B\3.00000000 B\6,RETURN V#{X}#{Y}#_#,!,!,/
t fadd
# first fractional, second integers
s/^VADD V#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)# B#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)#,/FADD B\1 B\4.00000000,SUBST X _,FADD B\2 B\5.00000000,SUBST Y _,FADD B\3 B\6.00000000,RETURN B#{X}#{Y}#_#,!,!,/
t fadd
b loop-error


:vsub
# VECTOR SUBTRACTION: (x1 y1 z1)-(x2 y2 z2) = (x1-x2 y1-y2 z1-z2)
s/^VSUB V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)# V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)#,/SUB B\1 B\4,SUBST X _,SUB B\2 B\5,SUBST Y _,SUB B\3 B\6,RETURN V#{X}#{Y}#_#,!,!,/
t sub
b loop-error


:scalarmult
# SCALAR MULTIPLICATION: s * (x y z) = (s*x s*y s*z)
# integer scalar
#   tuple of integers
s/^SCALARMULT B\(-\?[01]\+\) V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)#,/IMULT B\1 B\2,SUBST X _,IMULT B\1 B\3,SUBST Y _,IMULT B\1 B\4,RETURN V#{X}#{Y}#_#,!,!,/
t imult
#   tuple of fractions
s/^SCALARMULT B\(-\?[01]\+\) V#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#,/FMULT B\1.00000000 B\2,SUBST X _,FMULT B\1.00000000 B\3,SUBST Y _,FMULT B\1.00000000 B\4,RETURN V#{X}#{Y}#_#,!,!,/
t fmult
# fractional scalar
#   tuple of integers
s/^SCALARMULT B\(-\?[01]\+\.[01]\+\) V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)#,/FMULT B\1 B\2.00000000,SUBST X _,FMULT B\1 B\3.00000000,SUBST Y _,FMULT B\1 B\4.00000000,RETURN V#{X}#{Y}#_#,!,!,/
t fmult
#   tuple of fractions
s/^SCALARMULT B\(-\?[01]\+\.[01]\+\) V#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#,/FMULT B\1 B\2,SUBST X _,FMULT B\1 B\3,SUBST Y _,FMULT B\1 B\4,RETURN V#{X}#{Y}#_#,!,!,/
t fmult
b loop-error


:cross    # [CROSS PRODUCT algorithm]
# CROSS P1 P2: cross product; ok Ps are just tuples
# (x1 y1 z1)X(x2 y2 z2) = (y1z2-z1y2 z1x2-x1z2 x1y2-y1x2)
# integers
s/^CROSS V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)# V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)#,/IMULT B\2 B\6,SUBST T1 _,IMULT B\3 B\5,ISUB {T1} _,!,SUBST X _,IMULT B\3 B\4,SUBST T2 _,IMULT B\1 B\6,ISUB {T2} _,SUBST Y _,IMULT B\1 B\5,SUBST T3 _,IMULT B\2 B\4,ISUB {T3} _,RETURN V#{X}#{Y}#_#,!,!,!,!,/
t imult
# fractions
s/^CROSS V#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)# V#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#B\(-\?[01]\+\.[01]\+\)#,/FMULT B\2 B\6,SUBST T1 _,FMULT B\3 B\5,SUB {T1} _,!,SUBST X _,FMULT B\3 B\4,SUBST T2 _,FMULT B\1 B\6,FSUB {T2} _,SUBST Y _,FMULT B\1 B\5,SUBST T3 _,FMULT B\2 B\4,FSUB {T3} _,RETURN V#{X}#{Y}#_#,!,!,!,!,/
t fmult
b loop-error


:unit
# UNIT P
# u = (Px/m Py/m Pz/m), m=mag(P)
s/^UNIT V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)#,/MAG V#B\1#B\2#B\3#,SUBST M _,FDIV B\1.00000000 {M},SUBST X _,FDIV B\2.00000000 {M},SUBST Y _,FDIV B\3.00000000 {M},!,RETURN V#{X}#{Y}#_#,!,!,/
t mag
b loop-error


:mag
# MAG P = sqrt(Px*Px + Py*Py + Pz*Pz)
s/^MAG V#B\(-\?[01]\+\)#B\(-\?[01]\+\)#B\(-\?[01]\+\)#,/IMULT B\1 B\1,SUBST T1 _,IMULT B\2 B\2,IADD {T1} _,!,SUBST T2 _,IMULT B\3 B\3,IADD {T2} _,!,SQRT _.00000000,/
t imult
b loop-error


#    :gcd    # [GCD algorithm]
#    # gcd A 0 => A
#    # gcd A B => gcd B (A mod B)
#    # This could really use a faster MOD
#    s/^GCD B\([01]\+\) B0,/RETURN B\1,/
#    t return
#    s/^GCD B\([01]\+\) B\([01]\+\),/REM B\1 B\2,GCD B\2 _,/
#    t rem
#    b loop-error


################
## DTOB, BTOD ##
################


:dtob
s/^DTOB \(-\?[0-9]\+\),/IDTOB \1,/
t idtob
s/^DTOB \(-\?[0-9]\+\)\.\([0-9]\+\),/IDTOB \1,SUBST A _,FDTOB \2,RETURN {A}._,!,/
t idtob
b loop-error


:idtob    # convert decimal to binary
s/^IDTOB -\([0-9]\+\),/IDTOB \1,NEG _,/
:idtob+
s/^IDTOB \([0-9]\+\),/IDTOB2 \1 B1,/
t idtob2
b loop-error
:idtob2
s/^IDTOB2 \([0-9]\) B1,/IDTOB1 \1,/
t idtob1
s/^IDTOB2 \([0-9]\) B\([01]\+\),/IDTOB1 \1,IMULT _ B\2,/
t idtob1
s/^IDTOB2 \([0-9]\+\)\([0-9]\) B\([01]\+\),/IDTOB1 \2,IMULT _ B\3,SUBST T1 _,IMULT B\3 B1010,IDTOB2 \1 _,ADD {T1} _,!,/
t idtob1
b loop-error
:idtob1
s/^IDTOB1 \([0-9]\),/SELECT \1 (0 0) (1 1) (2 10) (3 11) (4 100) (5 101) (6 110) (7 111) (8 1000) (9 1001) (FAIL),RETURN B_,/
t select
b loop-error

:fdtob
# 0.D => FDTOB D => A=0.0, M=B1010, FDTOB2 D A M
# FDTOB2 dD A M
#     if dD=blank, return A
#     q=(lookup d), A=A+(q/M), M=M*B1010, FDTOB2 D A M
# This has some round off error (0.25=B0.00111111)
# Maybe A can be seeded differently?  Compute to 9 digits and round?
s/^FDTOB \([0-9]\+\),/FDTOB2 \1 B0.00000000 B1010,/
t fdtob2
b loop-error
:fdtob2
s/^FDTOB2  B0.\([01]\+\) B[01]\+,/RETURN \1,/
t return
s/^FDTOB2 \([0-9]\)\([0-9]*\) B\(0.[01]\+\) B\([01]\+\),/SELECT \1 (0 0) (1 1) (2 10) (3 11) (4 100) (5 101) (6 110) (7 111) (8 1000) (9 1001) (FAIL),FDIV B_.00000000 B\4.00000000,FADD B\3 _,SUBST A _,MULT B\4 B1010,FDTOB2 \2 {A} _,!,/
t select
b loop-error


:pow10
# POW10 N => 10**n
s/^POW10 B0,/RETURN B1,/
s/^POW10 B1,/RETURN B1010,/
s/^POW10 B10,/RETURN B1100100,/
s/^POW10 B11,/RETURN B1111101000,/
s/^POW10 B100,/RETURN B10011100010000,/
s/^POW10 B101,/RETURN B11000011010100000,/
s/^POW10 B110,/RETURN B11110100001001000000,/
s/^POW10 B111,/RETURN B100110001001011010000000,/
t return
s/^POW10 B1000,/RETURN B101111101011110000100000000,/
s/^POW10 B1001,/RETURN B111011100110101100101000000000,/
s/^POW10 B1010,/RETURN B1001010100000010111110010000000000,/
s/^POW10 B1011,/RETURN B1011101001000011101101110100000000000,/
s/^POW10 B1100,/RETURN B1110100011010100101001010001000000000000,/
s/^POW10 B1101,/RETURN B10010001100001001110011100101010000000000000,/
s/^POW10 B1110,/RETURN B10110101111001100010000011110100100000000000000,/
s/^POW10 B1111,/RETURN B11100011010111111010100100110001101000000000000000,/
t return
s/^POW10 B10000,/RETURN B100011100001101111001001101111110000010000000000000000,/
s/^POW10 B10001,/RETURN B101100011010001010111100001011101100010100000000000000000,/
s/^POW10 B10010,/RETURN B110111100000101101101011001110100111011001000000000000000000,/
s/^POW10 B10011,/RETURN B1000101011000111001000110000010010001001111010000000000000000000,/
s/^POW10 B10100,/RETURN B1010110101111000111010111100010110101100011000100000000000000000000,/
s/^POW10 B10101,/RETURN B1101100011010111001001101011011100010111011110101000000000000000000000,/
s/^POW10 B10110,/RETURN B10000111100001100111100000110010011011101010110010010000000000000000000000,/
s/^POW10 B10111,/RETURN B10101001011010000001011000111111000010100101011110110100000000000000000000000,/
s/^POW10 B11000,/RETURN B11010011110000100001101111001110110011001110110110100001000000000000000000000000,/
s/^POW10 B11001,/RETURN B100001000101100101010001011000010100000000010100100001001010000000000000000000000000,/
s/^POW10 B11010,/RETURN B101001010110111110100101101110011001000000011001101001011100100000000000000000000000000,/
s/^POW10 B11011,/RETURN B110011101100101110001111001001111111010000100000000011110011101000000000000000000000000000,/
s/^POW10 B11100,/RETURN B1000000100111111001110010111100011111000100101000000100110000100010000000000000000000000000000,/
s/^POW10 B11101,/RETURN B1010000110001111000001111101011100110110101110010000101111100101010100000000000000000000000000000,/
s/^POW10 B11110,/RETURN B1100100111110010110010011100110100000100011001110100111011011110101001000000000000000000000000000000,/
s/^POW10 B11111,/RETURN B1111110001101111011111000100000001000101100000010010001010010110010011010000000000000000000000000000000,/
t return
s/^POW10 B\([01]\+\),/DECR B\1,POW10 _,IMULT B1010 _,/
t decr


:btod
s/^BTOD B\(-\?[01]\+\),/IBTOD B\1,/
t ibtod
s/^BTOD B\(-\?[01]\+\)\.\([01]\+\),/IBTOD B\1,FBTOD _. B0.\2,/
t ibtod
b loop-error


:ibtod
# integer binary to decimal (base 10)
# IBTOD
#     0 => return 0
#    -A => -(IBTOD A)
#     A => M=0, IBTOD2 A M
# IBTOD2 A M => m=POW10 M, C=(m<A), IBTOD3 C A M
# IBTOD3
#     LT A M => (m<A): M=M+1, IBTOD2 A M
#     EQ A M => (m=A): M=M-1, R=blank, IBTOD4 <1:0:> M R
#     GT A M => (m>A): M=M-1, IBTOD5 A M R
# IBTOD4 <q:r:> M R => Q=(lookup q), R=RQ, A=r, IBTOD5 A M R
# IBTOD5
#     A -1 R => return R
#     A M R => m=POW10 M, <q:r:>=A/m, M=M-1, IBTOD4 <q:r:> M R
s/^IBTOD B0,/RETURN 0,/
t return
s/^IBTOD B-\([01]\+\),/IBTOD2 B\1 B0,NEG _,/
s/^IBTOD B\([01]\+\),/IBTOD2 B\1 B0,/
t ibtod2
b loop-error
:ibtod2
s/^IBTOD2 B\([01]\+\) B\([01]\+\),/POW10 B\2,ICMP _ B\1 (LT) (EQ) (GT),IBTOD3 _ B\1 B\2,/
t pow10
b loop-error
:ibtod3
s/^IBTOD3 LT B\([01]\+\) B\([01]\+\),/INCR B\2,IBTOD2 B\1 _,/
t incr
s/^IBTOD3 EQ B[01]\+ B\([01]\+\),/DECR B\1,IBTOD4 <B1:B0:> _ ,/
t decr
s/^IBTOD3 GT B\([01]\+\) B\([01]\+\),/DECR B\2,IBTOD5 B\1 _ ,/
t decr
b loop-error
:ibtod4
s/^IBTOD4 <B\([01]\+\):B\([01]\+\):> B\(-\?[01]\+\) \([0-9]*\),/SELECT \1 (0 0) (1 1) (10 2) (11 3) (100 4) (101 5) (110 6) (111 7) (1000 8) (1001 9) (FAIL),IBTOD5 B\2 B\3 \4_,/
t select
b loop-error
:ibtod5
s/^IBTOD5 B[01]\+ B-1 \([0-9]\+\),/RETURN \1,/
t return
s/^IBTOD5 B\([01]\+\) B\([01]\+\) \([0-9]*\),/POW10 B\2,SUBST m _,DECR B\2,SUBST M _,DIVREM B\1 {m},IBTOD4 _ {M} \3,!,!,/
t pow10
b loop-error


:fbtod
# R.A => FBTOD R. A =>
#     if R=n.nnnnnnnn, return R
#     A=A*10=C.D, FBTOD2 R C.D
# FBTOD2 R C.D =>
#     if C.D=0.0, return R
#     R=RC, A=D, FBTOD R A
s/^FBTOD \(-\?[0-9]\+\.[0-9]\{8\}\) B[01.]\+,/RETURN \1,/
t return
s/^FBTOD \(-\?[0-9]\+\.[0-9]*\) B\([01]\+\.[01]\+\),/FMULT B1010.00000000 B\2,FBTOD2 \1 _,/
t fmult
b loop-error
:fbtod2
s/^FBTOD2 \(-\?[0-9]\+\.[0-9]\+\) B0.00000000,/RETURN \1,/
t return
s/^FBTOD2 \(-\?[0-9]\+\.[0-9]*\) B\([01]\+\)\.\([01]\+\),/SELECT \2 (0 0) (1 1) (10 2) (11 3) (100 4) (101 5) (110 6) (111 7) (1000 8) (1001 9) (FAIL),FBTOD \1_ B0.\3,/
t select
b loop-error


#################################
## PRINT, PRINTOBJ, PRINTSTACK ##
#################################


:print
#s/^PRINT[^,]*,//
#t cycle
# a simple, dumb, print
h
s/^PRINT \?\([^,]*\),.*/\1/
p
g
# The idea was to allow PRINT to act as an identity filter,
# with printing as a side effect, but I never use it that way.
# I want a label.  I'm removing it because I want RETURN to
# become a PRINT when the stack ends, and when PRINT then becomes
# a RETURN, havoc ensues.
#s/^PRINT \?\([^,]*\),/RETURN \1,/
s/^PRINT \?\([^,]*\),//
t cycle
b loop-error


:printobj
# NUMBER
# without label
s/^PRINTOBJ B\(-\?[01]\+\(\.[01]\+\)\?\),/BTOD B\1,PRINT _ = B\1,RETURN B\1,/
t print
#
# with label
s/^PRINTOBJ \([^ ]*\) B\(-\?[01]\+\(\.[01]\+\)\?\),/BTOD B\2,PRINT \1 = _ = B\2,RETURN B\2,/
t btod
# 3 TUPLE
# without label
s/^PRINTOBJ V#B\(-\?[01]\+\(\.[01]\+\)\?\)#B\(-\?[01]\+\(\.[01]\+\)\?\)#B\(-\?[01]\+\(\.[01]\+\)\?\)#,/BTOD B\1,SUBST X _,BTOD B\3,SUBST Y _,BTOD B\5,SUBST Z _,PRINT ({X} {Y} {Z}) = (B\1 B\3 B\5),!,!,!,RETURN V#B\1#B\3#B\5#,/
t btod
# with label
s/^PRINTOBJ \([^ ,]\+\) V#B\(-\?[01]\+\(\.[01]\+\)\?\)#B\(-\?[01]\+\(\.[01]\+\)\?\)#B\(-\?[01]\+\(\.[01]\+\)\?\)#,/BTOD B\2,SUBST X _,BTOD B\4,SUBST Y _,BTOD B\6,SUBST Z _,PRINT \1 = ({X} {Y} {Z}) = (B\2 B\4 B\6),!,!,!,RETURN V#B\2#B\4#B\6#,/
t btod
# 5 TUPLE
# without label
s/^PRINTOBJ V#B\(-\?[01]\+\(\.[01]\+\)\?\)#B\(-\?[01]\+\(\.[01]\+\)\?\)#B\(-\?[01]\+\(\.[01]\+\)\?\)#B\(-\?[01]\+\(\.[01]\+\)\?\)#\([0-9A-F]\+\)#,/BTOD B\1,SUBST A _,BTOD B\3,SUBST B _,BTOD B\5,SUBST C _,BTOD B\7,SUBST D _,PRINT ({A} {B} {C} {D} \9) = (B\1 B\3 B\5 B\7 \9),!,!,!,!,RETURN V#B\1#B\3#B\5#B\7#\9#,/
t btod
# with label
# this would require >9 back references, so play a little loosely
s/^PRINTOBJ \([^ ,]\+\) V#B\([^#]*\)#B\([^#]*\)#B\([^#]*\)#B\([^#]*\)#\([0-9A-F]\+\)#,/BTOD B\2,SUBST A _,BTOD B\3,SUBST B _,BTOD B\4,SUBST C _,BTOD B\5,SUBST D _,PRINT \1 = ({A} {B} {C} {D} \6) = (B\2 B\3 B\4 B\5 \6),!,!,!,!,RETURN V#B\2#B\3#B\4#B\5#\6#,/
t btod
b loop-error


:printstack
p
s/^PRINTSTACK,//
t cycle
b loop-error


##############################################
## READCMD, ERROR, LOOP-ERROR, ASSERT, EXIT ##
##############################################


:readcmd    # read a new instruction (command) from the input stream
s/^READCMD,//
h
d


:error    # display error message and quit
/^ERROR /{
    s/^ERROR \([^,]\+\),.*/ERROR -- \1/
    q
}
b loop-error


:loop-error    # coding error; display stack and quit
#p
i\
** LOOP ERROR **
q


:assert
# ASSERT <expected-result> <expression>
s/^ASSERT \([^ ]\+\) \([^,]\+\),/\2,ASSERT2 \1 _ \2,/
t cycle
:assert2
# ASSERT2 <expected-result> <actual-result> <expression>
# expected = actual; move on
s/^ASSERT2 \([^ ]\+\) \1 [^,]\+,//
t cycle
# expected != actual; print message and fail
s/^ASSERT2 \([^ ]\+\) \([^ ]\+\) \([^,]\+\),/ERROR Assertion failed: \3 = \2; not \1,/
t error
b loop-error


:exit
s/^EXIT \([^,]*\),/PRINT \1,EXIT,/
t print
s/.*//
q


## END ########################################################################