--- /dev/null
+## Author: Heiner Marxen
+## Origin: http://www.drb.insel.de/~heiner/BB/hmMMsimu.awk
+##
+## Emulate a TM by a macro machine, and print a list of configs.
+## This program demonstrates/implements some of the acceleration
+## techniques used by Marxen/Buntrock to attack the 5-state busy beaver,
+## especially the construction of macro machines.
+##
+##=============================================================================
+## $Id: basics.awk,v 1.1 2010/05/06 17:24:17 heiner Exp $
+## Author: Heiner Marxen
+## Module: basics
+##
+## Some basic functions of general value.
+##
+## Prefixes: str_
+##-----------------------------------------------------------------------------
+
+function str_rep(str,cnt)
+{
+ if( cnt <= 2 ) {
+ if( cnt < 1 ) return ""
+ if( cnt < 2 ) return str
+ return str str
+ }
+ return str_rep(str str, int(cnt/2)) ((cnt%2) ? str : "")
+ ## FFS: iterate instead of recursion
+}
+
+function str_00(len)
+{
+ if( len <= 0 ) return ""
+ if( len <= 99 ) { ## sprintf is not safe for large lengthes!
+ return sprintf("%0*d", len, 0)
+ }
+ return str_rep("0", len)
+}
+
+function str_bb(len)
+{
+ if( len <= 0 ) return ""
+ if( len <= 99 ) { ## sprintf is not safe for large lengthes!
+ return sprintf("%*s", len, "")
+ }
+ return str_rep(" ", len)
+}
+
+function str_bi(len,str ,d)
+{
+ ## Insert blanks to reach length at least "len"
+ ## Can replace an sprintf("%*s", len,str)
+ d = len - length(str)
+ if( d <= 0 ) return str
+ return str_bb(d) str
+}
+
+function str_ba(len,str ,d)
+{
+ ## Append blanks to reach length at least "len"
+ ## Can replace an sprintf("%-*s", len,str)
+ d = len - length(str)
+ if( d <= 0 ) return str
+ return str str_bb(d)
+}
+
+##=============================================================================
+## Module initialization
+
+function basics_init() # module constructor
+{
+ if( basics_inited_ ) return # did it, already
+ basics_inited_ = 1 # committed to do so
+ # note our identity
+ g_Id[++g_IdCnt] = "$Id: basics.awk,v 1.1 2010/05/06 17:24:17 heiner Exp $"
+}
+
+##-----------------------------------------------------------------------------
+## End of module "basics"
+##=============================================================================
+##=============================================================================
+## $Id: bignum.awk,v 1.34 2010/05/06 17:58:14 heiner Exp $
+## Author: Heiner Marxen
+## Module: bignum
+##
+## Compute with arbitrarily large signed/unsigned numbers,
+## coded as a string the same way as AWK does in its automatic conversion.
+##
+## Prefixes: bix, bn
+##-----------------------------------------------------------------------------
+## awk problems:
+## - "nawk" of Solaris "optimizes": x="0" ; y = -x ;
+## such that y == "-0". It appears to not happen for: y = 0-x
+## - "nawk" of Solaris "optimizes": 0 * -1 --> -0
+## We try to compensate that bug, but we are not really sure.
+## Better avoid Solaris "nawk", if possible. FFS.
+## - "gawk" (at least) 3.1.4 and 3.1.5 sometimes does not invalidate the
+## numeric interpretation of a variable which is appended to. gawk-fix
+## On 4-Jul-2006 Arnold D. Robbins has provided a fix, which I assume
+## will make it into the next gawk release, but I'll leave in my workaround
+## for others which happen to have the buggy version.
+## - "mawk" (using version mawk-1.3.4-20100419): had a bug in "array.c".
+## Fix provided by Thomas Dickey 2010-05-05.
+## - "mawk" limits the length of the result of "sprintf()". With option
+## "-W sprintf=nnnnn" we can increase the limit, but it is still there.
+##-----------------------------------------------------------------------------
+## We expect to operate on strings representing (signed) integers.
+## The acceptable forms are: "0", "123" and "-567". We do not accept
+## leading zeros, plus signs, fractional parts or exponential forms.
+## Not acceptable forms are: "00", "-0", "+0", "+123", "0123", "7.0", "1E9".
+##
+## When a numeric value is passed, it must auto-convert into an acceptable
+## form. We recommend to not pass numeric values greater that 2^31-1.
+## An CONVFMT="%.0f" can help for larger numbers, but the author has
+## seen buggy awk-implementations for numbers around 2^31 or 2^32.
+##
+## The values returned from (most of) the functions here are strings.
+## NB: Do not compare such (string) results with the builtin operators!
+##
+## NB: Normally, numeric operations on uninitialized awk-variables
+## would auto-convert the ""-value into a numeric zero.
+## Currently, we here do not accept "" as a numeric zero. FFS.
+##-----------------------------------------------------------------------------
+## Internally of this "bignum" module we operate with unsigned numbers
+## of arbitrary length (U-numbers). U-numbers are coded as strings,
+## where a zero is coded as "", and no leading "0"s are permitted.
+## The functions that operate on these ## numbers shall not be called
+## from outside this module.
+##-----------------------------------------------------------------------------
+## This module is designed for correctness & maintainability, not for speed.
+## It can be used such that a lot of tests are performed during runtime
+## to assure proper operation.
+##
+## Section Overview:
+## - Condition Setup
+## - Calibration Setup
+## - Checking Support
+## - Internal Helper Functions
+## - Unsigned Arithmetic
+## - Signed Arithmetic
+## - Exported Arithmetic
+##
+## This module is an "experimental" first version. Therefore it is larger
+## than necessary. A correct version could be much shorter.
+##-----------------------------------------------------------------------------
+## Export Interface:
+##
+## bignum_init()
+## bnSetup(maxl)
+## Initializes module. Required. "maxl" is the number of decimal
+## digits this module shall use at most for the builtin arithmetic.
+## Typical is 9, 0 translates into the safe default 4.
+##
+## bn0(x) converts "" (uninitialized var) into a valid bignum zero
+## bnIsOk(x) whether is a valid bignum (quiet test)
+## bnCheck(x) whether is a valid bignum (non-quiet error reporting)
+##
+## bnSign(x) numeric: x<0: -1 x==0: 0 x>0: 1
+## bnCmp(x,y) numeric: x<y: -1 x==y: 0 x>y: 1
+##
+## bnLT(x,y) bool: x < y
+## bnLE(x,y) bool: x <= y
+## bnGE(x,y) bool: x >= y
+## bnGT(x,y) bool: x > y
+## bnEQ(x,y) bool: x == y
+## bnNE(x,y) bool: x != y
+##
+## bnLT0(x) bool: x < 0
+## bnLE0(x) bool: x <= 0
+## bnGE0(x) bool: x >= 0
+## bnGT0(x) bool: x > 0
+## bnEQ0(x) bool: x == 0
+## bnNE0(x) bool: x != 0
+##
+## bnNeg(x) -x
+## bnAbs(x) |x|
+##
+## bnInc(x) x + 1
+## bnDec(x) x - 1
+## bnAdd(x,y) x + y
+## bnSub(x,y) x - y
+## bnMul(x,y) x * y
+## bnMul2(x) x * 2
+## bnSqr(x) x * x
+##
+## bnDiv(x,y) floor(x / y)
+## bnMod(x,y) x mod y
+## bnGcd(x,y) gcd(x,y) FFS: precise definition
+##
+## bnLog10(x) log10(|x|)
+
+##=============================================================================
+## Condition Setup Section
+##
+## During runtime we will (conditionally) check many things:
+## s: Check missing Setup
+## S: (quiet) Self test after calibration
+## y: check sanity even at start of internal functions
+## I: check validity of input numbers at export interface
+## O: check validity of output numbers at export interface
+## i: check validity of input numbers at internal functions
+## o: check validity of output numbers at internal functions
+## f: errors are fatal
+## All: "fsSyioIO"
+##
+## By default, all error checking is enabled. To gain more speed,
+## error conditions can be turned off.
+##
+## bixDont[] contains the strings (indices), for which checking is *disabled*
+## FFS: export interface
+
+function bix_CheckAll( k) ## enable all checks/conditions
+{
+ for( k in bixDont ) delete bixDont[k]
+}
+
+function bix_CheckNot(conds ,c) ## diable some checks/conditions
+{
+ while( conds ) {
+ c = substr(conds,1,1)
+ conds = substr(conds,2)
+ bixDont[c] = 1
+ }
+}
+
+function bix_C_signature( r,k,a,i)
+{
+ r = "DONT:"
+ a = "fsSyioIO" ## serialize like this
+ for( i=1 ; i<=length(a) ; ++i ) {
+ k = substr(a,i,1)
+ if( k in bixDont ) {
+ r = r " " k
+ }
+ }
+ for( k in bixDont ) {
+ if( (length(k) != 1) || ! index(a,k) ) {
+ r = r " " k
+ }
+ }
+ return r ";"
+}
+
+##=============================================================================
+## Calibration Setup Section
+##
+## - Determine the largest absolute value which can be incremented and
+## decremented exactly by the builtin functions.
+## - Determine a length for operands of addition/multiplication, s.t.
+## the builtin operation is exact.
+
+function bix_conv_ok(s) # does not change when auto-converted
+{
+ if( ! ( (s "") == ("" (0 + s)) ) ) {
+ if( bix_optv ) printf("\nbix_conv: %s --> %s\n", s, "" (0+s))
+ return 0
+ }
+ return 1
+}
+
+function bix_incdec_ok(s ,d,t,u,v) # ++/-- are exact by builtin
+{
+ if( ! bix_conv_ok(s) ) return 0
+ if( ! bix_conv_ok(s+1) ) return 0
+ if( ! bix_conv_ok(s-1) ) return 0
+ for( d= -1 ; d<2 ; ++d ) {
+ t = s + d
+ u = t - s
+ v = s - t
+ if( ("" u) != ("" d) ) return 0
+ if( ("" v) != ("" 0-d) ) return 0
+ }
+ return 1
+}
+
+function bix_id_declen(neg,upto ,s,len,lenok)
+{
+ #FFS: further arg to prepend "1" instead of "": add-length
+ s = (neg ? "-" : "")
+ lenok = 0
+ for( len=1 ; len<=upto ; ++len ) {
+ s = "" s "9" ## gawk-fix
+ if( ! bix_incdec_ok(s) ) break
+ lenok = len
+ }
+ return lenok + (neg ? 1 : 0)
+}
+
+##-----------------------------------------------------------------------------
+## bixLinc length, ok for ++/-- builtin operands
+## bixLUinc length, ok for ++/-- builtin nonnegative operands
+## bixLadd length, ok for +/- builtin operands
+## bixLUadd length, ok for +/- builtin nonnegative operands
+## Addition with same sign: length can increase by one (carry)
+## Addition with different sign: magnitude cannot increase, but an
+## additional sign can be necessary, e.g.: "100" - "200" = "-100"
+## bixLmul length, ok for * builtin operands
+## bixLUmul length, ok for * builtin nonnegative operands
+
+function bix_L_signature()
+{
+ return sprintf("S++:%d U++:%d S+:%d U+:%d S*:%d U*:%d", \
+ bixLinc, bixLUinc, bixLadd, bixLUadd, bixLmul, bixLUmul)
+}
+
+function bix_E_setup(msg)
+{
+ printf("\n*** bignum: setup error: %s\n", msg)
+}
+
+function bn_Setup(maxl, ok,posl,negl,s,i,t1,t2) # --> whether ok
+{
+ ok = 1
+ maxl = int(maxl)
+ if( maxl <= 0 ) maxl = 4 ## safe default
+ posl = bix_id_declen(0,20)
+ negl = bix_id_declen(1,20)
+ bixLUinc = posl
+ bixLinc = ((posl < negl) ? posl : negl)
+ if( bixLUinc > 9 ) bixLUinc = 9 ## avoid 32-bit implementation bug
+ if( bixLinc > 9 ) bixLinc = 9 ## avoid 32-bit implementation bug
+ if( bixLUinc > maxl ) bixLUinc = maxl
+ if( bixLinc > maxl ) bixLinc = maxl
+ if( bixLUinc < 2 ) { bixLUinc = 2; ok = 0 }
+ if( bixLinc < 2 ) { bixLinc = 2; ok = 0 }
+ bixLUadd = bixLUinc - 1
+ bixLadd = bixLinc - 1
+ bixLUmul = int(bixLUadd/2)
+ bixLmul = int(bixLadd /2)
+ if( bixLUmul < 1 ) bixLUmul = 1
+ if( bixLmul < 1 ) bixLmul = 1
+ ## Check some builtin properties which we need ...
+ s = "0013" # must not be octal, but decimal
+ if( (0 + s) != 13 ) {
+ ok = 0
+ bix_E_setup("0013 is not converted decimally")
+ }
+ ## Check length-1 string compare to be compatible with arithmetic
+ s = "0123456789"
+ for( i=1 ; i<10 ; ++i ) {
+ t1 = substr(s,i ,1)
+ t2 = substr(s,i+1,1)
+ if( ! (t1 < t2) ) {
+ ok = 0
+ bix_E_setup( sprintf("\"%s\" < \"%s\" failed",t1,t2) )
+ }
+ }
+ if( ! ((s s s s "6") < (s s s s "7")) ) {
+ ok = 0
+ bix_E_setup( "long string compare failed" )
+ }
+ if( sprintf("%0*d", 9, 0) != "000000000" ) {
+ ok = 0
+ bix_E_setup( "sprintf(\"%0*d\") failed" )
+ }
+ if( ! ok ) ++bix_errs
+ bix_cal_done = 1
+ ## self test ...
+ if( ! bix_Test(0) ) ok = 0
+ return ok
+}
+
+function bnSetup(maxl) ## --> whether ok; tries to die on error
+{
+ if( ! bn_Setup(maxl) ) {
+ printf("\n*** bignum: bnSetup failed.\n")
+ bix_show_setup()
+ ## This error is always fatal!
+ do_exit(1); exit(1)
+ return 0 ## huh, did not really exit
+ }
+ return 1 ## ok
+}
+
+##=============================================================================
+## Self Test Section
+
+function bix__tst_add(a,b,c ,x,y,ok,z)
+{
+ bnCheck(a); bnCheck(b); bnCheck(c)
+ if( bix_vtst ) printf("Calling bnAdd(%s,%s)\n", a,b)
+ x = bnAdd(a,b)
+ y = "" (a+b)
+ ok = (x==c)
+ if( ! ok || bix_vtst ) {
+ printf("%s + %s\n --> %s [%s]\n ==> %s %s\n", a,b, x,y,c, ok?"OK":"BAD")
+ }
+ if( ! ok ) ++bix_errs
+ z = bnSub(a, bnNeg(b))
+ ok = (z==c)
+ if( ! ok ) ++bix_errs
+ if( ! ok && bix_vtst ) {
+ printf("bnSub(%s, bnNeg(%s))\n --> %s BAD\n", a,b,z)
+ }
+}
+
+function bix__tstAdd( s,s0,t)
+{
+ bix__tst_add("1234", "769", "2003")
+ bix__tst_add( 1234 , 769 , 2003 )
+ bix__tst_add( "9999999", "99999999", "109999998")
+ bix__tst_add("-9999999", "99999999", "90000000")
+ bix__tst_add( "9999999", "-99999999", "-90000000")
+ s = "9999999999"
+ s0 = "0000000000"
+ t = "999999999"
+ bix__tst_add( s s s, s s s, "1" s s t "8")
+ bix__tst_add( s s s, "-" s s s, "0")
+ bix__tst_add("-" s s s, "-" s s s, "-1" s s t "8")
+ bix__tst_add("-" s s s, "-" s s, "-1" s0 s t "8")
+ bix__tst_add( s s s, "-" s s, s s0 s0)
+ bix__tst_add( s s "3", "-" s s "1", "2")
+ bix__tst_add( s s "1", "-" s s "3", "-2")
+ bix__tst_add("1" s0 s0, "-" s s , "1")
+ bix__tst_add("1" s0 s0 "34", "-" s s "32", "102")
+}
+
+function bix__tst_mul(a,b,c ,x,y,ok)
+{
+ bnCheck(a); bnCheck(b); bnCheck(c)
+ if( bix_vtst ) printf("Calling bnMul(%s,%s)\n", a,b)
+ x = bnMul(a,b)
+ y = "" (a*b)
+ ok = (x==c)
+ if( ! ok || bix_vtst ) {
+ printf("%s * %s\n --> %s [%s]\n ==> %s %s\n", a,b,x,y, c, ok?"OK":"BAD")
+ }
+ if( ! ok ) ++bix_errs
+}
+
+function bix__tst_4mul(a,b,c)
+{
+ bix__tst_mul( a , b , c )
+ bix__tst_mul( a ,bnNeg(b),bnNeg(c))
+ bix__tst_mul(bnNeg(a), b ,bnNeg(c))
+ bix__tst_mul(bnNeg(a),bnNeg(b), c )
+}
+
+function bix__tstMul( s,s0,t)
+{
+ bix__tst_4mul(1,1,1)
+ bix__tst_4mul(0,0,0)
+ bix__tst_4mul(0,1,0)
+ bix__tst_4mul(0,999,0)
+ bix__tst_4mul(1,10,10)
+ bix__tst_4mul(9,1,9)
+ bix__tst_4mul(8,8,64)
+ bix__tst_4mul("1024","1024","1048576")
+ bix__tst_4mul("1234", "769", "948946")
+ s = "9999999999"
+ s0 = "0000000000"
+ t = "999999999"
+ bix__tst_4mul(2 s0, 9 s0, 18 s0 s0)
+ bix__tst_4mul("3468795187645","8273582372564", "28699362718534504679771780")
+ bix__tst_4mul(1, s s, s s)
+}
+
+function bix__tst_sqr(a,c ,x,y,ok)
+{
+ bnCheck(a); bnCheck(c)
+ if( bix_vtst ) printf("Calling bnSqr(%s)\n", a)
+ x = bnSqr(a)
+ y = "" (a*a)
+ ok = (x==c)
+ if( ! ok || bix_vtst ) {
+ printf("sqr(%s)\n --> %s [%s]\n ==> %s %s\n", a,x,y, c, ok?"OK":"BAD")
+ }
+ if( ! ok ) ++bix_errs
+}
+
+function bix__tst_2sqr(a,c)
+{
+ bix__tst_sqr( a , c)
+ bix__tst_sqr(bnNeg(a), c)
+}
+
+function bix__tstSqr( s,t)
+{
+ bix__tst_2sqr(0,0)
+ bix__tst_2sqr(1,1)
+ bix__tst_2sqr(5,25)
+ bix__tst_2sqr(99,9801)
+ bix__tst_2sqr(100,"10000")
+ t = "0000"
+ bix__tst_2sqr(1 t t, 1 t t t t)
+ t = t t
+ bix__tst_2sqr(1 t t, 1 t t t t)
+ bix__tst_2sqr("3468795187645", "12032540053829110760646025")
+ bix__tst_2sqr("8273582372564", "68452165275601747299934096")
+ s = "9999999999"
+ bix__tst_2sqr( s, "99999999980000000001")
+ bix__tst_2sqr(s s, "9999999999999999999800000000000000000001")
+}
+
+function bix__tst_divmod(a,b,c,d ,x,y,ok)
+{
+ bnCheck(a); bnCheck(b); bnCheck(c); bnCheck(d)
+ if( bix_vtst ) printf("Calling bnDiv/bnMod(%s,%s)\n", a,b)
+ x = bnDiv(a,b)
+ y = "" int(a/b)
+ ok = (x==c)
+ if( ! ok || bix_vtst ) {
+ printf("%s / %s\n --> %s [%s]\n ==> %s %s\n", a,b,x,y, c, ok?"OK":"BAD")
+ }
+ if( ! ok ) ++bix_errs
+ #
+ x = bnMod(a,b)
+ y = "" int(a%b)
+ ok = (x==d)
+ if( ! ok || bix_vtst ) {
+ printf("%s / %s\n --> %s [%s]\n ==> %s %s\n", a,b,x,y, d, ok?"OK":"BAD")
+ }
+ if( ! ok ) ++bix_errs
+}
+
+function bix__tstDivMod( s,s0,t)
+{
+ s = "9999999999"
+ s0 = "0000000000"
+ t = "999999999"
+ bix__tst_divmod(0,1, 0,0)
+ bix__tst_divmod(0,s, 0,0)
+ bix__tst_divmod(s,1, s,0)
+ bix__tst_divmod(20,3, 6,2)
+ bix__tst_divmod(-20,3, -7,1)
+ bix__tst_divmod(-21,3, -7,0)
+ bix__tst_divmod(s,t, 10,9)
+ bix__tst_divmod("-" s,t, -11,substr(t,2)"0")
+ bix__tst_divmod("959682618430873425",2, "479841309215436712",1)
+ bix__tst_divmod("48937568913476009","10092784", "4848768081","5848505")
+ #bix__tst_divmod(-1,-1, 1,0)
+}
+
+function bix__tst_gcd(a,b,c ,x,ok)
+{
+ bnCheck(a); bnCheck(b); bnCheck(c)
+ if( bix_vtst ) printf("Calling bnGcd(%s,%s)\n", a,b)
+ x = bnGcd(a,b)
+ ok = (x==c)
+ if( ! ok || bix_vtst ) {
+ printf("gcd(%s, %s)\n --> %s\n ==> %s %s\n", a,b,x, c, ok?"OK":"BAD")
+ }
+ if( ! ok ) ++bix_errs
+}
+
+function bix__tst_4gcd(a,b,c)
+{
+ bix__tst_gcd(a,b,c)
+ if( (a!="0") && (b=="0") ) bix__tst_gcd("-" a, b, "-" c)
+ if( (a!="0") && (b!="0") ) bix__tst_gcd("-" a, b, c)
+ if( b!="0" ) bix__tst_gcd( a, "-" b, "-" c)
+ if( (a!="0") && (b!="0") ) bix__tst_gcd("-" a, "-" b, "-" c)
+}
+
+function bix__tst_8gcd(a,b,c)
+{
+ bix__tst_4gcd(a,b,c)
+ if( a != b ) {
+ bix__tst_4gcd(b,a,c)
+ }
+}
+
+function bix__tstGcd( s)
+{
+ bix__tst_8gcd("33","0","33")
+ bix__tst_8gcd("1","1","1")
+ bix__tst_8gcd("7","7","7")
+ bix__tst_8gcd("6","35","1")
+ bix__tst_8gcd("6","14","2")
+ bix__tst_8gcd("60","100","20")
+ s = "999999999"
+ bix__tst_8gcd(s, "27", "27")
+ bix__tst_8gcd(s s, s, s)
+ ## two fibo-numbers
+ bix__tst_8gcd( "308061521170129", "498454011879264", "1")
+}
+
+function bix_Test(verbose ,ov,oe) # --> whether ok
+{
+ oe = bix_errs
+ ov = bix_vtst; bix_vtst = bix_vtst || verbose
+ bix__tstAdd()
+ bix__tstMul()
+ bix__tstSqr()
+ bix__tstDivMod()
+ bix__tstGcd()
+ bix_vtst = ov
+ return bix_errs == oe
+}
+
+##=============================================================================
+## Checking Support Section
+
+function bix_show_setup( k)
+{
+ printf("# bix_cal_done = %d\n", bix_cal_done)
+ printf("# bixLUinc = %d\n", bixLUinc)
+ printf("# bixLinc = %d\n", bixLinc )
+ printf("# bixLUadd = %d\n", bixLUadd)
+ printf("# bixLadd = %d\n", bixLadd )
+ printf("# bixLUmul = %d\n", bixLUmul)
+ printf("# bixLmul = %d\n", bixLmul )
+ printf("# bix_errs = %d\n", bix_errs)
+ for( k in bixDont ) {
+ printf("# bixDont[ %s ]\n", k)
+ }
+}
+
+function bix_CalLok(len)
+{
+ return (len >= 1) && (len <= 9) && (len == int(len))
+}
+
+function bix_CalOk()
+{
+ if( bix_cal_done != 1 ) return 0
+ if( ! bix_CalLok(bixLUinc) ) return 0
+ if( ! bix_CalLok(bixLinc ) ) return 0
+ if( ! bix_CalLok(bixLUadd) ) return 0
+ if( ! bix_CalLok(bixLadd ) ) return 0
+ if( ! bix_CalLok(bixLUmul) ) return 0
+ if( ! bix_CalLok(bixLmul ) ) return 0
+ return 1
+}
+
+function bix_Sanity() ## check Setup & internal calibration
+{
+ if( ! bix_CalOk() ) {
+ printf("\n*** bignum: missing or bad setup\n")
+ bix_show_setup()
+ ## This error is always fatal!
+ do_exit(1); exit(1)
+ return 0 ## huh, did not really exit
+ }
+ return 1 ## ok
+}
+
+##.............................................................................
+## Checks for permissable values
+##
+## bn_IsI whether it is a valid bignum (I)
+## bn_IsU whether it is a valid (internal) unsigned bignum (U)
+## bn_IsC whether it is a valid comparison result
+
+function bn_IsI(x) # --> whether is a valid bignum
+{
+ if( x == "0" ) return 1
+ if( x ~ /^-?[1-9][0-9]*$/ ) return 1
+ return 0
+}
+
+function bn_IsU(x) # --> whether is a valid internal unsigned
+{
+ if( x == "" ) return 1
+ if( x ~ /^[1-9][0-9]*$/ ) return 1
+ return 0
+}
+
+function bn_IsC(x) # --> whether is a valid comparison result
+{
+ if( x == "-1" ) return 1
+ if( x == "1" ) return 1
+ if( x == "0" ) return 1
+ return 0
+}
+
+function bix_Err(fatal) # --> 0 (not ok)
+{
+ ++bix_errs
+ if( fatal || ! ("f" in bixDont) ) {
+ printf("*** bignum: terminating.\n")
+ do_exit(1); exit(1)
+ }
+ return 0 ## sorry
+}
+
+function bix_Cry(what, fname, v)
+{
+ printf("\n*** bignum: %s of \"%s\" is not valid: \"%s\"\n", what, fname, v)
+}
+
+##.............................................................................
+## Checks for valid input operands of external (interface) functions
+## Can be disabled by bixDont["I"]
+## Return on error only for "f" in bixDont[].
+
+function bix_I1(fname,x) # --> ok
+{
+ if( ! bn_IsI(x) ) { bix_Cry("arg", fname, x); return bix_Err(0) }
+ return 1 ## ok
+}
+
+function bix_I2(fname,x,y) # --> ok
+{
+ if( ! bn_IsI(x) ) { bix_Cry("arg1", fname, x); return bix_Err(0) }
+ if( ! bn_IsI(y) ) { bix_Cry("arg2", fname, y); return bix_Err(0) }
+ return 1 ## ok
+}
+
+##.............................................................................
+## Checks for valid input operands of internal functions
+## Can be disabled by bixDont["i"]
+## Return on error only for "f" in bixDont[].
+
+function bix_i1(fname,x) # --> ok
+{
+ if( ! bn_IsU(x) ) { bix_Cry("arg", fname, x); return bix_Err(0) }
+ return 1 ## ok
+}
+
+function bix_i2(fname,x,y) # --> ok
+{
+ if( ! bn_IsU(x) ) { bix_Cry("arg1", fname, x); return bix_Err(0) }
+ if( ! bn_IsU(y) ) { bix_Cry("arg2", fname, y); return bix_Err(0) }
+ return 1 ## ok
+}
+
+function bix_iI(fname,x,y) # --> ok
+{
+ if( ! bn_IsU(x) ) { bix_Cry("arg1", fname, x); return bix_Err(0) }
+ if( ! bn_IsI(y) ) { bix_Cry("arg2", fname, y); return bix_Err(0) }
+ return 1 ## ok
+}
+
+##.............................................................................
+## Function for output operand validation.
+##
+## Arithmetic functions are obliged to use them at every "return" statement.
+## Instead of: return x
+## They call e.g.: return bix_O1("myname", x)
+##
+## For external (interface) functions, "O" in bixDont disables the check.
+## For internal functions, "o" in bixDont disables the check.
+
+function bix_OC(fname,x) # --> x
+{
+ if( ! ("O" in bixDont) ) {
+ if( ! bn_IsC(x) ) { bix_Cry("res", fname, x); bix_Err(0) }
+ }
+ return x ## anyhow
+}
+
+function bix_O1(fname,x) # --> x
+{
+ if( ! ("O" in bixDont) ) {
+ if( ! bn_IsI(x) ) { bix_Cry("res", fname, x); bix_Err(0) }
+ }
+ return x ## anyhow
+}
+
+function bix_O01(fname,arr)
+{
+ if( ! ("O" in bixDont) ) {
+ if( ! bn_IsI(arr[0]) ) { bix_Cry("res[0]", fname, arr[0]); bix_Err(0) }
+ if( ! bn_IsI(arr[1]) ) { bix_Cry("res[1]", fname, arr[1]); bix_Err(0) }
+ }
+}
+
+function bix_o1(fname,x) # --> x
+{
+ if( ! ("o" in bixDont) ) {
+ if( ! bn_IsU(x) ) { bix_Cry("res", fname, x); bix_Err(0) }
+ }
+ return x ## anyhow
+}
+
+function bix_o01(fname,arr)
+{
+ if( ! ("o" in bixDont) ) {
+ if( ! bn_IsU(arr[0]) ) { bix_Cry("res[0]", fname, arr[0]); bix_Err(0) }
+ if( ! bn_IsU(arr[1]) ) { bix_Cry("res[1]", fname, arr[1]); bix_Err(0) }
+ }
+}
+
+##=============================================================================
+## Non-exported help functions
+
+function bn_nR(x) ## conversion of small numeric results
+{
+ #return "" x ## everything is a string
+ return x ## small values stay numeric
+ #return (length(x) <= bixLinc) ? 0+x : "" x
+}
+
+function bn_U(x) # numeric --> string (U-number)
+{
+ return x ? ("" x) : "" # convert numeric 0 to ""
+}
+
+function bn_mksc(sign) # --> make sign character from numeric value
+{
+ return (sign < 0) ? "-" : ""
+}
+
+function bn_repstr(str,cnt)
+{
+ if( cnt <= 2 ) {
+ if( cnt < 1 ) return ""
+ if( cnt < 2 ) return str
+ return str str
+ }
+ return bn_repstr(str str, int(cnt/2)) ((cnt%2) ? str : "")
+ ## FFS: iterate instead of recursion
+}
+function bn_00(len)
+{
+ if( len <= 0 ) return ""
+ if( len <= 99 ) { ## sprintf is not safe for large lengthes!
+ return sprintf("%0*d", len, 0)
+ }
+ return bn_repstr("0", len)
+}
+
+function bn_0fill(x,len) ## --> string, at least length len, 0-filled
+{
+ if( len > length(x) ) return bn_00(len - length(x)) x
+ return "" x
+}
+
+function bn_conc(hi,lo,lolen) ## --> string: lo + hi*10^lolen
+{
+ ## hi is string, may be ""
+ ## lo is numeric, at most of length lolen
+ ## Result is "" for zero
+ if( hi == "" ) return bn_U(lo)
+ return hi bn_0fill(lo, lolen)
+}
+
+function bn_hinum(num,len) ## --> numeric hi part (split size len)
+{
+ ## num is string (or numeric), of length at most 2*len
+ if( length(num) <= len ) return 0
+ return 0 + substr(num, 1, length(num)-len)
+}
+
+function bn_lostr(num,len) ## --> string lo part, zero filled
+{
+ ## num is string (or numeric), of length at most 2*len
+ if( length(num) <= len ) return bn_0fill(num, len)
+ return substr(num, 1+ length(num)-len)
+}
+
+##-----------------------------------------------------------------------------
+## Helper functions for "digit" arrays
+##
+## When we build up a loner result it happens that we collect parts
+## of it in an array. Increasing indexes indicated higher significance,
+## i.e. larger indexes end up at the left end of the number.
+## Array contents may be strings of numbers. Logically the are some
+## fixed amount of decimal digits apart, fixed per array.
+##
+## Indexes "lox" and "hix", selecting a part of an array, are both inclusive.
+## FFS: 0-fill in highest element, also?
+
+function bn__arr_0fill(arr,lox,hix,dL ,i,v,vL)
+{
+ ## Force all elements to have size at least "dL" by filling high zeros
+ for( i=lox ; i<=hix ; ++i ) {
+ v = arr[i]
+ vL = length(v)
+ if( vL < dL ) {
+ arr[i] = bn_00(dL-vL) v
+ }
+ }
+}
+
+function bn__arr_pairs(arr,lox,hix ,i,z) ## --> new hix
+{
+ if( lox < hix ) { ## at least one pair present
+ z = lox - 1
+ for( i=lox ; i<hix ; i+=2 ) { ## have another pair
+ arr[++z] = arr[i+1] arr[i]
+ }
+ if( i == hix ) { ## single rest stays single
+ arr[++z] = arr[i]
+ }
+ ## FFS: delete arr[z+1 .. hix]
+ return z
+ }
+ return hix
+}
+
+function bn__arr_conc_triv(arr,lox,hix ,i,r) ## --> concatenated value
+{
+ if( lox > hix ) return ""
+ if( lox == hix ) return arr[lox]
+ i = lox
+ r = arr[i]
+ while( ++i <= hix ) {
+ r = arr[i] r ## prepend (!) next arr-elem
+ }
+ return r
+}
+
+function bn__arr_concL_triv(arr,lox,hix,dL ,i,r)
+{
+ r = ""
+ for( i=lox ; i<=hix ; ++i ) {
+ r = bn_0fill(arr[i],dL) r ## prepend (!) next arr-elem
+ }
+ return r
+}
+
+function bn__arr_conc(arr,lox,hix)
+{
+ if( ! (lox <= hix) ) return ""
+ if( (hix-lox+1) >= 4 ) { ## FFS: more than 4?
+ while( lox < hix ) {
+ hix = bn__arr_pairs(arr,lox,hix)
+ }
+ return arr[lox]
+ }
+ return bn__arr_conc_triv(arr,lox,hix)
+}
+
+function bn__arr_concL(arr,lox,hix,dL)
+{
+ bn__arr_0fill( arr,lox,hix,dL)
+ return bn__arr_conc(arr,lox,hix )
+}
+
+##-----------------------------------------------------------------------------
+
+function bn_Uadd(x,y,carry ,xL,yL,lo,xhi,yhi,r) # --> "" x+y+carry
+{
+ ## x,y are strings, 0 comes in (and out) as ""
+ ## 0 <= carry <= 1 comes in numeric
+ #printf("bn_Uadd(%s,%s,%s)\n",x,y,carry)
+ if( ! ("y" in bixDont) ) bix_Sanity()
+ if( ! ("i" in bixDont) ) bix_i2("bn_Uadd",x,y)
+
+ xL = length(x)
+ if( !xL ) return bix_o1("bn_Uadd", carry ? bn_Uadd(y, "" carry, 0) : y)
+ yL = length(y)
+ if( !yL ) return bix_o1("bn_Uadd", carry ? bn_Uadd(x, "" carry, 0) : x)
+ ## neither is empty
+ if( (xL <= bixLUadd) && (yL <= bixLUadd) ) {
+ return bix_o1("bn_Uadd", bn_U(x + y + carry))
+ }
+ ## At least one operand is long
+ if( 0 ) { ## old version
+ ## Split in lo and hi part, summing lo part by builtin
+ lo = carry
+ if( xL > bixLUadd ) {
+ lo += substr(x, 1+ xL-bixLUadd)
+ xhi = substr(x, 1, xL-bixLUadd)
+ }else {
+ lo += x
+ xhi = ""
+ }
+ if( yL > bixLUadd ) {
+ lo += substr(y, 1+ yL-bixLUadd)
+ yhi = substr(y, 1, yL-bixLUadd)
+ }else {
+ lo += y
+ yhi = ""
+ }
+ #printf("bn_Uadd(%s,%s,%s):xhi=%s,yhi=%s,lo=%s\n",x,y,carry,xhi,yhi,lo)
+ #lo = "" lo
+ r = bn_Uadd(xhi, yhi, bn_hinum(lo,bixLUadd)) bn_lostr(lo,bixLUadd)
+ return bix_o1("bn_Uadd", r)
+ }else { ## new version
+ ## FFS: sequentially building up the result can still be quadratic
+ r = ""
+ while( (xL >= bixLUadd) && (yL >= bixLUadd) ) {
+ ## both have a full digit
+ lo = carry
+ lo += substr(x, 1+ xL-bixLUadd, bixLUadd)
+ lo += substr(y, 1+ yL-bixLUadd, bixLUadd)
+ carry = bn_hinum(lo,bixLUadd)
+ r = bn_lostr(lo,bixLUadd) r ## prepend to result r
+ xL -= bixLUadd
+ yL -= bixLUadd
+ }
+ ## At least one length is smaller than bixLUadd
+ while( (xL && yL) || (carry && (xL || yL)) ) {
+ lo = carry
+ if( xL > bixLUadd ) {
+ lo += substr(x, 1+ xL-bixLUadd, bixLUadd)
+ xL -= bixLUadd
+ }else if( xL ) {
+ lo += substr(x, 1, xL)
+ xL = 0
+ }
+ if( yL > bixLUadd ) {
+ lo += substr(y, 1+ yL-bixLUadd, bixLUadd)
+ yL -= bixLUadd
+ }else if( yL ) {
+ lo += substr(y, 1, yL)
+ yL = 0
+ }
+ if( xL || yL || (length(lo) > bixLUadd) ) { ## more hi stuff
+ carry = bn_hinum(lo,bixLUadd)
+ r = bn_lostr(lo,bixLUadd) r
+ }else {
+ r = lo r
+ carry = 0
+ }
+ }
+ ## At most one of xL, yL and carry is left non-zero
+ if( xL ) {
+ r = substr(x, 1, xL) r
+ }else if( yL ) {
+ r = substr(y, 1, yL) r
+ }else if( carry ) {
+ r = carry r
+ }
+ return bix_o1("bn_Uadd", r)
+ }
+}
+
+function bn_Umul2(x,carry ,xL,lo,r) # --> "" x+x+carry
+{
+ ## 0 <= carry <= 1 comes in numeric
+ ## We implement (2*x) as (x+x): this is a reduced version of bn_Uadd.
+ #printf("bn_Umul2(%s,%s)\n",x,carry)
+ if( ! ("y" in bixDont) ) bix_Sanity()
+ if( ! ("i" in bixDont) ) bix_i1("bn_Umul2",x)
+
+ xL = length(x)
+ if( !xL ) return bix_o1("bn_Umul2", bn_U( carry))
+ if( xL <= bixLUadd ) return bix_o1("bn_Umul2", bn_U(x + x + carry))
+ ## x is long
+ ## FFS: sequentially building up the result can still be quadratic
+ r = ""
+ while( xL >= bixLUadd ) { ## yet another full digit there
+ lo = carry + 2*substr(x, 1+ xL-bixLUadd, bixLUadd)
+ carry = bn_hinum(lo,bixLUadd)
+ r = bn_lostr(lo,bixLUadd) r ## prepend to result "r"
+ xL -= bixLUadd
+ }
+ ## 0 <= xL < bixLUadd
+ if( xL ) {
+ r = (carry + 2*substr(x, 1, xL)) r
+ }else if( carry ) {
+ r = carry r
+ }
+ return bix_o1("bn_Umul2", r)
+}
+
+function bn_Usub(x,y,carry ,xL,yL,lo,hi,xhi,yhi,k,arr) # --> "" x-y-carry
+{
+ ## x,y are strings, 0 comes in (and out) as ""
+ ## 0 <= carry <= 1 comes in numeric
+ ## Caller guarantees, that the result is non-negative, i.e. x >= y+carry
+ if( ! ("y" in bixDont) ) bix_Sanity()
+ if( ! ("i" in bixDont) ) bix_i2("bn_Usub",x,y)
+ #printf("bn_Usub(%s,%s,%s)[%d]\n",x,y,carry,bixLUadd)
+
+ xL = length(x)
+ yL = length(y) # xL >= yL
+ if( ! yL ) return bix_o1("bn_Usub", carry ? bn_Usub(x, "" carry, 0) : x)
+ if( !carry && (xL == yL) && (x == y) ) {
+ return bix_o1("bn_Usub", "") # x-x --> ""
+ }
+ if( xL <= bixLUadd ) return bix_o1("bn_Usub", bn_U(x-y-carry))
+ ## At least one is long
+ if( 0 ) { ## old version
+ ## Split in lo and hi part, summing lo part by builtin
+ lo = 0-carry
+ if( xL > bixLUadd ) {
+ lo += substr(x, 1+ xL-bixLUadd)
+ xhi = substr(x, 1, xL-bixLUadd)
+ }else {
+ lo += x
+ xhi = ""
+ }
+ if( yL > bixLUadd ) {
+ lo -= substr(y, 1+ yL-bixLUadd)
+ yhi = substr(y, 1, yL-bixLUadd)
+ }else {
+ lo -= y
+ yhi = ""
+ }
+ #printf("bn_Usub(%s,%s,%s):xhi=%s,yhi=%s,lo=%s\n",x,y,carry,xhi,yhi,lo)
+ carry = (lo < 0)
+ if( carry ) lo += ("1" bn_00(bixLUadd)) ## FFS
+ hi = bn_Usub(xhi, yhi, carry)
+ #printf("bn_Usub(%s,%s,%s) --> %s\n",xhi,yhi,carry,hi)
+ return bix_o1("bn_Usub", bn_conc(hi, lo, bixLUadd) )
+ }else { ## new version
+ k = 0
+ while( yL >= bixLUadd ) {
+ lo = 0-carry
+ lo += substr(x, 1+ xL-bixLUadd, bixLUadd)
+ lo -= substr(y, 1+ yL-bixLUadd, bixLUadd)
+ carry = (lo < 0)
+ if( carry ) lo += ("1" bn_00(bixLUadd)) ## FFS
+ arr[++k] = lo ## save numeric "digit"
+ xL -= bixLUadd
+ yL -= bixLUadd
+ }
+ ## 0 <= yL < bixLUadd; yL <= xL
+ while( (xL > 0) && (carry || (yL > 0)) ) {
+ lo = 0-carry
+ if( xL > bixLUadd ) {
+ lo += substr(x, 1+ xL-bixLUadd, bixLUadd)
+ xL -= bixLUadd
+ }else if( xL > 0 ) {
+ lo += substr(x, 1, xL)
+ xL = 0
+ }
+ if( yL > 0 ) {
+ lo -= substr(y, 1, yL)
+ yL = 0
+ }
+ carry = (lo < 0)
+ if( carry ) lo += ("1" bn_00(bixLUadd)) ## FFS
+ arr[++k] = lo ## save numeric "digit"
+ }
+ if( carry ) { bix_Cry("carry", "bn_Usub", carry); return bix_Err(1) }
+ ## Append hi part of x (if any left) as a "super-digit"
+ if( xL ) {
+ arr[++k] = substr(x, 1, xL)
+ }else {
+ ## Skip hi zeros
+ while( k && !arr[k] ) --k
+ ## Any digits left?
+ if( ! k ) return bix_o1("bn_Usub", "")
+ }
+ ## Now arr[k] is the highest part, and is non-zero
+ return bix_o1("bn_Usub", arr[k] bn__arr_concL(arr,1,k-1,bixLUadd) )
+ }
+}
+
+##-----------------------------------------------------------------------------
+## Multiplication
+## We do it absolutely straigthforward, like with paper & pencil.
+
+function bn_UmulDig(x,y,carry ,xL,xlo,xhi,lo,hi,dL,r)
+{ ## --> x*y + carry
+ ## x is string ("" for 0)
+ ## y is numeric, a multiplication digit (length <= bixLUmul)
+ ## carry is numeric (length <= bixLUmul)
+ if( ! ("y" in bixDont) ) bix_Sanity()
+ if( ! ("i" in bixDont) ) bix_iI("bn_UmulDig",x,y)
+ #printf("bn_UmulDig(%s,%s,%s)[%d]\n",x,y,carry,bixLUmul)
+
+ xL = length(x)
+ if( !xL || !y ) return bix_o1("bn_UmulDig", bn_U(carry) )
+ ## Neither operand is zero.
+ if( xL <= bixLUmul ) return bix_o1("bn_UmulDig", bn_U(x * y + carry) )
+ if( y == 1 ) return bix_o1("bn_UmulDig", bn_Uadd(x, "", carry) )
+ if( 0 ) { ## old version
+ xlo = substr(x, 1+ xL-bixLUmul)
+ xhi = substr(x, 1, xL-bixLUmul)
+ lo = xlo * y + carry
+ hi = bn_UmulDig(xhi, y, bn_hinum(lo,bixLUmul))
+ return bix_o1("bn_UmulDig", hi bn_lostr(lo,bixLUmul) )
+ }else { ## new version
+ dL = 2*bixLUmul - length(y)
+ ## We bite off pieces of length dL from x, and build the result
+ ## in pieces of the same size.
+ r = ""
+ while( xL > dL ) {
+ lo = substr(x, 1+ xL-dL, dL) * y + carry
+ carry = bn_hinum(lo, dL)
+ r = bn_lostr(lo, dL) r
+ xL -= dL
+ }
+ return bix_o1("bn_UmulDig", (substr(x, 1, xL) * y + carry) r )
+ }
+}
+
+function bn__arrADD(za,zx,z,dL ,zL) # --> (numeric) carry
+{
+ z += za[zx]
+ zL = length(z)
+ if( zL <= dL ) {
+ za[zx] = z
+ return 0
+ }else { # zL > dL
+ za[zx] = 0 + substr(z, 1+ zL-dL)
+ return 0 + substr(z, 1, zL-dL)
+ }
+}
+
+function bn__UmulADA(xa,xhix,y,za,zlox ,c,i)
+{
+ ## za[zlox..] += xa[0..xhix] * y
+ if( ! y ) return
+ c = 0
+ for( i=0 ; i<=xhix ; ++i ) {
+ c = bn__arrADD(za,zlox+i, xa[i]*y + c, bixLUmul)
+ }
+ while( c ) {
+ c = bn__arrADD(za,zlox+i, c, bixLUmul)
+ ++i
+ }
+}
+
+function bn__arr_nsplit(x,dL,xarr,xlox ,xL) # --> (new) xhix
+{
+ xL = length(x)
+ while( xL >= dL ) {
+ xarr[xlox++] = 0 + substr(x, 1+ xL-dL, dL)
+ xL -= dL
+ }
+ if( xL ) {
+ xarr[xlox++] = 0 + substr(x, 1, xL)
+ }
+ return xlox - 1
+}
+
+function bn_Umul(x,y ,xL,yL,ylo,yhi,lo,hi,xa,xhix,za,zx)
+{
+ if( ! ("y" in bixDont) ) bix_Sanity()
+ if( ! ("i" in bixDont) ) bix_i2("bn_Umul",x,y)
+ #printf("bn_Umul(%s,%s)[%d]\n",x,y,bixLUmul)
+
+ xL = length(x)
+ yL = length(y)
+ if( !xL || !yL ) return bix_o1("bn_Umul", "" )
+ if( xL <= bixLUmul ) return bix_o1("bn_Umul", bn_UmulDig(y, 0+x, 0) )
+ if( yL <= bixLUmul ) return bix_o1("bn_Umul", bn_UmulDig(x, 0+y, 0) )
+ if( xL < yL ) return bix_o1("bn_Umul", bn_Umul(y,x) )
+ if( 0 ) { ## old version
+ ylo = substr(y, 1+ yL-bixLUmul)
+ yhi = substr(y, 1, yL-bixLUmul)
+ hi = bn_Umul(x, yhi) bn_0fill(0,bixLUmul)
+ lo = bn_UmulDig(x, 0+ylo, 0)
+ #printf("bn_Umul(%s,%s)[%d] --> %s + %s\n",x,y,bixLUmul, hi,lo)
+ return bix_o1("bn_Umul", bn_Uadd(hi, lo, 0) )
+ }else { ## new version
+ #printf("UMUL{ xL=%d yL=%d\n", xL, yL)
+ xa[0] = 0 ## arrify
+ za[0] = 0 ## arrify
+ xhix = bn__arr_nsplit(x, bixLUmul, xa, 0)
+ #printf("UMUL s.\n")
+ #printf("bixLUmul=%d x=%s xhix=%d y=%s\n", bixLUmul, x, xhix, y)
+ #for( j=0 ; j<=xhix ; ++j) printf(" xa[%d] = %4s\n", j,xa[j])
+ for( zx=0 ; yL >= bixLUmul ; yL-=bixLUmul ) {
+ bn__UmulADA(xa,xhix, 0+substr(y, 1+ yL-bixLUmul, bixLUmul), za,zx++)
+ #for( j=0 ; (j in za) ; ++j) printf(" za[%d] = %4s\n", j,za[j])
+ #printf("UMUL y.\n")
+ }
+ if( yL ) {
+ bn__UmulADA(xa,xhix, 0+substr(y, 1, yL), za,zx++)
+ #printf("UMUL y\n")
+ }
+ #printf("UMUL Y\n")
+ zx += xhix
+ if( (zx+1) in za ) { bix_Cry("zx", "bn_Usub", zx); return bix_Err(1) }
+ ## Skip hi zeros
+ while( (zx>=0) && !za[zx] ) --zx
+ lo = bn__arr_concL(za,0,zx-1, bixLUmul)
+ #printf("UMUL c.\n")
+ #printf("}UMUL\n")
+ return bix_o1("bn_Umul", za[zx] lo )
+ }
+}
+
+function bn_Usqr(x ,xL,xa,xhix,za,zhix,j,k,c,xaj)
+{
+ if( ! ("y" in bixDont) ) bix_Sanity()
+ if( ! ("i" in bixDont) ) bix_i1("bn_Usqr",x)
+ #printf("bn_Usqr(%s)[%d]\n",x,bixLUmul)
+
+ xL = length(x)
+ if( !xL ) return bix_o1("bn_Usqr", "" )
+ if( xL <= bixLUmul ) return bix_o1("bn_Usqr", bn_U(x*x) )
+
+ #printf("USQR{ xL=%d\n", xL)
+ xa[0] = 0 ## arrify
+ za[0] = 0 ## arrify
+ xhix = bn__arr_nsplit(x, bixLUmul, xa, 0)
+ #printf("USQR s.\n")
+ #printf("bixLUmul=%d xhix=%s x=%s\n", bixLUmul, xhix, x)
+ #for( j=0 ; j<=xhix ; ++j) printf(" xa[%d] = %4s\n", j,xa[j])
+
+ ## for(j) for(k) with j<k: += x[j]*x[k]
+ for( j=0 ; j<=xhix ; ++j ) {
+ xaj = xa[j]
+ if( xaj ) {
+ c = 0
+ for( k=j+1 ; k<=xhix ; ++k ) {
+ c = bn__arrADD(za,j+k, xaj*xa[k] + c, bixLUmul)
+ }
+ while( c ) {
+ c = bn__arrADD(za,j+k, c, bixLUmul); ++k
+ }
+ }
+ }
+ ## Duplicate, what we have up to now...
+ c = 0
+ zhix = 2*xhix + 1 ## 2*(xhix+1) - 1
+ for( j=0 ; j<=zhix ; ++j ) {
+ c = bn__arrADD(za,j, za[j]+c, bixLUmul)
+ }
+ if( c ) { bix_Cry("c1", "bn_Usqr", c); return bix_Err(1) }
+ ## Add diagonal: squares, just once
+ c = 0
+ for( j=0 ; j<=xhix ; ++j ) {
+ k = j+j
+ c = bn__arrADD(za,k , xa[j]*xa[j]+c, bixLUmul)
+ c = bn__arrADD(za,k+1, c, bixLUmul)
+ }
+ if( c ) { bix_Cry("c2", "bn_Usqr", c); return bix_Err(1) }
+ #printf("USQR Y\n")
+ if( (zhix+1) in za ) { bix_Cry("zhix", "bn_Usqr", zhix); return bix_Err(1) }
+ ## Skip hi zeros
+ while( (zhix>=0) && !za[zhix] ) --zhix
+ c = bn__arr_concL(za,0,zhix-1, bixLUmul)
+ #printf("}USQR\n")
+ return bix_o1("bn_Usqr", za[zhix] c )
+}
+
+##-----------------------------------------------------------------------------
+## Compare
+
+function bn_Ucmp(x,y ,xL,yL) # numeric -1|0|1: x<y | x==y | x>y
+{
+ if( ! ("y" in bixDont) ) bix_Sanity()
+ if( ! ("i" in bixDont) ) bix_i2("bn_Ucmp",x,y)
+
+ xL = length(x)
+ yL = length(y)
+ if( xL < yL ) return -1
+ if( xL > yL ) return 1
+ if( ("" x) < ("" y) ) return -1
+ if( ("" x) > ("" y) ) return 1
+ return 0
+}
+
+##-----------------------------------------------------------------------------
+## Div and Mod
+##
+## Again, we use the most simple technique, like with paper & pencil.
+## Also, we compute div and mod at the same time.
+## We return both values in two components of an array.
+##
+## Trivial method:
+## Ignoring the various lengths for addition and subtraction,
+## we shift left the second operand by one decimal place,
+## until the operation becomes trivial.
+## Then we use Ucmp() and Usub() to find the next digit of the quotient.
+## This does not impose any requirements on the built-in arithmetic.
+##
+## What do we need from the built-in arithmetic?
+## Since we can multiply (bixLUmul * bixLUmul) --> 2*bixLUmul,
+## we want the converse div/mod operation also to be possible and exact.
+
+function bn__divtrunc(x,y ,z) ## --> floor(x/y)
+{
+ ## Both, x and y are non-negative built-in numbers, and y>0
+ ## We try to compensate for possible (?) inexact division and/or
+ ## int(...) results.
+
+ x = 0 + x ## numerify
+ y = 0 + y ## numerify
+ z = int( x / y )
+ if( y*z > x ) --z
+ if( y*z > x ) {
+ bix_Cry("x/y=z", "bn__divtrunc", x "/" y "=" z); return bix_Err(1)
+ }
+ if( (x - (y*z)) >= y ) {
+ bix_Cry("x/y=z(%)", "bn__divtrunc", x "/" y "=" z); return bix_Err(1)
+ }
+ return z
+}
+
+function bn__Uapxdivtrunc(x,y ,xL,yL,nL) ## --> floor(x/y)--
+{
+ ## We divide aproximately (from the most significant digits).
+ ## The (numeric) result may be too small, but never too large.
+ ## We may: x--, y++
+ ## We may not: x++, y--
+ ## We expect y to be smaller, but not by more than length 2*bixLUmul-1.
+ ## FFS: check I/O
+ xL = length(x)
+ yL = length(y)
+ if( xL < yL ) return 0
+ ## yL <= xL
+ if( xL <= (2*bixLUmul) ) { ## both are small
+ return bn__divtrunc(x,y) ## we can do it exactly
+ }
+ nL = xL - (2*bixLUmul) ## >0, so many low digits ignored
+ x = 0 + substr(x, 1, xL-nL) ## may x--
+ if( yL <= nL ) { ## y vanishes by truncating nL
+ return x ## so we use y=1: x/1 == x
+ }
+ y = substr(y, 1, yL-nL) ## y--: needs ++
+ if( length(y) > bixLUinc ) { ## FFS: can that happen?
+ x = bn__divtrunc(x, 2)
+ y = bn__divtrunc(y, 2)
+ }
+ return bn__divtrunc(x, y+1)
+}
+
+function bn_UdivmodDig(x,y,res ,xL,yL,dL,nL,quot,rem,rx,j,r)
+{ ## res[0] = x/y, res[1] = x%y
+ ## res[0] gets the div value, res[1] gets the mod value
+ ## "y" is numeric, length <= bixLUmul
+ if( ! ("y" in bixDont) ) bix_Sanity()
+ if( ! ("i" in bixDont) ) bix_iI("bn_UdivmodDig",x,y)
+ #printf("bn_UdivmodDig(%s,%s)\n",x,y)
+ ## FFS: check y==1
+
+ xL = length(x)
+ yL = length(y)
+ dL = 2*bixLUmul - yL ## so many from x into quot in each step
+ nL = xL - yL
+ nL -= nL % dL ## rounded down to multiple of dL
+
+ quot = ""
+ rx = xL+1 - nL ## fetch next dL-digit from x here
+ rem = 0 + substr(x, 1, xL-nL) ## numeric!
+ for(;;) {
+ j = int( rem / y )
+ r = rem % y
+ if( (j*y + r) != rem ) {
+ bix_Cry("x/%y=q r", "bn_UdivmodDig", rem "/%" y "=" j " " r)
+ return bix_Err(1)
+ }
+ rem = r
+ quot = bn_conc(quot, j, dL)
+ if( rx > xL ) break
+ j = substr(x, rx, dL); rx += dL
+ if( rem || (0+j) ) rem = 0 + (rem j)
+ }
+ res[0] = quot
+ res[1] = bn_U(rem)
+
+ #printf("bn_UdivmodDig(%s,%s) --> (%s,%s)\n",x,y, res[0],res[1])
+ bix_o01("bn_UdivmodDig", res)
+}
+
+function bn_Udivmod(x,y,res ,xL,yL,rem,rx,j,quot,qa,qx)
+{ ## res[0] = x/y, res[1] = x%y
+ ## res[0] gets the div value, res[1] gets the mod value
+ if( ! ("y" in bixDont) ) bix_Sanity()
+ if( ! ("i" in bixDont) ) bix_i2("bn_Udivmod",x,y)
+ #printf("bn_Udivmod(%s,%s)\n",x,y)
+ ## The typical (important) case is: x is much larger than y.
+ ## FFS: concatenating the result via an array appears to be more expensive.
+
+ xL = length(x)
+ yL = length(y)
+
+ if( (xL < yL) || (bn_Ucmp(x,y) < 0) ) { # x<y : trivial part
+ res[0] = "" ## quot 0
+ res[1] = x
+ }else if( yL <= bixLUmul ) {
+ bn_UdivmodDig(x, 0+y, res)
+ }else { ## xL>=yL, x>=y, quot>=1
+ ## FFS: completely native for small x & y?
+ rem = substr(x, 1, yL)
+ rx = yL+1 ## index of next digit to use from x
+ #qx = 0
+ quot = ""
+ for(;;) {
+ ## find j, next digit of quot, and reduce rem by (j*y)
+ if( bn_Ucmp(rem,y) < 0 ) {
+ j = 0 ## sure, exact
+ }else {
+ j = 0 + bn__Uapxdivtrunc(rem,y)
+ if( j>9 ) {
+ bix_Cry("large q-digit", "bn_Udivmod", rem "/" y "=" j)
+ return bix_Err(1)
+ }
+ if( j == 1 ) {
+ rem = bn_Usub(rem, y)
+ }else if( j ) {
+ rem = bn_Usub(rem, bn_UmulDig(y,j,0))
+ }
+ while( bn_Ucmp(rem,y) >= 0 ) {
+ rem = bn_Usub(rem, y); ++j
+ if( j>9 ) {
+ bix_Cry("large q-digit", "bn_Udivmod", j ":" rem "/" y)
+ return bix_Err(1)
+ }
+ }
+ }
+ ## Append j as next digit to quot
+ #if( qx || j ) qa[--qx] = j
+ if( quot || j ) quot = quot j
+ ## Check for more digits to fetch from x
+ if( rx > xL ) break
+ ## Fetch next digit from x and append it to rem
+ j = substr(x,rx,1); ++rx
+ if( rem || (0+j) ) rem = rem j
+ }
+ #res[0] = bn__arr_conc(qa, qx, -1)
+ res[0] = quot
+ res[1] = rem
+ }
+
+ #printf("bn_Udivmod(%s,%s) --> (%s,%s)\n",x,y, res[0],res[1])
+ bix_o01("bn_Udivmod", res)
+}
+
+function bn_Umod(x,y ,res)
+{
+ ## We delegate checking to "bn_Udivmod".
+ ## FFS: we should tell it, that we need the mod value, only
+ res[0] = 0 ## arrify (FFS: indicate mod only)
+ bn_Udivmod(x,y,res)
+ return res[1]
+}
+
+## FFS: operation "xdiv": only quotient wanted, remainder must be 0,
+## i.e. the division is an exact one (otherwise error).
+
+##-----------------------------------------------------------------------------
+## GCD = greatest common divisor
+##
+## The most important application of this is to normalize rational numbers.
+## Therefore gcd(0,0) will not be needed to be defined very sensibly.
+## Either we make it an error, or we make it 1.
+##
+## We use Euclids algorithm, of course.
+## The binary gcd algorithm is not very attractive, since our numbers are
+## represented in base 10, not in base 2, what makes halving expensive.
+##
+## On U-numbers (non-negative) we define gcd as usual:
+## gcd(0,0) 1 (or error)
+## gcd(x,y) gcd(y,x) [x < y]
+## gcd(x,0) x [x > 0]
+## gcd(x,y) gcd(y,x%y) [x >= y]
+##
+## For x>=1: gcd(x,1) = gcd(1,0) = 1
+## NB: GCD result is never zero.
+
+function bn_Ugcd(x,y ,xL,yL,t)
+{
+ if( ! ("y" in bixDont) ) bix_Sanity()
+ if( ! ("i" in bixDont) ) bix_i2("bn_Ugcd",x,y)
+ #printf("bn_Ugcd(%s,%s)\n",x,y)
+
+ xL = length(x)
+ yL = length(y)
+
+ if( !xL && !yL ) return bix_o1("bn_Ugcd", bn_U(1) )
+ if( !xL ) return bix_o1("bn_Ugcd", y) ## gcd(0,y) = y
+ if( !yL ) return bix_o1("bn_Ugcd", x) ## gcd(x,0) = x
+
+ t = bn_Ucmp(x,y)
+ if( t == 0 ) return bix_o1("bn_Ugcd", x) ## gcd(x,x) = x
+ if( t < 0 ) { ## x < y
+ t = x ; x = y ; y = t ## exchange x and y
+ t = xL; xL = yL ; yL = t ## exchange xL and yL
+ }
+ ## x > y ; xL >= yL FFS: native arith.
+ while( yL ) {
+ t = bn_Umod(x,y) ## x > y > t
+ ## (x,y) <-- (y,t) ...
+ x = y; xL = yL
+ y = t; yL = length(y) ## x > y
+ }
+ return bix_o1("bn_Ugcd", x)
+}
+
+##=============================================================================
+## Non-exported arithmetic functions, non-checked I/O
+
+function bn__Sign(x) # numeric value: -1 | 0 | 1
+{
+ if( length(x) <= bixLinc ) { ## use builtin arithmetic
+ if( (0+x) > 0 ) return 1
+ if( (0+x) < 0 ) return -1
+ return 0
+ }
+ if( substr(x,1,1) == "-" ) return -1
+ return 1
+}
+
+function bn__Neg(x) # --> -x
+{
+ if( length(x) <= bixLadd ) return bn_nR( 0-x )
+ if( substr(x,1,1) == "-" ) return substr(x,2)
+ return "-" x
+}
+
+function bn__Abs(x) # --> |x|
+{
+ if( length(x) <= bixLadd ) return bn_nR( ((0+x) < 0) ? 0-x : 0+x )
+ if( substr(x,1,1) == "-" ) return substr(x,2)
+ return x
+}
+
+function bn__Cmp(x,y ,xL,yL,xs,ys) # numeric -1|0|1: sign(x-y)
+{
+ xL = length(x)
+ yL = length(y)
+ if( (xL <= bixLinc) && (yL <= bixLinc) ) {
+ if( (0+x) < (0+y) ) return -1
+ if( (0+x) > (0+y) ) return 1
+ return 0
+ }
+ xs = bn__Sign(x)
+ ys = bn__Sign(y)
+ if( xs < ys ) return -1
+ if( xs > ys ) return 1
+ if( !xs ) return 0-ys
+ if( xL < yL ) return 0-xs
+ if( xL > yL ) return xs
+ ## same length: use string compare
+ x = "" x ; y = "" y
+ if( x < y ) return 0-xs
+ if( x > y ) return xs
+ return 0
+}
+
+function bn__Add(x,y ,xL,yL,xs,ys) # --> x+y
+{
+ xL = length(x)
+ yL = length(y)
+ if( (xL <= bixLadd) && (yL <= bixLadd) ) return bn_nR( x+y )
+ if( ! (xs = bn__Sign(x)) ) return y # 0+y = y
+ if( ! (ys = bn__Sign(y)) ) return x # x+0 = x
+ ## both are non-zero: strip off signs; assert x,y strings
+ if( xs < 0 ) { x = substr(x,2); --xL }else { x = "" x }
+ if( ys < 0 ) { y = substr(y,2); --yL }else { y = "" y }
+ if( xs == ys ) return bn_mksc(xs) bn_Uadd(x,y,0)
+ ## different sign, subtract: if x>y: xs*(x-y)
+ ## if x<y: -xs*(y-x)
+ if( xL < yL ) return bn_mksc(0-xs) bn_Usub(y,x,0)
+ if( xL > yL ) return bn_mksc( xs) bn_Usub(x,y,0)
+ ### equal length: use string compare
+ if( x < y ) return bn_mksc(0-xs) bn_Usub(y,x,0)
+ if( x > y ) return bn_mksc( xs) bn_Usub(x,y,0)
+ return bn_nR(0) # x + (-x) = 0
+}
+
+function bn__Mul2(x ,xL,xs) # --> 2*x
+{
+ xL = length(x)
+ if( xL <= bixLadd ) return bn_nR( x*2 )
+ if( ! (xs = bn__Sign(x)) ) return "" # 2*0 = 0
+ ## non-zero: strip off sign; assert x string
+ if( xs > 0 ) return bn_Umul2("" x, 0)
+ return bn_mksc(xs) bn_Umul2(substr(x,2), 0)
+}
+
+function bn__Sub(x,y ,xL,yL,xs,ys) # --> x-y
+{
+ xL = length(x)
+ yL = length(y)
+ if( (xL <= bixLadd) && (yL <= bixLadd) ) return bn_nR( x-y )
+ if( ! (ys = bn__Sign(y)) ) return x # x-0 = x
+ if( ! (xs = bn__Sign(x)) ) return bnNeg(y) # 0-y = -y
+ ## both are non-zero: strip off signs; assert x,y strings
+ if( xs < 0 ) { x = substr(x,2); --xL }else { x = "" x }
+ if( ys < 0 ) { y = substr(y,2); --yL }else { y = "" y }
+ if( xs != ys ) return bn_mksc(xs) bn_Uadd(x,y,0) # x-(-y) = x+y
+ ## same sign, subtract: if x>y: xs*(x-y)
+ ## if x<y: -xs*(y-x)
+ if( xL < yL ) return bn_mksc(0-xs) bn_Usub(y,x,0)
+ if( xL > yL ) return bn_mksc( xs) bn_Usub(x,y,0)
+ ## equal length: use string compare
+ if( x < y ) return bn_mksc(0-xs) bn_Usub(y,x,0)
+ if( x > y ) return bn_mksc( xs) bn_Usub(x,y,0)
+ return bn_nR( 0 ) # x - x = 0
+}
+
+function bn__Inc(x)
+{
+ if( length(x) <= bixLinc ) return bn_nR( x+1 )
+ return bnAdd(x,1)
+}
+
+function bn__Dec(x)
+{
+ if( length(x) <= bixLinc ) return bn_nR( x-1 )
+ return bnSub(x,1)
+}
+
+function bn__Mul(x,y ,xL,yL,xs,ys)
+{
+ xL = length(x)
+ yL = length(y)
+ if( (xL <= bixLmul) && (yL <= bixLmul) ) return bn_nR( 0+ x*y )
+ ## at least one operand is large
+ if( ! (ys = bn__Sign(y)) ) return bn_nR(0) # x*0 = 0
+ if( ! (xs = bn__Sign(x)) ) return bn_nR(0) # 0*y = 0
+ ## both are non-zero: strip off signs; assert x,y strings
+ if( xs < 0 ) { x = substr(x,2); --xL }else { x = "" x }
+ if( ys < 0 ) { y = substr(y,2); --yL }else { y = "" y }
+ return bn_mksc(0+ xs*ys) bn_Umul(x,y)
+}
+
+function bn__Sqr(x ,xL,xs)
+{
+ xL = length(x)
+ if( xL <= bixLmul ) return bn_nR( 0+ x*x )
+ ## the operand is large
+ xs = bn__Sign(x)
+ if( ! xs ) return bn_nR(0) # 0*0 = 0
+ ## the result is always non-negative
+ if( xs < 0 ) return bn_Usqr(substr(x,2))
+ return bn_Usqr("" x )
+}
+
+function bn__DivMod(x,y,res ,xs,ys) ## res[0] = floor(x/y), res[1] = x%y
+{
+ ## When y is negative, the division must be exact & res[1]==0
+ ## x /% 0 --> error
+ ## 0 /% y --> ( 0,0)
+ ## x /% 1 --> ( x,0) [important common case]
+ ## x /% -1 --> (-x,0)
+
+ if( y == "1" ) { res[0] = x ; res[1] = 0 ; return } ## x/%1 = (x,0)
+
+ ys = bn__Sign(y)
+ if( ys == 0 ) { bix_Cry("arg2", "bn__DivMod", y); return bix_Err(1) }
+ xs = bn__Sign(x)
+ if( xs == 0 ) { res[0] = x ; res[1] = 0 ; return } ## 0/%y = (0,0)
+
+ ## both are non-zero: strip off signs; assert x,y strings
+ if( xs < 0 ) { x = substr(x,2) }else { x = "" x }
+ if( ys < 0 ) { y = substr(y,2) }else { y = "" y }
+
+ if( y == "1" ) { ## x/%(ys) = (ys*x,0)
+ res[0] = bn_mksc(ys*xs) x ; res[1] = 0 ; return
+ }
+
+ bn_Udivmod(x,y,res)
+ res[0] = bn0(res[0])
+ res[1] = bn0(res[1])
+
+ ## If y originally was negative, the division must be exact (zero res[1])
+ if( (ys < 0) && bnNE0(res[1]) ) {
+ x = bn_mksc(xs) x
+ y = bn_mksc(ys) y
+ bix_Cry("rem", "bn__DivMod", (x " /% " y ": rem=" res[1]))
+ return bix_Err(0)
+ }
+ ## If x & y had equal sign, the result is not negative, and already ok.
+ ## If x & y had different signs, the result is negative, but the
+ ## unsigned operation truncated the quotient towards 0.
+ ## We must correct this...
+ if( xs != ys ) {
+ res[0] = bnNeg(res[0]) ## quot = -quot
+ if( bnNE0(res[1]) ) { ## implies ys>0
+ res[0] = bnDec(res[0]) ## quot-- towards -Inf (floor)
+ res[1] = bnSub(y, res[1]) ## rem = y-rem
+ }
+ }
+}
+
+function bn__Div(x,y ,res) ## --> floor(x/y)
+{
+ res[0] = 0 ## arrify
+ #FFS: y<0: -x/-y und Korrektur fuer floor
+ bn__DivMod(x,y,res)
+ return res[0]
+}
+
+function bn__Mod(x,y ,res) ## --> x mod y
+{
+ res[0] = 0 ## arrify
+ bn__DivMod(x,y,res)
+ return res[1]
+}
+
+function bn__Gcd(x,y ,z,xs,ys) ## --> gcd(x,y)
+{
+ ## We want y / gcd(x,y) to be non-negative (for rationals).
+ ## Hence, the result shall be negated for negative y.
+
+ ys = bn__Sign(y)
+ xs = bn__Sign(x)
+ if( !xs && !ys ) return bn_nR(1) ## cf. bn_Ugcd
+ if( !xs ) return y ## sign ok: y/y > 0
+ if( !ys ) return x ## sign ok: 0/x >= 0
+
+ ## both are non-zero: strip off signs; assert x,y strings
+ if( xs < 0 ) { x = substr(x,2) }else { x = "" x }
+ if( ys < 0 ) { y = substr(y,2) }else { y = "" y }
+
+ z = bn_Ugcd(x,y) ## non-zero ==> valid bigint!
+ if( ys < 0 ) z = "-" z
+ return z
+}
+
+function bn__Log10(x ,xL,i,y) ## --> log10(|x|)
+{
+ x = bnAbs(x)
+ if( bnLE0(x) ) return "-Inf"
+ ## x > 0
+ y = 0
+ xL = length(x)
+ ## From low (right) digits to high (left) digits collect the upper,
+ ## significant part of the matissa as a floating point value
+ ## with the decimal point left of the mantissa: y = 0.xxxx, 0 <= y < 1.
+ ## FFS: do it with more digits at once.
+ i = xL ## index of least significant (decimal) digit
+ if( i > 60 ) i = 60 ## collect at most 60 significant digits
+ for( ; i>0 ; --i ) {
+ y += substr(x,i,1) ## "prepend" digit
+ y /= 10 ## shift it down bhind the decimal point
+ }
+ return xL + (log(y) / log(10))
+}
+
+##=============================================================================
+## Export section, checking I/O
+## Many functions here are just wrappers for the previous section.
+
+function bn0(x) # convert "" into bignum zero
+{
+ return (("" x) == "") ? 0 : x
+}
+
+function bnIsOk(x) # --> whether is a valid bignum
+{
+ return bn_IsI(x)
+}
+
+function bnCheck(x) # --> whether is a valid bignum
+{
+ if( bnIsOk(x) ) return 1
+ ++bix_errs #FFS: andere var dafuer? bix_NaNs?
+ printf("\n*** bignum: not a number: \"%s\"\n", x)
+ return 0
+}
+
+function bnSign(x) # numeric value: -1 | 0 | 1
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I1("bnSign",x)
+ return bix_OC("bnSign", bn__Sign(x))
+}
+
+function bnNeg(x) # --> -x
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I1("bnNeg",x)
+ return bix_O1("bnNeg", bn__Neg(x))
+}
+
+function bnAbs(x) # --> |x|
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I1("bnAbs",x)
+ return bix_O1("bnAbs", bn__Abs(x))
+}
+
+function bnCmp(x,y) # numeric -1|0|1: sign(x-y)
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I2("bnCmp",x,y)
+ return bix_OC("bnCmp", bn__Cmp(x,y))
+}
+
+function bnLT(x,y) { return bnCmp(x,y) < 0 }
+function bnLE(x,y) { return bnCmp(x,y) <= 0 }
+function bnGE(x,y) { return bnCmp(x,y) >= 0 }
+function bnGT(x,y) { return bnCmp(x,y) > 0 }
+function bnEQ(x,y) { return bnCmp(x,y) == 0 }
+function bnNE(x,y) { return bnCmp(x,y) != 0 }
+
+function bnLT0(x) { return bnSign(x) < 0 }
+function bnLE0(x) { return bnSign(x) <= 0 }
+function bnGE0(x) { return bnSign(x) >= 0 }
+function bnGT0(x) { return bnSign(x) > 0 }
+function bnEQ0(x) { return bnSign(x) == 0 }
+function bnNE0(x) { return bnSign(x) != 0 }
+
+function bnAdd(x,y) # --> x+y
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I2("bnAdd",x,y)
+ return bix_O1("bnAdd", bn__Add(x,y) )
+}
+
+function bnSub(x,y) # --> x-y
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I2("bnSub",x,y)
+ return bix_O1("bnSub", bn__Sub(x,y) )
+}
+
+function bnInc(x) # --> x+1
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I1("bnInc",x)
+ return bix_O1("bnInc", bn__Inc(x) )
+}
+
+function bnDec(x) # --> x-1
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I1("bnDec",x)
+ return bix_O1("bnDec", bn__Dec(x) )
+}
+
+function bnMul(x,y) # --> x*y
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I2("bnMul",x,y)
+ return bix_O1("bnMul", bn__Mul(x,y) )
+}
+
+function bnMul2(x) # --> x*2
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I1("bnMul2",x)
+ return bix_O1("bnMul2", bn__Mul2(x) )
+}
+
+function bnSqr(x) # --> x*x
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I1("bnSqr",x)
+ return bix_O1("bnSqr", bn__Sqr(x) )
+}
+
+function bnDivMod(x,y,res) ## res[0] = floor(x/y), res[1] = x%y
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I2("bnDivMod",x,y)
+ bn__DivMod(x, y, res)
+ bix_O01("bnDivMod", res)
+}
+
+function bnDiv(x,y) # --> floor(x/y)
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I2("bnDiv",x,y)
+ return bix_O1("bnDiv", bn__Div(x,y) )
+}
+
+function bnMod(x,y) # --> x mod y
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I2("bnMod",x,y)
+ return bix_O1("bnMod", bn__Mod(x,y) )
+}
+
+function bnGcd(x,y) # --> gcd(x,y)
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I2("bnGcd",x,y)
+ return bix_O1("bnGcd", bn__Gcd(x,y) )
+}
+
+function bnLog10(x) # --> log10(|x|)
+{
+ bix_Sanity(); if( ! ("I" in bixDont) ) bix_I1("bnLog10",x)
+ return bn__Log10(x)
+}
+
+##-----------------------------------------------------------------------------
+## Printing reduced numbers
+##
+## When a number becomes excessively long, we may want to show (print)
+## only a "relevant" portion of it:
+## - some initial (high) digits
+## - some trailing (low) digits
+## - an indicator of how many digits are not shown in the middle
+## Example: 12345[99]54321
+## Example: 12345{99}54321
+##
+## For the middle part to reduce anything at all, it must have length
+## at least 4, since we add 2 additional delimiter chars.
+## FFS: should the reduced part also have some minimal length > 4 ?
+##
+## bn__opt_H 0 | (half length) length of remaining prefix and suffix part
+## bn__opt_L 0 | minimum length for a number to be "long" (shortable)
+## bn__didShort whether we did shorten some number
+
+function bnGetOptLongHalf()
+{
+ return 0 + bn__opt_H
+}
+
+function bnSetOptLongHalf(h)
+{
+ h += 0
+ if( h < 1 ) {
+ bn__opt_L = 0
+ bn__opt_H = 0
+ }else {
+ bn__opt_H = h
+ bn__opt_L = bn__opt_H*2 + 4 ##FFS
+ #bn__opt_L = bn__opt_H*2 + 9 ##FFS
+ }
+ #printf("{{H=%s}}", bn__opt_H)
+ return bn__opt_H
+}
+
+function bnSetOptLong(lo)
+{
+ lo += 0
+ if( lo < 1 ) {
+ bnSetOptLongHalf(h)
+ }else if( ! bnSetOptLongHalf( int((lo+1-4)/2) ) ) {
+ bnSetOptLongHalf(2) ##FFS
+ }
+ return bn__opt_H
+}
+
+function bnIsLong(x)
+{
+ ## NB: input validity explicitly only quiet!
+ return bn__opt_L && (length(x) >= bn__opt_L) && bn_IsI(x)
+}
+
+function bnShort(x ,xL,lo,mid,hi)
+{
+ #printf("{{SH(L=%d) %d}}", length(x), bn__opt_L)
+ if( ! bnIsLong(x) ) {
+ return x
+ }
+ bn__didShort = 1
+ xL = length(x)
+ ##FFS: the leading minus sign should be taken out of the lengths
+ lo = substr(x, 1, bn__opt_H)
+ hi = substr(x, 1+xL-bn__opt_H)
+ mid = sprintf("[%d]", xL-2*bn__opt_H)
+ return lo mid hi
+}
+
+function bnDidShort() ## --> "" | the brackets used in shortening
+{
+ return bn__didShort ? "[]" : ""
+}
+
+##=============================================================================
+## Test Section
+
+function bix_show_incdec_num(base,upto ,x,p,y,ok,pok,s)
+{
+ x = 1
+ pok = 0
+ for( p=1 ; p<=upto ; ++p ) {
+ y = x * base
+ if( (y / x) != base ) {
+ printf("%d^%d is not exact: %s\n", base, p, y)
+ break
+ }
+ x = y
+ printf("%d^%-3d = %13.6e can ++/-- exact: ", base, p, x)
+ ok = bix_incdec_ok(x)
+ s = substr(substr(x,1,1) substr(x,1,1) substr(x,2), 2)
+ printf("%d [str: %d]\n", ok, bix_incdec_ok(s))
+ if( ok && (pok == (p-1)) ) pok = p
+ }
+ printf("Last ok was %d^%d\n", base, pok)
+}
+
+function bix_show_incdec_str(upto ,s,len,ok,lenok)
+{
+ s = ""
+ lenok = 0
+ for( len=1 ; len<=upto ; ++len ) {
+ s = "9" s "" ## gawk-fix
+ printf("Decimal length %3d can ++/-- exact: ", len)
+ ok = bix_incdec_ok(s)
+ printf("%d\n", ok)
+ if( ok && (lenok == (len-1)) ) lenok = len
+ printf("str %s =?= num %d %.0f\n", s, s, s)
+ }
+ printf("Last ok was decimal length %d\n", lenok)
+}
+
+function bix_show_test( s,ll,ok) # --> whether ok
+{
+ bix_optv = 1
+ bix_errs = 0
+ #OFMT = "%.0f"
+ #CONVFMT = "%.0f"
+ bix_show_incdec_num( 2,55)
+ #exit(0)
+ #bix_show_incdec_num( 2,55)
+ #bix_show_incdec_num(10,20)
+ #bix_show_incdec_str( 20)
+ printf("Declen ok pos = %2d\n", bix_id_declen(0,50))
+ printf("Declen ok neg = %2d\n", bix_id_declen(1,50))
+ #exit(0)
+ s = "0013"
+ printf("Conv '%s' = '%s' (%d)\n", s, 0+s, s)
+ #
+ for( ll=2 ; ll<10 ; ++ll ) {
+ ok = bn_Setup(ll) # test small length
+ #
+ printf("Setup(%d) = %d: %d\n", ll, bixLinc, ok)
+ bix_show_setup()
+ if( !ok ) continue
+ bix_Test(1)
+ }
+ bix_optv = 0
+ printf("Error count = %d\n", bix_errs)
+ return bix_errs == 0
+}
+
+function bix_try_test()
+{
+ #if( "012" ~ /[0-9]/ ) print "rats"
+ #if( "012" ~ /[0-9]+/ ) print "fine"
+ #if( "012" ~ /^[0-9]$/ ) print "2*rats"
+ #if( "012" ~ /^[0-9]+$/ ) print "2*fine"
+ ## NB: awk patterns match partially, and must be anchored explicitly,
+ ## if wanted
+ #bix_vtst = 1
+ if( ! bixLinc ) bn_Setup()
+ #bix_show_setup()
+ #res[0] = 0
+ #bn_Udivmod("100000","33",res)
+ bix_show_test()
+}
+
+##=============================================================================
+
+function bnSignature()
+{
+ return "bignum signature: LEN={" bix_L_signature() "} " bix_C_signature()
+}
+
+function bignum_app_signature()
+{
+ if( bignum_signature_ ) return
+ bignum_signature_ = 1
+ g_Id[++g_IdCnt] = bnSignature()
+}
+
+##=============================================================================
+## Module initialization
+
+function bignum_init() # module constructor
+{
+ if( bignum_inited_ ) return # did it, already
+ bignum_inited_ = 1 # committed to do so
+ g_Id[++g_IdCnt] = "$Id: bignum.awk,v 1.34 2010/05/06 17:58:14 heiner Exp $" # note our identity
+
+ #bnSetup()
+ bnSetup(9)
+}
+
+##-----------------------------------------------------------------------------
+## End of module "bignum"
+##=============================================================================
+##=============================================================================
+## $Id: varLI.awk,v 1.11 2005/01/15 21:01:29 heiner Exp $
+## Athor: Heiner Marxen
+## Module: varLI -- Variables, linear integer
+## Prefix: liv
+##=============================================================================
+
+## Variables are denoted "V(%d)", i.e. they are numbered.
+## Variables for this module are always integral and non-negative !!
+##
+## Expressions built with variables are strings, containing a fixed prefix
+## bignum integer, and an optional, "+" separated list of linear variable
+## terms. Each variable term consists of an optional integral bignum
+## factor connected by "*" with a single Variable.
+##
+## Example: 1+7*V(1) == 7V(1) + 1
+## 0+V(1) == V(1)
+## -2+-1*V(1)+V(2) == V(2) - V(1) - 2
+## Variables with a factor 0 are omitted.
+## A factor 1 is omitted.
+##
+## Note, that non of the special chars "+", "*" and "V()" can occur
+## in a legal bignum.
+##
+## We only use normalized expressions (NEs). (FFS: "" for "0")
+## In an NE the constant part is always first, and never omitted.
+## The variable parts are sorted by increasing variables (dictionary order).
+## No variable appears twice.
+
+## Internally we represent the linear form in a single array, which maps
+## - "" to the constant part
+## - each contained variable (name!) to its factor
+
+## To split an expression into its components, we could use "split()".
+## Since we have had portability problems with it, we do not do so,
+## i.e. we avoid split(). We do not expect a significant speed penalty.
+
+##-----------------------------------------------------------------------------
+## Export interface:
+##
+## livVname(vx) --> var name of var index
+## livIsOk(x) --> whether x is a valid operand
+## livCheck(x) returns for valid operand, only
+## otherwise complains and exits program
+##
+## livIsConst(x) --> whether guaranteed to be constant
+## livEQ0(x) --> whether guaranteed to be == 0
+## livGT0(x) --> whether guaranteed to be > 0
+## livEQ(x,y) --> whether guaranteed identical
+## livGEconst(x,y) --> x>=y guaranteed, and y is const
+##
+## livNeg(x) --> -x
+## livAdd(x,y) --> x + y
+## livSub(x,y) --> x - y
+## livMul(x,y) --> x * y
+##
+## livSubst(vnval,x) --> x with var V mapped to vnval[V]
+## livSum0lt(vn,lim,x) --> sum(vn=0,vn<lim,++vn: x)
+
+##=============================================================================
+## Split Expression into Array
+
+function liv_F(fv ,p) ## --> factor part
+{
+ if( p = index(fv,"*") ) {
+ return substr(fv,1,p-1) ## factor part
+ }
+ return 1 ## implicit factor 1
+}
+
+function liv_V(fv ,p) ## --> variable part (name)
+{
+ if( p = index(fv,"*") ) {
+ return substr(fv,1+p) ## variable part
+ }
+ return fv ## must be variable, completely
+}
+
+function liv_A_addFV(arr,f,v) ## --> whether new
+{ ## arr += f*v in internal array format
+ if( bnEQ0(f) ) return 0 ## nothing to be added
+ if( v in arr ) {
+ arr[v] = bnAdd(arr[v], f)
+ return 0 ## no, that entry was present, already
+ }
+ arr[v] = f
+ return 1 ## yes, new entry created
+}
+
+function liv_fv_toA(fv,arr) ## --> whether new
+{
+ return liv_A_addFV(arr, liv_F(fv), liv_V(fv))
+}
+
+function liv_toA(x,arr ,p,n) ## --> #(vars) in there
+{
+ if( p = index(x,"+") ) {
+ arr[""] = substr(x,1,p-1) ## constant part
+ n = 0
+ x = substr(x,p+1)
+ while( p = index(x,"+") ) {
+ n += liv_fv_toA(substr(x,1,p-1),arr)
+ x = substr(x,p+1)
+ }
+ n += liv_fv_toA(x,arr)
+ return n
+ }
+ arr[""] = x
+ return 0 ## no var in there
+}
+
+##-----------------------------------------------------------------------------
+## Construct Expression from Array
+
+function liv_mk_fv(f,v)
+{
+ return bnEQ(f,1) ? v : (f "*" v)
+}
+
+function liv_mk(arr ,v,n,var,z,i,j,b)
+{ ## --> string from array
+ n = 0
+ for( v in arr ) {
+ if( v == "" ) continue ## skip non-var part
+ if( bnEQ0(arr[v]) ) continue ## ignore factor 0
+ var[++n] = v
+ }
+ z = arr[""] ## fixed non-var part
+
+ ## simple selection sort...
+ for( i=1 ; i <= n ; ++i ) {
+ b = i ## smallest, so far
+ for( j=i+1 ; j <= n ; ++j ) {
+ if( var[j] < var[b] ) b = j ## even smaller
+ }
+ z = z "+" liv_mk_fv(arr[var[b]], var[b])
+ ## We do not really sort, but rather skip slot var[i].
+ ## If that did not contain the just used b, we must
+ ## save the not use [i] in the still to use [b].
+ if( b != i ) var[b] = var[i]
+ }
+ return z
+}
+
+##=============================================================================
+## Checking Validity
+
+function liv__isV(x) ## whether is a var name
+{
+ return x ~ /^V[(][0-9]+[)]$/
+}
+
+function liv__is_fv(x ,p) ## whether is a "f*v" | "v"
+{
+ if( p = index(x,"*") ) {
+ if( ! bnIsOk(substr(x,1,p-1)) ) return 0
+ x = substr(x,1+p)
+ }
+ return liv__isV(x)
+}
+
+function liv__is_fvl(x ,p) ## whether is a f*v list
+{
+ while( p = index(x,"+") ) {
+ if( ! liv__is_fv(substr(x,1,p-1)) ) return 0
+ x = substr(x,1+p)
+ }
+ return liv__is_fv(x)
+}
+
+##=============================================================================
+## Exported Functions
+
+function livVname(vx) ## --> var name of var index
+{
+ return "V(" vx ")"
+}
+
+function livIsOk(x ,p) ## --> whether x is a valid operand
+{
+ if( p = index(x,"+") ) {
+ if( ! liv__is_fvl(substr(x,1+p)) ) return 0 ## variable part
+ x = substr(x,1,p-1) ## constant part
+ }
+ return bnIsOk(x)
+}
+
+function livCheck(x,fname) ## returns for valid operand, only
+{
+ if( ! livIsOk(x) ) {
+ printf("\n*** varLI: %s: illegal value \"%s\"\n", fname, x)
+ do_exit(1)
+ return 0 ## sorry
+ }
+ return 1 ## ok
+}
+
+##-----------------------------------------------------------------------------
+## The following tests must also work on any illegal operand.
+## They just have not the asked for property.
+
+function livIsConst(x) ## --> whether guaranteed to be constant
+{
+ ## This effectively reduces to: it is a valid bignum...
+ return bnIsOk(x)
+}
+
+function livEQ0(x) ## --> whether guaranteed to be == 0
+{
+ ## NB: we allow normalized representations, only
+ return x == "0"
+}
+
+function livGT0(x ,arr,v) ## --> whether guaranteed to be > 0
+{
+ if( ! livIsOk(x) ) return 0 ## sorry, not even a valid input
+ arr[""] = 0 ## arrify & assign 0
+ liv_toA(x, arr)
+ for( v in arr ) {
+ if( v == "" ) {
+ if( ! bnGT0(arr[v]) ) return 0 ## sorry
+ }else {
+ if( ! bnGE0(arr[v]) ) return 0 ## sorry
+ }
+ }
+ return 1 ## yes
+}
+
+function livEQ(x,y) ## --> whether guaranteed identical
+{
+ ## NB: we allow normalized representations, only
+ ## FFS: should non-operands compare different, like NaN's ?
+ #return (("" x) == ("" y)) && livIsOk(x)
+ return ("" x) == ("" y)
+}
+
+function livGEconst(x,y ,arr,v) ## --> x>=y guaranteed, and y is const
+{
+ if( ! livIsOk(x) ) return 0 ## sorry, not even a valid input
+ if( ! bnIsOk(y) ) return 0
+ arr[""] = 0 ## arrify & assign 0
+ liv_toA(x, arr)
+ for( v in arr ) {
+ if( v == "" ) {
+ if( ! bnGE(arr[v], y) ) return 0 ## sorry
+ }else {
+ if( ! bnGE0(arr[v]) ) return 0 ## sorry
+ }
+ }
+ return 1 ## yes
+}
+
+##-----------------------------------------------------------------------------
+## The following arithmetic operations only accept valid operands
+
+function livNeg(x ,arr,v) ## --> -x
+{
+ livCheck(x,"livNeg")
+ arr[""] = 0 ## arrify & assign 0
+ liv_toA(x, arr)
+ for( v in arr ) {
+ arr[v] = bnNeg(arr[v])
+ }
+ return liv_mk(arr)
+}
+
+function livAdd(x,y ,xa,ya,v) ## --> x + y
+{
+ livCheck(x,"livAdd"); livCheck(y,"livAdd")
+ xa[""] = 0 ## arrify & assign 0
+ ya[""] = 0 ## arrify & assign 0
+ liv_toA(x, xa)
+ liv_toA(y, ya)
+ for( v in ya ) liv_A_addFV(xa, ya[v], v)
+ return liv_mk(xa)
+}
+
+function livSub(x,y) ## --> x - y
+{
+ #livCheck(x,"livSub"); livCheck(y,"livSub")
+ return livAdd(x, livNeg(y))
+}
+
+function livMul(x,y ,xa,ya,v,xnv,ynv) ## --> x * y
+{
+ livCheck(x,"livMul"); livCheck(y,"livMul")
+ xa[""] = 0 ## arrify & assign 0
+ ya[""] = 0 ## arrify & assign 0
+ xnv = liv_toA(x, xa)
+ ynv = liv_toA(y, ya)
+ if( ! xnv ) {
+ for( v in ya ) ya[v] = bnMul(ya[v], x)
+ return liv_mk(ya)
+ }else if( ! ynv ) {
+ for( v in xa ) xa[v] = bnMul(xa[v], y)
+ return liv_mk(xa)
+ }else {
+ printf("\n*** varLI: Cannot (%s)*(%s)\n", x, y)
+ do_exit(1)
+ return "NaN"
+ }
+}
+
+function livSubst(vnval,x ,xa,ya,za,v,w)
+{ ## --> x with (V --> vnval[V])
+ ## A linear combination of linear combinations is a linear combination,
+ ## again, so we should be able to do this always!
+ livCheck(x,"livSubst")
+ #pr("livSubst: " x "\n")
+ #for(v in vnval) pr("with " v " --> " vnval[v] "\n")
+
+ xa[""] = 0 ## arrify & assign 0
+ liv_toA(x, xa)
+ ## We build the new expression in a separate array.
+ ## That avoids confusion with substituted vars.
+ ## FFS: we do not re-subst (repeat)
+ za[""] = 0 ## arrify & assign 0
+ for( v in xa ) {
+ #pr("livSubst: (" xa[v] ")*(" v "):\n")
+ if( v == "" ) { ## must never be substituted (FFS)
+ liv_A_addFV(za, xa[v], v)
+ }else if( v in vnval ) { ## known var, to be substituted
+ ## We break the value into parts, multiply the parts with
+ ## the current factor, and collect.
+ #pr("livSubst: use " vnval[v] ":\n")
+ livCheck(vnval[v],"livSubst(v)")
+ ya[""] = 0 ## arrify (portability work around)
+ for( w in ya ) delete ya[w] ## delete old garbage
+ ya[""] = 0 ## arrify & assign 0
+ liv_toA(vnval[v], ya)
+ for( w in ya ) {
+ #pr("livSubst: add in (" xa[v] "*" ya[w] ")*" w ":\n")
+ liv_A_addFV(za, bnMul(xa[v],ya[w]), w)
+ }
+ }else { ## unknown var: leave alone
+ liv_A_addFV(za, xa[v], v)
+ }
+ }
+ return liv_mk(za)
+}
+
+function livSum0lt(vn,lim,x ,v,xa) ## --> sum(vn=0,vn<lim: x)
+{
+ livCheck(lim,"livSum0lt"); livCheck(x,"livSum0lt")
+ if( ! liv__isV(vn) ) {
+ printf("\n*** varLI: livSum0lt: illegal var \"%s\"\n", vn)
+ do_exit(1)
+ }
+ ## If the limit is not const, we had to express a quadratic term...
+ ## FFS: Strictly, a constant body would work, but we do not need that.
+ if( ! livIsConst(lim) ) {
+ printf("\n*** varLI: livSum0lt: not const: \"%s\"\n", lim)
+ do_exit(1)
+ }
+
+ ## With lim==1 we have exactly one summand with value 0.
+ ## Hence, lim is also the number of terms in the sum.
+ if( bnLE0(lim) ) {
+ return 0 ## empty sum --> 0
+ }
+
+ ## All terms in x, which are not just for the var "vn", can be
+ ## treated as constant with respect to the summation,
+ ## and are just multiplied by the number of sum-terms, i.e. "lim".
+ ## Of the vn-term f*vn we can extract the f before the sum,
+ ## leaving us with the elementary quadratic term:
+ ## sum(vn=0, vn<lim: vn)
+ ## == (lim-1) * lim/2
+ ## == ((lim-1)*lim) / 2
+ ## which we can compute exactly, since we have given lim as integral
+ ## constant...
+
+ xa[""] = 0 ## arrify & assign 0
+ liv_toA(x, xa)
+ ## We build the new expression in place.
+ ## The only cross term effect is with the vn-term, which adds
+ ## to the "" term. But that is done last...
+ for( v in xa ) {
+ if( v == vn ) continue
+ xa[v] = bnMul(xa[v], lim) ## *= lim
+ }
+
+ if( vn in xa ) {
+ #v = bnDiv(bnMul(lim, bnDec(lim)), 2) ## (lim * (lim-1)) / 2
+ ## lim * (lim-1) == lim^2 - lim == (lim-1)^2 + (lim-1)
+ v = bnDiv( bnSub( bnSqr(lim), lim ), 2)
+ xa[""] = bnAdd(xa[""], bnMul(v, xa[vn]))
+ delete xa[vn] ## just replaced
+ }
+ return liv_mk(xa)
+}
+
+##=============================================================================
+## Module initialization
+
+function varLI_init()
+{
+ if( varLI_inited_ ) return ## already done
+ varLI_inited_ = 1 ## committed to do so
+ ## register us globally:
+ g_Id[++g_IdCnt] = "$Id: varLI.awk,v 1.11 2005/01/15 21:01:29 heiner Exp $"
+
+ bignum_init()
+ ## no own initialization, here
+}
+
+##-----------------------------------------------------------------------------
+## End of module "varLI"
+##=============================================================================
+##=============================================================================
+## $Id: htSupp.awk,v 1.14 2010/07/06 19:48:32 heiner Exp $
+## Author: Heiner Marxen
+## Module: htSupp -- HTML Quoting & Attributes
+
+##=============================================================================
+## - html whether producing HTML (otherwise plain ASCII text)
+##
+## - hQ[x] how to quote character "x" within HTML
+## - hA[x] HTML tag for attribute character "x"
+## - hCC[x] named color for index x (1 <= x <= hccCnt)
+## - hccCnt #(named colors in cycle)
+
+function ht_init() # to be called from BEGIN
+{
+ hQ["&"] = "&"
+ hQ["\""] = """
+ hQ["<"] = "<"
+ hQ[">"] = ">"
+ # #xxxxxx numeric color font color=#xxxxxx
+ # {yyy} named color
+ hA["i"] = "I" # italic
+ hA["b"] = "B" # bold
+ hA["B"] = "BIG" # big
+ hA["s"] = "SMALL" # small
+ hA["u"] = "U" # underlined
+ hA["^"] = "SUP" # superscript
+ hA["v"] = "SUB" # subscript
+ # cyclic color for states
+ hccCnt = 0
+ hCC[++hccCnt] = "blue"
+ hCC[++hccCnt] = "red"
+ hCC[++hccCnt] = "fuchsia"
+ hCC[++hccCnt] = "green"
+ hCC[++hccCnt] = "maroon"
+ hCC[++hccCnt] = "navy"
+ hCC[++hccCnt] = "purple"
+ hCC[++hccCnt] = "olive"
+ hCC[++hccCnt] = "black"
+ #
+ hCAhead = "b" # color attribute for tape head
+ hCAhalt = "b" # color attribute for halt state
+ #
+ htVtSelf = "Simulation is done "
+ htVtOthr = "The same TM "
+ #
+ htPref = "sim"
+}
+
+function htQ(s ,l,i,t,r) # --> "s" quoted for HTML
+{ # FFS: speed up
+ if( ! html ) return s # no quoting needed
+ r = 0
+ for( t in hQ ) {
+ if( r = index(s, t) ) break
+ }
+ if( r == 0 ) return s
+ r = ""
+ l = length(s)
+ for( i=1 ; i<=l ; ++i ) {
+ t = substr(s,i,1)
+ r = r ((t in hQ) ? hQ[t] : t)
+ }
+ return r
+}
+
+function htA(attrs,be ,z,a,r,s,c,i) # --> begin or end of HTML "attrs"
+{
+ #printf("[htA(%s,%s)]", attrs,be)
+ if( ! html ) return ""
+ if( length(attrs) <= 0 ) return ""
+ a = substr(attrs, 1, 1) # first char is attr itself
+ r = substr(attrs, 2)
+ s = "" ; c = ""
+ if( a in hA ) { # known (simple) attribute
+ s = hA[a]
+ }else if( a == "#" ) { # followed by 6 hex digits
+ c = "#" substr(r,1,6)
+ r = substr(r,7)
+ }if( a == "{" ) {
+ if( i = index(r, "}") ) {
+ c = "\"" substr(r,1,i-1) "\""
+ r = substr(r,i+1)
+ }
+ }
+ if( c ) s = "FONT" ((be=="/") ? "" : (" color=" c))
+ if( s ) z = "<" be s ">" ;else z = ""
+ r = htA(r, be) # recursively do others (rest)
+ return (be != "/") ? z r : r z # close others (inner) first
+}
+
+function htAaug(txt, attrs) # txt is quoted, already
+{
+ return htA(attrs, "") txt htA(attrs, "/")
+}
+
+function htCCaug(txt, num) # txt is quoted, already
+{
+ if( hccCnt < 1 ) return txt
+ if( num < 0 ) num = 0-num
+ num -= 1
+ if( num < 0 ) num += hccCnt
+ num %= hccCnt
+ num += 1
+ return htAaug(txt, "{" hCC[num] "}")
+}
+
+## FFS: color cycle for state names
+##-----------------------------------------------------------------------------
+## Printing Support
+
+function pr(s)
+{
+ printf("%s", htQ(s))
+}
+
+function prA(txt, attrs) # txt is not yet quoted
+{
+ printf("%s", htAaug(htQ(txt), attrs))
+}
+
+function prB(cnt)
+{ # we know, that blanks are safe
+ if( cnt > 0 ) printf("%*s", cnt, " ")
+}
+
+function nl()
+{ # we know, that newlines are safe
+ printf("\n") # while we are inside <PRE>
+}
+
+##-----------------------------------------------------------------------------
+## Comment line collection
+##
+## - htCmt[1] First comment line (to be printed)
+## - htNcmt #(collected htCmt), last used index
+
+function ht_CmtAdd(txt)
+{
+ if( txt ) {
+ htCmt[++htNcmt] = txt
+ }
+}
+
+function ht_CmtPr( i)
+{
+ if( htNcmt > 0 ) {
+ for( i=1 ; i<=htNcmt ; ++i ) {
+ pr("Comment: " htCmt[i]); nl()
+ #FFS each comment as paragraph, outside PRE
+ }
+ nl()
+ }
+}
+
+##-----------------------------------------------------------------------------
+## Print registered Ids & document saved input
+##
+## - htIL[] remember input (to document it in output)
+## - htILcnt
+
+function ht_Iline(line)
+{
+ htIL[ ++htILcnt ] = line
+}
+
+function ht_ByPr(preembed,inl,tnl ,pref,i)
+{
+ if( g_IdCnt || htILcnt ) {
+ if( inl ) nl()
+ if( ! html ) preembed = 0
+ if( preembed ) printf("<SMALL><PRE>")
+
+ if( htILcnt ) {
+ pr("Input to awk program:"); nl()
+ for( i=1 ; i<=htILcnt ; ++i ) {
+ prB(4); pr(htIL[i]); nl()
+ }
+ if( g_IdCnt ) {
+ if( html ) {
+ printf("<HR>")
+ }else {
+ nl()
+ }
+ }
+ }
+
+ if( g_IdCnt ) {
+ pref = "Constructed by: "
+ for( i=1 ; i<=g_IdCnt ; ++i ) {
+ pr(pref g_Id[i]); nl()
+ pref = sprintf("%*s", length(pref), "")
+ }
+ }
+
+ if( preembed ) printf("</PRE></SMALL>")
+ if( tnl ) nl()
+ }
+}
+
+##-----------------------------------------------------------------------------
+## Version registering
+## - htVers[1] first version (name)
+## - htVtxt[nam] maps version name to version text
+## - htVcnt #(versions)
+## - htPref URL prefix for page versions
+## - htViam my own version (name)
+## - htVtSelf
+## - htVtOthr
+
+function ht_AddVers(name, txt)
+{
+ ++htVcnt
+ htVers[htVcnt] = name
+ htVtxt[name ] = txt
+}
+
+function ht_VersPr(refonly ,i,c) # --> printed lines
+{
+ # refonly: leave out myself; also add <BR> as line breaks.
+ c = 0
+ for( i=1 ; i<= htVcnt ; ++i ) {
+ if( htVers[i] == htViam ) {
+ if( refonly ) continue
+ pr(sprintf("%s%s.\n", htVtSelf, htVtxt[htViam]))
+ ++c
+ }else if( html ) {
+ pr(htVtOthr)
+ printf("<A HREF=\"%s%s.html\">", htPref, htVers[i])
+ pr(htVtxt[htVers[i]])
+ printf("</a>.%s\n", refonly ? "<BR>" : "")
+ ++c
+ }
+ }
+ return c
+}
+
+##=============================================================================
+## Page embedding
+##
+## - htTit title
+## - htDate construction date
+## - htEdate date for end of computation (marker to substitute)
+## - HM whether embedding for Heiner's BB page
+
+function ht_BEG()
+{
+ if( html ) {
+ printf("<!DOCTYPE HTML PUBLIC \"-//W3C//DTD HTML 3.2//EN\">\n")
+ printf("<HTML><HEAD><TITLE>\n")
+ pr(htTit) ; nl()
+ printf("</TITLE></HEAD><BODY><H2>\n")
+ pr(htTit) ; nl()
+ printf("</H2><PRE>\n")
+ }else {
+ if( htTit != "" ) printf("%s\n\n", htTit)
+ }
+ #ht_ByPr(0,0,1)
+ ht_CmtPr()
+ #ht_VersPr(0)
+}
+
+function ht_ChkProfile( hasprof)
+{
+ ## Without aaP function profiling this function will not be called,
+ ## but "mawk" still wants to compile its body, and produces errors
+ ## for any unknown function like aaP_has_pr().
+ hasprof = 0
+ #{@}aaP:# if( aaP_has_pr() ) hasprof = 1
+ if( hasprof ) {
+ #{@}aaP:# aaP_freeze()
+ if( html ) {
+ printf("<PRE>")
+ nl()
+ }else {
+ nl()
+ printf("Function profiling ...")
+ nl()
+ }
+ #{@}aaP:# aaP_pr()
+ if( html ) {
+ printf("</PRE>")
+ }else {
+ nl()
+ }
+ #{@}aaP:# aaP_unfreeze()
+ }
+}
+
+function ht_END()
+{
+ if( html ) {
+ printf("</PRE><HR>\n")
+ if( ht_VersPr(1) ) printf("<HR>")
+ if( HM ) {
+ printf("To the <A HREF=\"bbsimtab.html\">\n")
+ printf("BB simulations page</A> of Heiner Marxen.<BR>\n")
+ printf("To the <A HREF=\"index.html\">busy beaver page</A>\n")
+ printf("of Heiner Marxen.<BR>\n")
+ printf("To the <A HREF=\"http://www.drb.insel.de/~heiner/\">\n")
+ printf("home page</A> of Heiner Marxen.<HR>\n")
+ }
+ ht_ByPr(1,0,0)
+ if( htDate ) {
+ printf("<SMALL>Start: %s</SMALL><BR>", htDate)
+ }
+ if( htEdate ) {
+ printf("<SMALL>Ready: %s</SMALL><BR>", htEdate)
+ }
+ #{@}aaP:# ht_ChkProfile()
+ printf("</BODY></HTML>\n")
+ }else {
+ ht_ByPr(0,1,0)
+ if( htDate || htEdate ) {
+ nl()
+ if( htDate ) {
+ printf("Start: %s\n", htDate)
+ }
+ if( htEdate ) {
+ printf("Ready: %s\n", htEdate)
+ }
+ }
+ #{@}aaP:# ht_ChkProfile()
+ }
+}
+
+##=============================================================================
+## Module initialization
+
+function htSupp_init()
+{
+ if( htSupp_inited_ ) return # already done
+ htSupp_inited_ = 1 # committed to do so
+ # register us globally:
+ g_Id[++g_IdCnt] = "$Id: htSupp.awk,v 1.14 2010/07/06 19:48:32 heiner Exp $"
+
+ ht_init()
+}
+
+##-----------------------------------------------------------------------------
+## End of module "htSupp"
+##=============================================================================
+#==============================================================================
+# $Id: mmSim.awk,v 1.34 2005/01/09 22:23:28 heiner Exp $
+# Athor: Heiner Marxen
+# Module: mmSim -- Multiple Machines Simulated
+#==============================================================================
+# Design principles for this software in decreasing order of importance:
+# - correctness
+# - clarity & maintainability of source
+# - portability
+# - (algorithmic) speed
+
+#------------------------------------------------------------------------------
+# In order to completely run a record machine, performing millions
+# of steps, we must not do it directly. Not even 1-macro machines suffice.
+# Hence we build higher machines on top of the basic machine.
+# The most important way to do this, are K-macro machines:
+# - The tape is cut into cells of size K. Each such cell contains
+# K basic symbols, which form one higher symbol.
+# - Repetitions of the same (higher) symbol are explicitly counted (exponents).
+# Such a tape is called a "macro tape".
+# - A macro state is built from a basic state, and a direction, into
+# which the macro head is looking. Hence, we have up to twice as
+# many macro-states as basic states.
+# The transitions of the higher machine are calculated as needed, by running
+# the basic machine on a single cell (K basic symbols), until it falls off.
+# Note, that not all basic configurations correspond to a macro configuration.
+# This is a feature, not a bug.
+
+# Of course, in order to construct a macro machine the lower machine does
+# not need to be a basic machine, it may itself be built upon another machine.
+
+# There are other ways to build higher machines, e.g.:
+# - Back-Context machine: the tape symbol behind the head is coded into
+# the state. Example: instead of
+# a b c <S x y z we now code
+# a b c <S(x) y z where "<S(x)" is the new state
+# This can help to encapsulate multiple crawling forth and back.
+
+# We assign IDs to the different machines built on top of the basic machine.
+# The basic machine has mID=1, the other machines get higher (numeric) mIDs.
+# - mIDtop currently active top machine
+# - mIDbas = 1, most basic machine
+
+# There are several mappings involved for higher machines:
+# - State mapping to/from the lower machine(s).
+# - Symbol mapping to/from basic symbol strings
+# - Transition mapping: those of the higher machine use some of the lower.
+# - Step counts and tape positions
+
+#------------------------------------------------------------------------------
+# Legend:
+# FFS for further study
+# TR transition
+# TRX transition index
+# TT transition table
+# LT limited tape
+# mID machine identification
+# ITR incomplete transition
+# Yspc symbol space
+# Sspc state space
+# RST run-state
+# GLT globally registered LT (neg symbol)
+#==============================================================================
+# Directions:
+# The direction of the tape head (into which it "looks") and relative tape
+# head movements internally are coded numerically:
+# -1 left ("<")
+# 0 none ("=")
+# 1 right (">")
+
+# Basic machines:
+# We always start with a 4-tuple or 5-tupel (4T/5T) machine.
+# We use a transition coding which fits for both:
+# - Head movement may be -1, 1 or 0 (left/right/stay).
+# - Non-printing 4T-transitions just re-print the input symbol.
+# - There is a global flag to allow/deny execution of a halting transition.
+
+# Recognition of 4T-machines:
+# Building a 5T on top of a 4T machine can help, if there is a transition,
+# which does not move and also does not go into a halt state.
+# The 4T machine is run on limited tapes of length 1.
+# Since cell count = 1, state and symbol mapping is trivially 1:1.
+
+#==============================================================================
+# Generally "steps" are counted from 0. A config/tape is shown with
+# the #steps that produced it, i.e. the initial config has steps=0.
+
+# Big numbers:
+# Running record machines easily can produce numbers which will not fit
+# into 32 bit, or even not into the typical 53 bit of a IEEE double mantissa.
+# Hence, we should use arbitrary length signed integer arithmetic.
+# Big numbers occur there:
+# - tape exponents
+# - tape position, tape head movement, tape scanning range
+# - TR-execution statistics
+# - step counts
+
+#------------------------------------------------------------------------------
+# For record machines we will have to compute with natural numbers much
+# larger than the numbers built in into "awk".
+# - native numbers: we expect "awk" to be able to exactly compute up to 10^9.
+# - big numbers: we use special functions to calculate with larger numbers
+# especially for steps and exponents. We limit ourselves to whole numbers.
+# - var numbers: when we start to prove "arithmetic steps", we will have
+# expressions containing polynomials with multiple variables in them,
+# where we otherwise would have pure numbers. Like big numbers,
+# we use awk-strings to code such expressions, and use functions to perform
+# arithmetic on them.
+
+# But even where we do not use big numbers (e.g. index numbers for states)
+# we try to detect overflow, by explicitly limiting all counters. FFS.
+
+#==============================================================================
+# General Machine mapping:
+# Building a higher machine on top of one or more lower machines, can include:
+# - building a new state space (Sspc)
+# - building a new symbol space (Yspc)
+# - computing transitions by running lower machines on a LT
+# Since not every configuration of the basic machine (especially a halting one)
+# corresponds to a configuration of the higher machine (head inside symbol),
+# we also need "run-states", incomplete transitions and augmented symbols.
+
+#------------------------------------------------------------------------------
+# States and State spaces:
+# A state is in one of the following classes:
+# - normal (corresponds to a basic state with transitions defined)
+# - halting (corresponds to a basic state without transitions defined)
+# - funny (result of an ITR)
+# One of the basic states is the initial state.
+# Higher machines can have states which carry more information.
+# They are constructed by augmenting a lower state.
+# Thus we construct state spaces by:
+# - initial: predefined basic states, normal | halting
+# - adding orientation (direction)
+# - adding backward symbol (must have direction, already)
+# State spaces are numbered (from 1)
+# - ssCnt #(state spaces) = last created Sspc
+# - ssMethod[ss] "ini" | "+dir" | "+bck"
+# - ssSubSspc[ss] 0 | state space below
+# - ssBckYspc[ss] 0 | symbol space (mID) of bck-symbols
+# The orientation attribute is inherited from lower state spaces:
+# - ssHasDir[ss] 0|1: whether it has Direction
+# This section does not contain any big/var numbers.
+
+function ss_new(meth, subss, bckys, hasdir ,ss) # --> "" | new ss
+{
+ if( subss == "" ) return ""
+ ss = ++ssCnt
+ ssMethod[ ss] = meth
+ ssSubSspc[ss] = subss
+ ssBckYspc[ss] = bckys
+ ssHasDir[ ss] = hasdir
+ return ss
+}
+
+function ss_find(meth, subss, subys ,ss) # --> 0 | ss
+{
+ for( ss=1 ; ss<=ssCnt ; ++ss ) {
+ if( (ssSubSspc[ss] == subss) && (ssBckYspc[ss] == subys) \
+ && (ssMethod[ss] == meth) ) {
+ return ss # found it
+ }
+ }
+ return 0 # not found
+}
+
+function ss_new_ini()
+{
+ return ss_new("ini", 0,0,0)
+}
+
+function ss_mk_dir(subss ,ss) # --> ss
+{
+ if( ssHasDir[subss] ) return subss
+ if( ss = ss_find("+dir", subss, 0) ) return ss
+ return ss_new("+dir", subss, 0, 1)
+}
+
+function ss_mk_bck(subss, subys ,ss) # --> ss
+{
+ if( ! ssHasDir[subss] ) do_exit(1)
+ if( ss = ss_find("+bck", subss, subys) ) return ss
+ return ss_new("+bck", subss, subys, ssHasDir[subss])
+}
+
+#------------------------------------------------------------------------------
+# Basic states are introduced by names. Our internal states are numeric > 0.
+# We use the prefix "sta" to indicate an internal numeric state.
+# - staCnt #(states): last number assigned to a state
+# - staSS[sta] state space of sta
+# - staNam[sta] string (to print), derived from basic state name
+# The mapping from staNam[sta] back to sta is unique only for basic states:
+# - basSta[nam] state number of basic state name
+# - staCla[sta] class: "normal" | "halting" | "funny"
+# - staBas[sta] basic state (internal sta)
+# - staSub[sta] 0 | lower state
+# - staDir[sta] 0 | direction/orientation
+# - staBck[sta] backward symbol (may be 0)
+# - staByDir[sta,dir] "" | state constructed by adding dir to sta
+# - staByBck[sta,bck] "" | state constructed by adding bck to sta
+# In order to enumerate states of a space:
+# - ssStaCnt[ss] #states in ss
+# - ssSta[ss,n] the n'th state in state space
+# FFS: detection of intermediate empty tape would want to recognize states
+# corresponding to start states (staBas==1?) and empty tape (staBck==0).
+# FFS: alternal state naming (Al, Ar) for printing outside the tape,
+# and indicating the head by just a "<" or ">". lrn = left/right name
+# This section does not contain any big/var numbers.
+
+function ss_new_sta(ss ,sta) # --> new sta, linked into ss
+{
+ sta = ++staCnt
+ staSS[sta] = ss
+ ssSta[ss, ++ssStaCnt[ss]] = sta
+ return sta
+}
+
+function sta_nam_dir(sta, dir) # --> nam
+{
+ if( dir < 0 ) return "<" staNam[sta]
+ if( dir > 0 ) return staNam[sta] ">"
+ return staNam[sta]
+}
+
+function sta_nam_bck(sta, bck, bckI ,n,s) # --> nam
+{
+ n = staNam[sta]
+ #s = bck
+ s = y_str(bckI, bck) #FFS: global flag
+ s = "(" s ")"
+ return (staDir[sta] <= 0) ? n s : s n
+}
+
+function sta_new_ini(ss,nam,cla ,sta) # --> new sta
+{
+ sta = ss_new_sta(ss)
+ basSta[nam] = sta # !
+ staNam[sta] = nam
+ staCla[sta] = cla
+ staBas[sta] = sta
+ staSub[sta] = 0
+ staDir[sta] = 0
+ #staByDir
+ #staByBck
+ return sta
+}
+
+function sta_new_funny(ss ,x) # --> new sta
+{
+ x = ss_new_sta(ss)
+ staNam[x] = "?"
+ staDir[x] = 0
+ staCla[x] = "funny"
+ staBas[x] = 0 # is in GLT, not here
+ staSub[x] = 0
+ return x
+}
+
+function sta_make(ss, sta, dir, bck ,m,x) # --> (old|new) sta
+{
+ if( staSS[sta] == ss ) return sta # direkt mapping
+ if( staSS[sta] != ssSubSspc[ss] ) do_exit(1)
+ m = ssMethod[ss]
+ if( m == "+dir" ) {
+ if( x = staByDir[sta,dir] ) return x # already constructed
+ staByDir[sta,dir] = x = ss_new_sta(ss)
+ staNam[x] = sta_nam_dir(sta, dir)
+ staDir[x] = dir
+ }else if( m == "+bck" ) {
+ if( x = staByBck[sta,bck] ) return x # already constructed
+ staByBck[sta,bck] = x = ss_new_sta(ss)
+ staNam[x] = sta_nam_bck(sta, bck, ssBckYspc[ss])
+ staDir[x] = staDir[sta]
+ staBck[x] = bck
+ }else do_exit(1) # unknown method | tries to create initial state
+ staCla[x] = staCla[sta]
+ staBas[x] = staBas[sta]
+ staSub[x] = sta
+ return x
+}
+
+# To start a higher machine we need the state corresponding to the basic
+# start state and an empty tape. The direction we use is arbitrary.
+
+function sta_make_empty(ss, bassta, dir ,m,sta) # --> (old|new) sta
+{
+ #printf("sta_make_empty(%s,%s,%s)\n",ss,bassta,dir)
+ m = ssMethod[ss]
+ if( m == "ini" ) return bassta
+ if( ss <= 0 ) do_exit(1)
+ if( (m == "+dir") || (m == "+bck") ) {
+ sta = sta_make_empty(ssSubSspc[ss], bassta, dir)
+ sta = sta_make(ss, sta, dir, 0)
+ }else do_exit(1)
+ return sta
+}
+
+function sta_ones(mID, sta ,ss,subss,submID,m,sum,bck,sumbck)
+{
+ #printf("sta_ones(%s,%s=%s)\n", mID,sta,staNam[sta])
+ if( staCla[sta] == "funny" ) return 0
+ if( mID <= mIDbas ) return 0
+ ## There is some machine below this one
+ submID = mYsubM[mID,1] # we need the first, only
+ ss = mSspc[mID]
+ subss = mSspc[submID]
+ if( subss == ss ) {
+ # Since the sub-machines share the state space, it is not
+ # important, which one we follow... we take the first one.
+ return sta_ones(submID, sta)
+ }
+ ## We are the lowest machine with this state space.
+ ## Follow the construction process to the lower state space:
+ m = ssMethod[ss]
+ if( m == "ini" ) return 0
+ if( (m == "+dir") || (m == "+bck") ) {
+ sum = 0
+ if( m == "+bck" ) {
+ bck = staBck[sta] # cannot be negative
+ sumbck = mYsum[submID,bck,0]
+ #printf("bck=%s (%s), sumbck=%s\n", bck, y_str(submID,bck), sumbck)
+ sum += sumbck
+ }
+ sum += sta_ones(submID, staSub[sta])
+ #printf("sta_ones(%s,%s=%s) --> %s\n", mID,sta,staNam[sta],sum)
+ return sum
+ }
+ do_exit(1)
+}
+
+#------------------------------------------------------------------------------
+# State mapping:
+# While we already know, how to (dis-)assemble states, we also need
+# to transfer states between lower machines running parts of a LT.
+# We avoid this mapping by forcing all lower machines to share the same Sspc.
+# [ This results from the limited way we use to "push" machines. ]
+
+#------------------------------------------------------------------------------
+# Symbol spaces:
+# Each new machine may construct a new symbol space, or may effectively
+# use the same symbols as a lower machine. Currently we do not construct
+# machines and symbol spaces separately: every mID introduces a new Yspc.
+# - Symbols are numbered locally to each mID, starting with 0.
+# - Symbol 0 is always the empty symbol (implicitly at the end of the tape)
+# - A symbol is constructed from a series of symbols of one or more basic mIDs.
+# This series is the LT used in transition computation.
+# - Symbols are constructed as needed. 0 is predefined.
+# - Basic machines have symbol numbers 0 and 1 (or more, cf. nBasSym).
+# - In our construction methods symbols within a Yspc are of fixed width.
+# The following builds a new Yspc for a machine and maps it to LTs
+# - mYcnt[mID] #symbols for mID (next number to define)
+# - mYnum[mID,conctap] symbol number for LT conctap (separator SUBSEP)
+# - mYsubY[mID,sym,1] first (left most) sub-symbol in LT for sym
+# - mYnsub[mID] #subsymbols in LTs for mID, 0 for initial
+# - mYnbas[mID] #(basic symbols) in a symbol for mID
+# - mYsubM[mID,1] mID for first subsymbol
+# - mYsum[mID,sym,dir] #nonzero in sym, at left / total / at right end
+# This is not the most general approach, but currently sufficient.
+# Note: bck-Machines have an additional symbol in their state for the LT.
+# This section does not contain any big/var numbers.
+
+function y_str(mID, sym ,n,i,str) # --> string (for printing)
+{
+ if( sym < 0 ) return gltTap[sym] # cf GLT-section below
+ n = mYnsub[mID]
+ if( n <= 0 ) return sym # directly basic symbol
+ str = ""
+ for( i=1 ; i<=n ; ++i ) {
+ str = str y_str(mYsubM[mID,i], mYsubY[mID,sym,i])
+ }
+ return str
+}
+
+#==============================================================================
+# Defining (pushing) higher machine definitions on top of lower,
+# according to wanted characteristics: This is mainly Yspc construction,
+# but also triggers state space construction.
+# - mType[mID] type of machine: "bas" | "mac" | "bck"
+# - mSspc[mID] state space for machine
+# For Macro construction we trace the factor and base
+# - mYsubK[mID] effective macro factor
+# - mYsubI[mID] on this machine
+# This section does not contain any big/var numbers.
+
+function mdef_0sym(mID ,n,i,key)
+{
+ n = mYnsub[mID]
+ if( n > 0 ) { # builds a Yspc
+ key = mID
+ for( i=1 ; i<=n ; ++i ) {
+ mYsubY[mID,0,i] = 0
+ key = key SUBSEP 0
+ }
+ mYnum[key] = 0
+ mYcnt[mID] = 1
+ for( i= -1 ; i<=1 ; ++i ) mYsum[mID,0,i] = 0 #FFS: unnecessary
+ }
+}
+
+function mdef_is_4T(mID ,trx) # --> bool FFS
+{
+ # exists a transition into a normal state, which does not move: yes
+ for( trx in trXmID ) {
+ if( trXmID[trx] == mID ) {
+ if( trAdd[trx] == 0 ) {
+ if( staCla[trSta[trx]] == "normal" ) {
+ return 1
+ }
+ }
+ }
+ }
+ return 0 # no, no 4T characteristic seen
+}
+
+function find_Ktyp(typ, K, onID ,i) # --> 0 | mID
+{
+ for( i=mIDbas ; i<=mIDtop ; ++i ) {
+ if( (mYsubI[i]==onID) && (mYsubK[i]==K) && (mType[i]==typ) ) return i
+ }
+ return 0
+}
+
+function push_macpair(onID1, onID2 ,mID) # combine 2 mac to mac
+{
+ if( !onID1 || !onID2 ) return 0 # cannot push onto undefined
+ if( mYsubI[onID1] != mYsubI[onID2] ) do_exit(1) #internal error
+ mID = ++mIDtop
+ mType[ mID] = "mac"
+ mSspc[ mID] = mSspc[onID1]
+ mYsubK[mID] = mYsubK[onID1] + mYsubK[onID2]
+ mYsubI[mID] = mYsubI[onID1]
+ mYnsub[mID] = 2
+ mYnbas[mID] = mYnbas[onID1] + mYnbas[onID2]
+ mYsubM[mID,1] = onID1
+ mYsubM[mID,2] = onID2
+ mdef_0sym(mID)
+ return mID
+}
+
+function mdef_K_typ(mID, typ, K, onID ,i) # --> echo mID
+{
+ if( mID ) {
+ mType[ mID] = typ
+ mYsubK[mID] = K
+ mYsubI[mID] = onID
+ mYnsub[mID] = K
+ mYnbas[mID] = K * mYnbas[onID]
+ for( i=1 ; i<=K ; ++i ) mYsubM[mID,i] = onID
+ mdef_0sym(mID)
+ }
+ return mID
+}
+
+function mdef_new_K_typ(ss, typ, K, onID ,mID) # --> new mID
+{
+ mID = ++mIDtop
+ mSspc[mID] = ss
+ return mdef_K_typ(mID, typ, K, onID)
+}
+
+function push5T(onID ,mID) # --> onID | 5T on 4T | 0
+{ #FFS: wieso machen wir das?
+ if( onID != mIDbas ) return onID
+ if( ! mdef_is_4T(onID) ) return onID
+ return mdef_new_K_typ(mSspc[onID], "bas", 1, onID)
+}
+
+function pushKmac(K, onID ,k1,k2,mID) # --> (old|new) mID
+{
+ if( !onID ) return 0
+ onID = push5T(onID) # FFS: always?
+ if( !onID ) return 0
+ if( mID = find_Ktyp("mac",K,onID) ) return mID # already there
+ if( K > 3 ) { # FFS: > 2
+ k1 = int(K/2)
+ k2 = K - k1
+ return push_macpair(pushKmac(k1,onID), pushKmac(k2,onID))
+ }
+ return mdef_new_K_typ( ss_mk_dir(mSspc[onID]), "mac", K, onID )
+}
+
+function pushFmac(K, onID) # F = flexible factor --> (old|new) mID
+{
+ if( !onID ) return 0
+ if( mType[onID] != "mac" ) return pushKmac(K, onID)
+ return pushKmac( K * mYsubK[onID], mYsubI[onID] )
+}
+
+function pushKbck(K, onID ,mID)
+{
+ if( !onID ) return 0
+ if( ! ssHasDir[mSspc[onID]] ) onID = pushKmac(1,onID)
+ if( mID = find_Ktyp("bck",K,onID) ) return mID # already there
+ return mdef_new_K_typ( ss_mk_bck(mSspc[onID], onID), "bck", K, onID )
+}
+
+function pushDup(onID ,mID) # --> (old|new) mID
+{
+ if( !onID ) return 0
+ mID = pushFmac(2, onID)
+ if( mType[onID] == "bck" ) mID = pushKbck(1, mID) # FFS
+ return mID
+}
+
+#FFS: before pushing: if( onID < mIDbas ) onID = mIDbas|mIDtop
+
+#==============================================================================
+# Transition coding:
+# Since we want to reference TRs of different machines, we create
+# a global registering of TRs with a numeric "transition index".
+# - trxCnt #TRs (last defined TRX)
+# - trInx[mId,sta,sym] 0 | TRX of transition
+# - trSym[trx] effective new symbol number
+# - trAdd[trx] effective head move (-1|1) [0 not for macro machines?]
+# - trSta[trx] effective new state (name)
+# - trRst[trx] "" | run state achieved by execution
+# Backwards:
+# - trXmID[trx] mID
+# - trXsta[trx] sta (state name)
+# - trXsym[trx] sym (local number)
+#
+# Running a LT may stop before "fell off tape", and fail to build an upper
+# proper state. Nevertheless we define an ITR for this, which shares the
+# properties of normal TRs (statitics etc), except:
+# - trSym[trx] is a GLT (created here)
+# - trAdd[trx] =0 (no head move)
+# - trSta[trx] a special new state of class "funny"
+# Currently, funny states and negative symbols are tightly coupled.
+# The GLT completely codes state and position on basic level, and the funny
+# state does not carry any information besides its Space and being funny.
+
+# Additional information (not crucial):
+# - trBasAdd[trx] basic head position movement
+# - trScnMin[trx] minimum basic tape position scanned
+# - trScnMax[trx] maximum ... both 0-based
+# - trBasStp[trx] basic steps simulated by transition (BIG)
+# With bck-machines the scan range can exceed the direct tape symbol.
+
+function trxNew(mID, sta, sym ,trx) # --> new TRX
+{
+ trx = ++trxCnt
+ trInx[mID,sta,sym] = trx
+ trXmID[trx] = mID
+ trXsta[trx] = sta
+ trXsym[trx] = sym
+ # Not yet defined while !(trx in trSym)
+ return trx
+}
+
+function trxSet(trx, nsym, nadd, nsta, rst)
+{
+ trSym[trx] = nsym
+ trAdd[trx] = nadd
+ trSta[trx] = nsta
+ trRst[trx] = rst
+}
+
+function trx_get(mID, sta, sym ,trx) # --> 0 | (old|new) TRX
+{
+ trx = trInx[mID, sta, sym]
+ if( ! trx ) {
+ trx = tr_defnew(mID, sta, sym) # implemented below
+ }
+ return trx
+}
+
+# Up to here we do not expect any big/var numbers, except in trBasStp[]
+#------------------------------------------------------------------------------
+# Wrapper functions for arithmetic with B/V (big numbers and variables)
+# We here also define some function we do not need ourselves, but are
+# needed by higher level machines.
+
+function bv_Vname(vx) # --> "V(vx)"
+{
+ return livVname(vx)
+}
+
+function bv_IsConst(x) # does not contain variable
+{
+ return livIsConst(x)
+}
+
+function bv_EQ0(x) # x==0 guaranteed
+{
+ return livEQ0(x)
+}
+
+function bv_GT0(x) # x>0 guaranteed
+{
+ return livGT0(x)
+}
+
+function bv_EQ(x,y) # x==y guaranteed
+{
+ return livEQ(x,y)
+}
+
+function bv_GEconst(x,y) # x>=y guaranteed
+{
+ return livGEconst(x,y)
+}
+
+function bv_add(x,y) # --> x + y
+{
+ return livAdd(x,y)
+}
+
+function bv_sub(x,y) # --> x - y
+{
+ return livSub(x,y)
+}
+
+function bv_mul(x,y) # --> x * y
+{
+ return livMul(x,y)
+}
+
+function bv_Subst(vnval,x) # --> x with vars V substituted by vnval[V]
+{
+ return livSubst(vnval,x)
+}
+
+function bv_Sum0lt(vn, lim, x) # --> sum(vn=0,vn<lim: x)
+{
+ return livSum0lt(vn, lim, x)
+}
+
+#..............................................................................
+
+function bv_OptShort(h)
+{
+ return bnSetOptLongHalf(h)
+}
+
+function bv_SHpush(h)
+{
+ bv__SHv[ ++bv__SHx ] = bnGetOptLongHalf()
+ return bnSetOptLongHalf(h)
+}
+
+function bv_SHpop()
+{
+ if( bv__SHx ) {
+ return bnSetOptLongHalf( bv__SHv[ bv__SHx-- ] )
+ }
+ return bnGetOptLongHalf()
+}
+
+function bv_P(x) ## --> print value (string)
+{
+ return bnShort(x)
+}
+
+#------------------------------------------------------------------------------
+# Transition statistics:
+# For non-basic TRs we keep track of which sub-TRs have been used how often
+# to build the TR (by running a LT). This data can be used to determine:
+# - the number of basic steps, one execution of this TR stands for
+# - the basic tape head movement, one execution of this TR stands for
+# FFS: isn't that possible more easily by counting macStep-types?
+#
+# A TR-statistics array (trcnt[]) is indexed by the TRXs.
+# Each non-basic TR has such a statistic, which must be enumerable:
+# - trUse[trx,subtrx] #uses of subtrx by execution of one trx B/V
+# - trUcnt[trx] #(different subtrx) for trx
+# - trUtrx[trx,1] first subtrx for trx
+
+function tr_add_use(trx, trcnt ,subtrx)
+{
+ # Add the tr-usage in "trcnt[]" into the statistic for "trx":
+ # trcnt[]: B/V
+ for( subtrx in trcnt ) {
+ if( ! bv_EQ0(trcnt[subtrx]) ) {
+ if( ! trUse[trx,subtrx] ) {
+ trUtrx[trx, ++trUcnt[trx]] = subtrx
+ }
+ trUse[trx,subtrx] = bv_add(bn0(trUse[trx,subtrx]), trcnt[subtrx])
+ }
+ }
+}
+
+# FFS: When running the top machine (on an unlimited tape), we also gather
+# such a statistic, which is used for printing steps and tape position.
+# For the top machine we gather two arrays, one with exponent multiplications
+# in it (effective execution), and one without (real executions).
+# - trCnt[trx] so often used explicitly (+= 1) B/V
+# - trSum[trx] so often used effectively (+= exponent) B/V
+# Each trSum may be explained by more elementary transitions.
+
+# Array trcnt[] is filled with (some) execution counts. Change it to
+# larger counts of more basic transitions, but not below mID.
+# Even mith multiple vars, we will not fail here.
+
+function trcnt_reduc(trcnt,mID ,ncnt,chg,trx,w,n,i,subtrx)
+{
+ # We avoid inserting into the scanned array, but delete current entry.
+ # trcnt[]: B/V
+ # NB: non-numeric and negative entries in trcnt are deleted also.
+ do {
+ chg = 0
+ for( trx in trcnt ) {
+ if( trx <= 0 ) { delete trcnt[trx]; continue }
+ w = trcnt[trx] # so often occurred trx
+ if( bv_EQ0(w) ) { delete trcnt[trx]; continue }
+ if( trXmID[trx] <= mID ) continue
+ n = trUcnt[trx] # so many subtrx for trx
+ if( n <= 0 ) continue # is a basic trx
+ delete trcnt[trx]
+ for( i=1 ; i<=n ; ++i ) {
+ subtrx = trUtrx[trx,i]
+ ncnt[subtrx] = bv_add(bn0(ncnt[subtrx]), \
+ bv_mul(w, trUse[trx,subtrx]))
+ chg = 1
+ }
+ }
+ if( ! chg ) break # nothing filled into ncnt[]
+ chg = 0
+ for( trx in ncnt ) {
+ w = ncnt[trx]; delete ncnt[trx]
+ if( w ) {
+ trcnt[trx] = bv_add(bn0(trcnt[trx]), w)
+ chg = 1
+ }
+ }
+ }while( chg )
+}
+
+function trcnt_basstp(trcnt ,sum,trx)
+{
+ sum = 0
+ for( trx in trcnt ) {
+ if( trx <= 0 ) continue
+ sum = bv_add(sum, bv_mul(trcnt[trx], trBasStp[trx]))
+ }
+ return sum
+}
+
+#------------------------------------------------------------------------------
+# Limited tape coding:
+# A LT is an array [1..n] of symbols with a state and position.
+# - limtap["len"] = n
+# - limtap["pos"] = 0 .. n+1, current symbol to read
+# - limtap["sta"] state (of machine to run this)
+# LTs can appear as the result of a transition: see GLT below.
+# FFS: scan range etc...
+# - limtap["baspos"]
+# - limtap["ipos"]
+# - limtap["ista"]
+# - limtap["scnmin"]
+# - limtap["scnmax"]
+
+#==============================================================================
+# Globally registered LTs (GLTs) as negative symbols:
+# The resulting LTs of an ITR are translated into special symbols with
+# global scope and negative value. They consist of basic symbols.
+# - gltCnt #(negative symbols)
+# - gltTap[y] string of basic symbols
+# - gltSta[y] (basic) state
+# - gltPos[y] position of head, starts with 1 (substr index)
+# - gltSum[y,dir] #nonzero at left / total / at right end
+# This section does not contain any big/var numbers.
+
+function glt_app_sym(y, mID, sym)
+{
+ if( sym < 0 ) {
+ gltSta[y] = gltSta[sym]
+ gltPos[y] = length(gltTap[y]) + gltPos[sym]
+ gltTap[y] = gltTap[y] gltTap[sym]
+ }else {
+ gltTap[y] = gltTap[y] y_str(mID,sym)
+ }
+}
+
+function glt_app_sta(y, sta, off ,m,bss)
+{
+ m = ssMethod[staSS[sta]]
+ if( m == "+bck" ) {
+ bss = ssBckYspc[staSS[sta]]
+ if( staDir[sta] >= 0 ) glt_app_sym(y, bss, staBck[sta])
+ glt_app_sta(y, staSub[sta], off)
+ if( staDir[sta] < 0 ) glt_app_sym(y, bss, staBck[sta])
+ }else if( m != "ini" ) {
+ glt_app_sta(y, staSub[sta], off)
+ }else {
+ gltSta[y] = sta
+ gltPos[y] = length(gltTap[y]) + off
+ }
+}
+
+function glt_app_stasym(y, mID, sta, sym)
+{
+ if( staDir[sta] >= 0 ) glt_app_sta(y, sta, 1)
+ glt_app_sym(y, mID, sym)
+ if( staDir[sta] < 0 ) glt_app_sta(y, sta, 0)
+}
+
+function glt_new(limtap, submid ,y,n,i,sum) # --> new neg sym
+{
+ ## NB: here we need basic symbols to be single chars.
+ y = -(++gltCnt)
+ n = limtap["len"]
+ for( i=1 ; i<=n ; ++i ) {
+ if( (i == limtap["pos"]) && (limtap[i] >= 0) ) {
+ glt_app_stasym(y, submid[i], limtap["sta"], limtap[i])
+ }else {
+ glt_app_sym(y, submid[i], limtap[i])
+ }
+ }
+ n = length(gltTap[y])
+ for( i=1 ; i<=n ; ++i ) {
+ gltSum[y, 0] += !! ( 0 + substr(gltTap[y], i, 1))
+ }
+ for( i=1 ; i<=n ; ++i ) {
+ gltSum[y,-1] += !! (sum = 0 + substr(gltTap[y], i, 1))
+ if( sum < 1 ) break
+ }
+ for( i=n ; i>0 ; --i ) {
+ gltSum[y, 1] += !! (sum = 0 + substr(gltTap[y], i, 1))
+ if( sum < 1 ) break
+ }
+ return y
+}
+
+#==============================================================================
+# Every state in every space can have run-state attributes (RST):
+# "" normal
+# "offtape" lt_1step: normal end of LT operation
+# "fail" lt_1step: could not retrieve/define transition
+# mt_DoStep: could not retrieve/define transition
+# "would halt" lt_1step: next steps would go into halt state (but !gohalt)
+# "stop" lt_1step: cannot continue with current state
+# mt_DoStep: cannot continue with current state
+# "E<=0: E" mt_DoStep: exponent not usable, anymore
+# "loop/insym" lt_run: loop inside LT detected
+# "loop/border" mt_PopN: pops Inf-exponent from border
+# A RST can also be filled in the initial trRst[].
+
+# FFS: remember last (normal) state/symbol?
+
+# - gohalt whether we really go into the halt state
+# otherwise we stop just before that transition
+
+#==============================================================================
+# Running some machine (mID > mIDbas) and a transition is not yet defined:
+# Go define it:
+# * set up limited tape and run submachine(s).
+# * find old or create new symbol from result tape
+# * construct my kind of state from substate, direction (and tape)
+# * define transition contents & return new trx
+# Go on to use the transition.
+#
+# Running a LT (and defining a new transition) never involves any Vars.
+# Big numbers may be involved for counting steps/transitions.
+
+function tr_defnew(mID, sta, sym \
+ ,trx,limtap,submid,trcnt,rst,nsym,nadd,nsta,dir,cut,bck) # --> 0|trx
+{
+ if(DBG) printf("> tr_defnew(%s, %s, %s) {\n", mID,sta,sym)
+ if( mID <= 1 ) return 0
+ if( sym < 0 ) return 0 #FFS: fatal?
+ # Check LIM(trx)
+ # Now will not fail, any more, but rather build an ITR
+ limtap[1] = 0 # mark it as array
+ submid[1] = 0 # mark it as array
+ trcnt[0] = 0 # mark it as array (FFS)
+ lt_build(mID, sta, sym, limtap, submid)
+ rst = lt_run(limtap, submid, trcnt)
+ if( DBG ) dump_arr(limtap, "res LT")
+ if( rst == "offtape" ) rst = "" # wants to go on, here
+ if( (rst == "") || lt_res_ok(limtap) ) { #FFS: nil-step pruefen
+ if( rst == "" ) {
+ nadd = lt_res_nadd(limtap) # new orientation
+ dir = nadd
+ cut = nadd
+ }else { # RST, but needs no GLT
+ nadd = 0 # == lt_res_nadd(limtap)
+ dir = staDir[sta]
+ cut = 0-dir
+ }
+ bck = 0
+ if( mType[mID] == "bck" ) bck = lt_cut_sym(limtap, cut)
+ nsta = sta_make(mSspc[mID], limtap["sta"], dir, bck)
+ nsym = lt_res_nsym(mID, limtap, submid)
+ }else {
+ nadd = 0
+ nsta = sta_new_funny(mSspc[mID])
+ nsym = glt_new(limtap, submid)
+ }
+ trx = trxNew(mID, sta, sym)
+ trxSet(trx, nsym, nadd, nsta, rst)
+ trBasAdd[trx] = limtap["baspos"] # was 0-based
+ trScnMin[trx] = limtap["scnmin"]
+ trScnMax[trx] = limtap["scnmax"]
+ trBasStp[trx] = trcnt_basstp(trcnt)
+ tr_add_use(trx, trcnt)
+ if(DBG) printf("< tr_defnew(%s, %s, %s) } %s (%s,%s,%s,%s)\n", mID,sta,sym,\
+ trx,nsym,nadd,nsta,rst)
+ return trx
+}
+
+function lt_build(mID, sta, sym, limtap, submid ,nbas)
+{
+ nbas = mYnbas[mID]
+ lt_fill(mID, sym, limtap, submid)
+ limtap["sta"] = ((mSspc[mID] == mSspc[submid[1]]) ? sta : staSub[sta])
+ limtap["pos"] = ((staDir[sta] < 0) ? mYnsub[mID] : 1)
+ if( mType[mID] == "bck" ) {
+ lt_ins_sym(limtap, staBck[sta], 0-staDir[sta])
+ submid[limtap["len"]] = submid[1]
+ nbas += mYnbas[ssBckYspc[mSspc[mID]]]
+ }
+ limtap["ipos"] = limtap["pos"]
+ limtap["ista"] = limtap["sta"]
+ limtap["baspos"] = 0
+ limtap["scnmin"] = nbas
+ limtap["scnmax"] = 0-nbas
+}
+
+function lt_fill(mID, sym, limtap, submid ,n,i)
+{
+ n = mYnsub[mID]
+ for( i=1 ; i<=n ; ++i ) {
+ limtap[i] = mYsubY[mID,sym,i]
+ submid[i] = mYsubM[mID, i]
+ }
+ limtap["len"] = n
+ if(DBG) { dump_arr(limtap, "fill LT"); dump_arr(submid, "fill subM") }
+}
+
+function lt_ins_sym(limtap,sym,dir ,pos)
+{
+ # NB: baspos & scan range are relative to the head and NOT changed here
+ if( dir == 0 ) return
+ pos = ++limtap["len"]
+ if( dir > 0 ) {
+ limtap[pos] = sym
+ }else { # dir < 0
+ limtap["pos"] += 1
+ while( --pos > 0 ) limtap[pos+1] = limtap[pos]
+ limtap[1] = sym
+ }
+}
+
+function lt_res_ok(limtap) # --> can build proper state
+{
+ if( ! (limtap["pos"] in limtap) ) return 1 # "offtape", normal
+ # try trAdd==0 special case.
+ # When orientation and position is the same as at the start,
+ # we can have a normal macro symbol, and a "stay" transition.
+ if( staCla[limtap["sta"]] != "funny" ) { # can build state
+ if( (limtap["pos"] == limtap["ipos"]) \
+ && (staDir[limtap["sta"]] == staDir[limtap["ista"]]) ) {
+ return 1
+ }
+ }
+ return 0
+}
+
+function lt_res_nadd(limtap ,pos) # --> mov
+{
+ pos = limtap["pos"]
+ if( pos < 1 ) return -1
+ if( pos in limtap ) return 0 #FFS: limtap["len"]
+ return 1
+}
+
+function lt_cut_sym(limtap,dir ,n,sym,i,xcut) # --> sym cut out
+{
+ # NB: baspos & scan range are relative to the head and NOT changed here
+ if( dir == 0 ) return 0 # FFS
+ n = limtap["len"]--
+ if( dir > 0 ) {
+ sym = limtap[xcut = n]
+ }else { # dir < 0
+ sym = limtap[xcut = 1]
+ for( i=2 ; i<=n ; ++i ) limtap[i-1] = limtap[i]
+ }
+ delete limtap[n]
+ if( limtap["pos"] > xcut ) limtap["pos"] -= 1
+ return sym
+}
+
+function lt_res_nsym(mID,limtap,submid ,key,i,n,nsym,sum) # --> (old|new) sym
+{
+ key = mID
+ n = limtap["len"]
+ for( i=1 ; i<=n ; ++i ) key = key SUBSEP limtap[i]
+ if( key in mYnum ) return mYnum[key] # already known symbol
+ mYnum[key] = nsym = mYcnt[mID]++
+ # NB: the following loops cannot be folded together
+ for( i=1 ; i<=n ; ++i ) {
+ mYsubY[mID,nsym,i] = limtap[i]
+ mYsum[mID,nsym, 0] += mYsum[submid[i],limtap[i],0] # native
+ }
+ for( i=1 ; i<=n ; ++i ) {
+ mYsum[mID,nsym,-1] += (sum = mYsum[submid[i],limtap[i],-1])
+ if( sum < mYnbas[submid[i]] ) break
+ }
+ for( i=n ; i>0 ; --i ) {
+ mYsum[mID,nsym, 1] += (sum = mYsum[submid[i],limtap[i], 1])
+ if( sum < mYnbas[submid[i]] ) break
+ }
+ return nsym
+}
+
+function lt_run(limtap, submid, trcnt, ltslow,k,n,rst)
+{
+ if(DBG) { dump_arr(limtap, "run LT"); dump_arr(submid, "run subM") }
+ # Loop detection: a "slow" tape does half as many steps. When looping
+ # both tape contents (syms+pos+sta) eventually will be equal.
+ # Worst case: the slow tape completes the cycle once.
+ # The slow tape contains only state, pos and symbols.
+ # To avoid problems with NIL-steps, we first do one normal step.
+ ltslow["pos"] = limtap["pos"]
+ ltslow["sta"] = limtap["sta"]
+ n = limtap["len"]
+ for( k=1 ; k<=n ; ++k ) ltslow[k] = limtap[k]
+ n = -1 # one addional initial step
+ while( "" == (rst = lt_1step(limtap, submid, trcnt)) ) {
+ if( ++n >= 2 ) {
+ n = 0
+ lt_slow_1step(ltslow, submid)
+ if( lt_slow_eq(ltslow, limtap) ) return "loop/insym"
+ }
+ }
+ return rst
+}
+
+function lt_1step(limtap, submid, trcnt ,pos,sym,sta,trx) # --> RST
+{
+ pos = limtap["pos"]
+ if( ! (pos in limtap) ) return "offtape" # fell off tape
+ sym = limtap[ pos ]
+ sta = limtap["sta"]
+ if( staCla[sta] != "normal" ) return "stop" # must not go on
+ # FFS check global credit (overall steps)
+ trx = trx_get(submid[pos], sta, sym)
+ if( ! trx ) return "fail"
+ lt_scn_update(limtap, trx)
+ if( !gohalt && staCla[trSta[trx]] == "halting" ) return "would halt"
+ trcnt[trx] = bv_add(bn0(trcnt[trx]), 1) # ++trcnt[trx]
+ limtap["baspos"] += trBasAdd[trx]
+ limtap[ pos ] = trSym[trx]
+ limtap["pos"] = trAdd[trx] + pos
+ limtap["sta"] = trSta[trx]
+ return trRst[trx]
+}
+
+function lt_slow_1step(limtap, submid ,pos,trx)
+{
+ # Guaranteed to be defined. No statistics.
+ pos = limtap["pos"]
+ trx = trInx[submid[pos], limtap["sta"], limtap[pos]]
+ limtap[ pos ] = trSym[trx]
+ limtap["pos"] = trAdd[trx] + pos
+ limtap["sta"] = trSta[trx]
+}
+
+function lt_slow_eq(ltslow, limtap ,k) # --> bool
+{
+ for( k in ltslow ) if( ltslow[k] != limtap[k] ) return 0
+ return 1
+}
+
+function lt_scn_update(limtap, trx ,p)
+{
+ # Since the positions in LTs are relative: only native numbers, here
+ if( trScnMin[trx] <= trScnMax[trx] ) {
+ p = limtap["baspos"] + trScnMin[trx]
+ if( p < limtap["scnmin"] ) limtap["scnmin"] = p
+ p = limtap["baspos"] + trScnMax[trx]
+ if( p > limtap["scnmax"] ) limtap["scnmax"] = p
+ }
+}
+
+#==============================================================================
+# Build internal machine (mID) from Input definitions:
+# Input machines have their TT implemented so:
+# - aSym[st,sym] effective new symbol (also 4T)
+# - aAdd[st,sym] effective head move (-1|0|1)
+# - aSta[st,sym] effective new state (name)
+# States above are always state names (not their numbers).
+# - nBasSym #(basic symbols) (minimum/default 2)
+# - nState #(numbered states) non-halt states
+# - stNum[st] ++nState; used to recognize halt states
+# - stNam[stnum] used for enumeration of states
+# We build:
+# - stSta[stnum] maps input state numbers to internal sta
+# stSta[stnum] == basSta[stNam[stnum]]
+
+function pushInitial( ss,mID,i,stnum,st,sym,nst,nsta,nsym,nadd,trx) # mIDbas
+{
+ mID = mIDbas = ++mIDtop
+ ss = ss_new_ini()
+ mType[mID] = "bas"
+ mSspc[mID] = ss
+ # add symbols...
+ mYcnt[ mID] = nBasSym #symbols for mID (next number to define)
+ mYnsub[mID] = 0 #subsymbols in LTs for mID, 0 for initial
+ mYnbas[mID] = 1 #(basic symbols) in a symbol for mID
+ for( sym=0 ; sym<nBasSym ; ++sym ) {
+ for( i= -1 ; i<=1 ; ++i ) mYsum[mID,sym,i] = !! sym
+ }
+ # add normal states...
+ for( stnum=1 ; stnum<=nState ; ++stnum ) {
+ stSta[stnum] = sta_new_ini(ss, stNam[stnum], "normal")
+ }
+ # add transitions & halting states...
+ for( stnum=1 ; stnum<=nState ; ++stnum ) {
+ st = stNam[stnum]
+ for( sym=0 ; sym<nBasSym ; ++sym ) {
+ nst = aSta[st,sym]
+ if( nst in basSta ) {
+ nsta = basSta[nst]
+ }else {
+ nsta = sta_new_ini(ss, nst, "halting")
+ }
+ nsym = aSym[st,sym] + 0
+ nadd = aAdd[st,sym]
+ trx = trxNew(mID, stSta[stnum], sym)
+ trxSet(trx, nsym, nadd, nsta, "")
+ trBasAdd[trx] = nadd
+ trScnMin[trx] = 0
+ trScnMax[trx] = 0
+ trBasStp[trx] = 1
+ }
+ }
+ return mID
+}
+
+#==============================================================================
+# Macro tapes:
+# They have exponents and consist of two stacks, e.g.:
+# [-1] [-2] [-3] <A [4] [3] [2] [1]
+# tape[ 1] The rightmost symbol
+# tape["E",1] The rightmost exponent B/V
+# tape[0] The orientation of the head (the active stack side): -1 | 1
+# tape["T",side] the index of the topmost elements in each stack
+# The side is empty, if topmost index == 0.
+# tape["mID"] The Yspc defining the symbols & the TM running this tape
+# tape["sta"] The current state
+# tape["rst"] The current run state
+
+# Exponents in tapes are always >= 0. Pure numeric exponents are > 0.
+# With Vars in exponents we avoid reading the symbol, when the exponent
+# may be zero, even if the operation is independent of the exponent. FFS.
+# With Vars in exponents, operations on macro tapes may fail, when the
+# exponents becomes too small. We concentrate the check for failure
+# in the function "mt_DoStep()".
+
+# Not used parts of the tape (popped off) are deleted. This is not strictly
+# necessary, but prepares easy copying of complete tapes.
+
+function mt_Init(tape, mID, sta)
+{
+ tape["T",-1] = 0
+ tape["T", 1] = 0
+ tape["mID"] = mID
+ tape["sta"] = sta
+ tape[0] = ((staDir[sta] < 0) ? -1 : 1)
+ # FFS: "rst"
+ #printf("mt_Init(tape,%s,%s): ori=%s\n", mID,sta, tape[0])
+}
+
+function mt_PopN(tape,dir ,x,expo) # --> popped expo B/V
+{
+ x = tape["T",dir]
+ if( x ) { # existing cell
+ expo = tape["E",x]
+ tape["T",dir] -= dir
+ delete tape[ x]
+ delete tape["E",x]
+ return expo
+ }else { # implicit border cell
+ tape["rst"] = "loop/border"
+ return 1000000 # sort of +Inf, FFS
+ }
+}
+
+function mt_Pop1(tape,dir ,x)
+{
+ x = tape["T",dir]
+ if( x ) { # existing cell
+ tape["E",x] = bv_add(tape["E",x], -1) # --tape["E",x]
+ if( bv_EQ0(tape["E",x]) ) {
+ tape["T",dir] -= dir
+ delete tape[ x]
+ delete tape["E",x]
+ }
+ }
+}
+
+function mt_PushN(tape,dir, sym,expo ,x)
+{
+ x = tape["T",dir]
+ if( (x ? tape[x] : 0) == sym ) { # unite old and new
+ if( x != 0 ) { # existing cell
+ tape["E",x] = bv_add(tape["E",x], expo)
+ }
+ }else { # need new tape cell
+ tape["T",dir] = (x += dir)
+ tape[ x] = sym
+ tape["E",x] = expo
+ }
+}
+
+function mt_Push1(tape,dir, sym)
+{
+ mt_PushN(tape,dir, sym,1)
+}
+
+function mt_Step(tape, dir,osym,osta, nsym,nsta,mov ,repcnt) #--> repcnt B/V
+{
+ if( mov != dir ) { # change dir | stay
+ if( nsym != osym ) {
+ mt_Pop1(tape,dir)
+ mt_Push1(tape,dir, nsym)
+ }
+ if( mov ) tape[0] = mov
+ repcnt = 1
+ }else { # step forward
+ if( nsta != osta ) { # step forward one
+ mt_Pop1(tape,dir)
+ mt_Push1(tape,0-dir, nsym)
+ repcnt = 1
+ }else { # step forward many
+ # here we could live with 0-exponents
+ # tape is empty at dir ==> border loop
+ repcnt = mt_PopN(tape,dir)
+ mt_PushN(tape,0-dir, nsym,repcnt)
+ }
+ }
+ tape["sta"] = nsta
+ # We don't do any statistics here. The caller does it his way.
+ return repcnt # so often executed transition
+}
+
+function mt_UpdStats(cinfo, trx, repcnt)
+{
+ ## Update config info according to repcnt executions of trx:
+ ## FFS: should we only update already existing components?
+ ## These computations tend to become really expensive for PA-CTRs.
+
+ cinfo["stpSoft"] = bv_add( bn0(cinfo["stpSoft"]), 1)
+ if( "stpEff" in cinfo ) {
+ cinfo["stpEff"] = bv_add(cinfo["stpEff" ], repcnt)
+ }
+ cinfo["stpBas" ] = bv_add( bn0(cinfo["stpBas" ]),
+ bv_mul(repcnt, trBasStp[trx]) )
+ cinfo["tposBas"] = bv_add( bn0(cinfo["tposBas"]),
+ bv_mul(repcnt, trBasAdd[trx]) )
+ #FFS: ScnMin, ScnMax
+ #FFS: trx usage
+ #FFS: attribute "intermediate empty"
+}
+
+function mt_TRX(tape ,dir,x,sym,sta) # --> 0 | known TRX to use
+{
+ if( tape["rst"] ) return 0
+ sta = tape["sta"]
+ if( staCla[sta] != "normal" ) return 0
+ dir = tape[0]
+ x = tape["T",dir]
+ # Here is the central point to check the exponent for usability.
+ if( x && ! bv_GT0(tape["E",x]) ) return 0
+ sym = (x ? tape[x] : 0)
+ return trInx[tape["mID"], sta, sym]
+}
+
+function mt_DoStep(tape,cinfo ,dir,x,sym,sta,trx,repcnt) # --> did it
+{
+ if( tape["rst"] ) return 0
+ sta = tape["sta"]
+ if( staCla[sta] != "normal" ) {
+ tape["rst"] = "stop"; return 0
+ }
+ dir = tape[0]
+ x = tape["T",dir]
+ # Here is the central point to check the exponent for usability.
+ if( x && ! bv_GT0(tape["E",x]) ) {
+ tape["rst"] = "E<=0: " tape["E",x]; return 0
+ }
+ sym = (x ? tape[x] : 0)
+ trx = trx_get(tape["mID"], sta, sym)
+ if( ! trx ) {
+ tape["rst"] = "fail"; return 0
+ }
+ repcnt = mt_Step(tape, dir,sym,sta, trSym[trx],trSta[trx],trAdd[trx])
+ tape["rst"] = trRst[trx]
+ mt_UpdStats(cinfo, trx, repcnt)
+ # Warning: the step may have just added an RST (without a real step)
+ return 1
+}
+
+function mt_Sum(tape ,sum,mID,dir,p,sym,ones) # --> non-0s in tape B/V
+{
+ sum = 0
+ mID = tape["mID"] # Yspc
+ ## Scan tape symbols, and multiply their ones with their exponent...
+ for( dir= -1 ; dir<=1 ; dir+=2 ) {
+ for( p=tape["T",dir] ; p ; p-=dir ) {
+ if( sym = 0 + tape[p] ) { # non-0 symbol
+ ones = ((sym < 0) ? gltSum[sym,0] : mYsum[mID,sym,0])
+ sum = bv_add(sum, bv_mul(tape["E",p], ones))
+ }
+ }
+ }
+ ## Do not forget the ones hidden in the "bck"-context of the state...
+ sum = bv_add(sum, sta_ones(tape["mID"], tape["sta"]))
+ return sum
+}
+
+#==============================================================================
+# Dumping: Print collected data low level (for debugging & non-believers :-)
+# Note: this section does not consider HTML. It just prints plain ASCII.
+
+function dump_arr_kv(arr, fmt ,k,kk)
+{
+ for( k in arr ) {
+ kk = k; gsub(SUBSEP, ",", kk)
+ printf(fmt, kk, arr[k])
+ }
+}
+
+function dump_arr_vk(arr, fmt ,k,kk)
+{
+ for( k in arr ) {
+ kk = k; gsub(SUBSEP, ",", kk)
+ printf(fmt, arr[k], kk)
+ }
+}
+
+function dump_arr(arr, txt)
+{
+ printf("ARR %s:", txt)
+ dump_arr_kv(arr, " [%s]=%s")
+ printf("\n")
+}
+
+#------------------------------------------------------------------------------
+
+function dump_1_SS(ss ,n,i)
+{
+ n = ssStaCnt[ss]
+ printf("SS %s: meth=%-4s subSS=%s bckYS=%-2s hasdir=%s #state=%s\n", \
+ ss, ssMethod[ss], ssSubSspc[ss], ssBckYspc[ss], ssHasDir[ss], n)
+ for( i=1 ; i<=n ; ++i ) {
+ printf("ss %s[%2d] = state %3s\n", ss, i, ssSta[ss,i])
+ }
+}
+
+function dump_all_SS( ss)
+{
+ for( ss=1 ; ss<=ssCnt ; ++ss ) dump_1_SS(ss)
+}
+
+function dump_1_ST(sta)
+{
+ printf("STA %2s %-7s: ss=%s cla=%-7s bas=%-2s sub=%-2s dir=%-2s bck=%s\n",\
+ sta, staNam[sta], staSS[sta], staCla[sta], staBas[sta], \
+ staSub[sta], staDir[sta], staBck[sta])
+}
+
+function dump_all_ST( sta)
+{
+ printf("STA CNT=%s\n", staCnt)
+ for( sta=1 ; sta<=staCnt ; ++sta ) dump_1_ST(sta)
+ dump_arr_kv(staByDir, "STA by dir %-7s = %s\n")
+ dump_arr_kv(staByBck, "STA by bck %-7s = %s\n")
+}
+
+function dump_1_M(mID ,n,i,y)
+{
+ printf("M %2s %3s ss=%-2s", mID, mType[mID], mSspc[mID])
+ printf(" Ycnt=%-3s Ynsub=%-3s Ynbas=%-3s subK=%2s subI=%s\n", \
+ mYcnt[mID], mYnsub[mID], mYnbas[mID], mYsubK[mID], mYsubI[mID])
+ n = mYnsub[mID]
+ if( n>=1 ) {
+ printf("\tsubM[1..%d] =", n)
+ for( i=1 ; i<=n ; ++i ) printf(" %3s", mYsubM[mID,i])
+ printf("\n")
+ for( y=0 ; y<mYcnt[mID] ; ++y ) {
+ printf("\tsym %6s =", y)
+ for( i=1 ; i<=n ; ++i ) printf(" %3s", mYsubY[mID,y,i])
+ printf("\t sums:")
+ for( i= -1 ; i<=1 ; ++i ) printf(" %3s", mYsum[mID,y,i])
+ printf("\n")
+ }
+ }
+}
+
+function dump_all_M( mID)
+{
+ printf("mIDbas=%s mIDtop=%s\n", mIDbas, mIDtop)
+ for( mID=mIDbas ; mID <= mIDtop ; ++mID ) dump_1_M(mID)
+ dump_arr_vk(mYnum, "Y %3s <-- %s\n")
+}
+
+function dump_1_TR(trx ,n,i,subtrx)
+{
+ printf("TRX %4s mID=%s sta=%-2s sym=%-3s --> ", \
+ trx, trXmID[trx], trXsta[trx], trXsym[trx])
+ printf("nsym=%-3s add=%2s nsta=%-2s rst=%s\n", \
+ trSym[trx], trAdd[trx], trSta[trx], trRst[trx])
+ printf("\t basadd=%3s scan=%3s..%-3s steps=%s\n", \
+ trBasAdd[trx], trScnMin[trx], trScnMax[trx], trBasStp[trx])
+ n = trUcnt[trx]
+ if( n>0 ) {
+ printf(" U %3d:", n)
+ for( i=1 ; i<=n ; ++i ) {
+ subtrx = trUtrx[trx,i]
+ printf(" %s*[%s]", trUse[trx,subtrx], subtrx)
+ }
+ printf("\n")
+ }
+}
+
+function dump_all_TR( trx)
+{
+ printf("TR CNT=%s\n", trxCnt)
+ for( trx=1 ; trx<=trxCnt ; ++trx ) dump_1_TR(trx)
+ dump_arr_vk(trInx, "trx %4s <-- %s\n")
+}
+
+function dump_1_GLT(y)
+{
+ printf("GLT %3s sta=%-2s pos=%-3s tap=%s\n", \
+ y, gltSta[y], gltPos[y], gltTap[y])
+}
+
+function dump_all_GLT( y)
+{
+ printf("GLT CNT = %s\n", gltCnt)
+ for( y=1 ; y<=gltCnt ; ++y ) dump_1_GLT(0-y)
+}
+
+function dump_ALL(what)
+{
+ printf("\nDUMP start: %s {\n", what)
+ dump_all_SS()
+ dump_all_ST()
+ dump_all_M()
+ dump_all_TR()
+ dump_all_GLT()
+ printf("DUMP complete} %s.\n", what)
+}
+
+#==============================================================================
+# Testing
+
+function mmTest( sta,trx)
+{
+ DBG=1
+ aSta[3,0] = "*"
+ pushInitial()
+ pushKmac(2, mIDbas)
+ dump_ALL("before")
+
+ sta = sta_make_empty(mSspc[mIDtop], stSta[1], 1)
+ trx = tr_defnew(mIDtop, sta, 0)
+ printf("sta = %s --> trx %s\n", sta, trx)
+ dump_ALL("after 1")
+
+ sta = trSta[trx]
+ trx = tr_defnew(mIDtop, sta, 0)
+ printf("sta = %s --> trx %s\n", sta, trx)
+ dump_ALL("after 2")
+}
+
+#==============================================================================
+# Auxiliary stuff for reading configurations (super monster 6)
+
+function ic__sym(str ,i)
+{
+ i = index(str, "^")
+ return i ? substr(str,1,i-1) : str
+}
+
+function ic__expo(str ,i)
+{
+ i = index(str, "^")
+ return i ? substr(str,i+1) : 1
+}
+
+function mk_symbol(mID, basstr ,submid,limtap,i,n,skiplen)
+{
+ #printf("mk_symbol(%s,%s)\n", mID,basstr)
+ if( mID <= mIDbas ) {
+ return int(basstr) # FFS
+ }
+ ## build a limtap by recursive calls
+ skiplen = 0
+ for( i=1 ; i<=mYnsub[mID] ; ++i ) {
+ submid[i] = mYsubM[mID,i]
+ n = mYnbas[submid[i]]
+ limtap[i] = mk_symbol(submid[i], substr(basstr, 1+skiplen, n))
+ skiplen += n
+ }
+ limtap["len"] = mYnsub[mID]
+ return lt_res_nsym(mID, limtap, submid)
+}
+
+function mk_state(mID, str ,ss,subss,m,dir,stsub,bck,i1,i2,s1,s2,s3)
+{
+ #printf("mk_state(%s,%s)\n", mID,str)
+ if( mID <= mIDbas ) return basSta[str]
+ ss = mSspc[mID]
+ subss = mSspc[mYsubM[mID,1]]
+ if( subss == ss ) {
+ return mk_state(mYsubM[mID,1], str)
+ }
+ ## we are the lowest machine with this state space
+ m = ssMethod[ss]
+ if( m == "ini" ) {
+ return basSta[str]
+ }
+ if( m == "+dir" ) {
+ if( str ~ "^<.*" ) {
+ dir = -1 ; str = substr(str,2)
+ }else if( str ~ ".*>$" ) {
+ dir = 1 ; str = substr(str,1,length(str)-1)
+ }else do_exit(1)
+ stsub = mk_state(subss, str)
+ return sta_make(ss, stsub, dir, 0)
+ }
+ if( m == "+bck" ) {
+ i1 = index(str, "(")
+ i2 = index(str, ")")
+ s1 = substr(str, 1, i1-1)
+ s2 = substr(str, i1+1, i2-i1-1)
+ s3 = substr(str, i2+1)
+ bck = mk_symbol(mYsubM[mID,1], s2)
+ stsub = mk_state(subss, s1 s3)
+ return sta_make(ss, stsub, dir, bck)
+ }
+ do_exit(1)
+}
+
+#==============================================================================
+# Module initialization
+
+function mmSim_init()
+{
+ if( mmSim_inited_ ) return # already done
+ mmSim_inited_ = 1 # committed to do so
+ # register us globally:
+ g_Id[++g_IdCnt] = "$Id: mmSim.awk,v 1.34 2005/01/09 22:23:28 heiner Exp $"
+
+ bignum_init()
+ varLI_init()
+}
+
+#------------------------------------------------------------------------------
+# End of module "mmSim"
+#==============================================================================
+##=============================================================================
+## $Id: tmJob.awk,v 1.34 2010/05/06 18:26:17 heiner Exp $
+## Author: Heiner Marxen
+## Module: tmJob -- read the job for TM simulation
+##
+##-----------------------------------------------------------------------------
+## The job we have to do is
+## - to simulate a Turing machine,
+## - to show (some of the) configurations
+## - to show summary information
+## - to show header & footer information
+## - to optionally present all that in HTML
+##
+## The job is coded as line oriented commands in the awk-input.
+## The job ist started/executed when the input is complete.
+##
+## NB: we support both, 4-tupel and 5-tupel machines (even mixes).
+##
+## Actions are added with their input coded positionally.
+## The action for A0 comes first, then A1, B0, B1 etc.
+## The action itself contains 2 or 3 items (characters).
+##
+## Arguments of type bool are coded numerically (0 == false, 1 == true).
+##
+## Input "commands" (linewise, awk fields/tokens):
+##
+## # comment
+## Is ignored.
+## / comment
+## Is ignored.
+## C comment
+## Recorded as comment in output.
+## : nonzerocount steps
+## Recorded as comment in output: "This TM produces ..."
+##
+## nbs count
+## Use count many basic symbols (at least 2)
+##
+## A stname act0 act1 ...
+## Add a new state, with actions on all basic symbols.
+## 5T tupel-list
+## Add more 5-tupel actions, states named A,B,C...
+## 5t tupel-list
+## Add more 5-tupel actions, states named 1,2,3...
+## 4T tupel-list
+## Add more 4-tupel actions, states named 1,2,3...
+## 4t tupel-list
+## Add more 4-tupel actions, states named 0,1,2...
+##
+## gohalt bool
+## Whether transition into halt state is executed.
+## M maxlines
+## The maximum number of lines (configurations) to print.
+## m maxsteps
+## The maximum number of steps (of basic TM) to simulate.
+## mtype factor ...
+## How to push macro machines. Factor == 0 is "bck" machine.
+## mmtyp N
+## How to run the top level macro machine:
+## 1 == normal
+## 2 == with pure additive config transitions (PA-CTRs)
+## 3 == also with PPA-CTRs
+## iniori dir
+## Start with an initial (state) orientation "dir".
+##
+## L leftoffset
+## Number of left tape cells to show at least (mini tapes only).
+## R leftround
+## Round up length of left tape part to multiple of this (mini tapes only)
+## r bool
+## Whether repetition reduction is done (mini machines only)
+## mac bool
+## Whether we do a 1-macro machine (old)
+##
+## H bool
+## Whether we do HTML output.
+## E bool
+## Whether exponents are shown as superscripts.
+## sympr style
+## Print symbols in tapes in that style ("i" == internal numbers)
+## show showvector...
+## List of numeric values which control which configs are shown.
+## short H
+## Shorten long numbers to H leading/trailing digits and a bracketed
+## middle length suppression indicator. Default H=0 suppresses this.
+##
+## T title ...
+## HTML title & main header
+## date ...
+## Date of page construction, documented in footer
+## edate ...
+## Marker to be replaced by final date, documented in footer
+## machv name text...
+## Add a named machine variant to the list.
+## iam name
+## Note our own machine variant.
+## pref URLprefix
+## Prefix for URLs, before the machine variant name.
+## HM bool
+## Whether we add links for Heiner's pages
+
+##=============================================================================
+## Helper functions
+
+function tj_ConcArgs(first ,r)
+{
+ r = $first
+ while( ++first <= NF ) r = r " " $first
+ return r
+}
+
+function ones_txt(suf)
+{
+ return ((nBasSym <= 2) ? "ones" : "nonzeros") suf
+}
+
+function tj_CmtRest(i) # whether field i and rest is comment
+{
+ if( substr($i,1,1) == "/" ) return 1 # yes, comment (ignored)
+ if( substr($i,1,1) == ":" ) {
+ if( (i+2) > NF ) return 1 # yes, comment (ignored)
+ ht_CmtAdd( sprintf("This TM produces %s %s in %s steps.",
+ $(i+1), ones_txt(), $(i+2)) )
+ return 1 # yes, comment registered
+ }
+ return 0 # no, nothing recognized
+}
+
+##=============================================================================
+## Construction of basic TM
+##
+## Here we accept various forms of TM construction commands, and fill some
+## variables and arrays coding all relevant data for "mmSim" to run the TM.
+##
+## State names as they appear in input, are directly used to index arrays
+## coding the basic transition table:
+## - aSym[st,sym] effective new symbol (0|1|...|9) (also 4T)
+## - aAdd[st,sym] effective head move (-1|0|1)
+## - aSta[st,sym] effective new state (name)
+##
+## - aIsy[st,sym] new symbol as coded in input
+## - aMov[st,sym] head move as coded in input
+## - aAct[st,sym] complete action as coded in input
+##
+## Those states, for which some transition is defined (cannot be a halt
+## transition), are numbered, starting with 1. All other states are
+## considered halt states.
+## - nBasSym #(basic symbols) (minimum/default 2)
+## - nState #(numbered states) non-halt states
+## - stNum[st] ++nState; used to recognize halt states
+## - stNam[stnum] state name (used for enumeration of states)
+
+## As support for 4-tupel machines we support '=' as code for head movement
+## (do not move) and 'i' or '.' for the new symbol (same as input symbol).
+## Normal left/right head moves are accepted as "LR", "-+" or "<>".
+
+## Note that we are not very strict in rejecting unexpected coding.
+
+function tj_SetNBS(nbs)
+{
+ nBasSym = 0 + nbs
+ if( nBasSym < 2 ) nBasSym = 2
+ if( nBasSym > 9 ) nBasSym = 9
+}
+
+function tj_isMovChr(c)
+{
+ return (length(c) == 1) && (0 != index("+-LR<>=", c))
+}
+
+function tj_mvAdd(mov) # map char --> head move (-1|0|+1)
+{
+ if( mov == "-" ) return -1
+ if( mov == "+" ) return 1
+ if( mov == "L" ) return -1
+ if( mov == "R" ) return 1
+ if( mov == "<" ) return -1
+ if( mov == ">" ) return 1
+ if( mov == "=" ) return 0
+ if( mov == -1 ) return -1 # allowes repeated mapping
+ if( mov == 1 ) return 1 # allowes repeated mapping
+ return 0
+}
+
+function tj_effMvAdd(mov) # map char --> effective head move (-1|0|+1)
+{
+ mov = tj_mvAdd(mov)
+ if( tjMapLR ) mov = 0-mov # for testing/debugging, only
+ return mov
+}
+
+function tj_effNewSym(osym,nsym ,y) # map char --> symbol (0|1|...)
+{
+ y = index("0123456789", nsym) - 1
+ if( (y >= 0) && (y < nBasSym) ) return nsym ## accept digits
+ if( nsym == "_" ) return "0" # accept "_" as blank (0)
+ if( nsym == " " ) return "0" # accept " " as blank (0)
+ if( nsym == "." ) return osym # accept "." as "same" (input)
+ if( nsym == "i" ) return osym # accept "i" as "same" (input)
+ if( nsym == "=" ) return osym # accept "=" as "same" FFS
+ return ""
+}
+
+function tj_mk_A(stnam,osym, nsym,nmov,nsta)
+{
+ aIsy[stnam,osym] = nsym # remember input code
+ aMov[stnam,osym] = nmov # remember input code
+
+ aSym[stnam,osym] = tj_effNewSym(osym, nsym)
+ aAdd[stnam,osym] = tj_effMvAdd(nmov)
+ aSta[stnam,osym] = nsta
+}
+
+## As action we support two different 5-tupel formats:
+## sym/mov/sta ala Green
+## sta/sym/mov ala MaBu
+## The reduced 4-tupel formats would be:
+## sym/sta ala Green
+## mov/sta ala Green
+## sta/sym ala MaBu
+## sta/mov ala MaBu
+## Since both, states and symbols, can be coded numerically, these can not
+## automatically be parsed correctly. Therefore, for the 4-tupel actions
+## (containing only 2 chars) we only support the MaBu coding.
+## Of course, an explicit 5-tupel code with embedded '='s is always ok.
+
+function tj_mkA5(stnam,sym,act ,nmov,nsym,nsta)
+{
+ aAct[stnam,sym] = act
+ if( length(act) < 3 ) do_exit(1)
+ if( tj_isMovChr(substr(act,2,1)) ) { # ala Green: sym/mov/sta
+ nsym = substr(act,1,1)
+ nmov = substr(act,2,1)
+ nsta = substr(act,3)
+ }else if( tj_isMovChr(substr(act,3,1)) ) { # ala MaBu: sta/sym/mov
+ nsta = substr(act,1,1)
+ nsym = substr(act,2,1)
+ nmov = substr(act,3)
+ }else {
+ do_exit(1)
+ }
+ tj_mk_A(stnam,sym, nsym,nmov,nsta)
+}
+
+function tj_mkA4(stnam,sym,act ,nmov,nsym,nsta)
+{
+ aAct[stnam,sym] = act
+ if( length(act) < 2 ) do_exit(1)
+ if( tj_isMovChr(substr(act,2,1)) ) { # sta/mov
+ nsta = substr(act,1,1)
+ nmov = substr(act,2,1)
+ nsym = sym
+ }else { # sta/sym
+ nsta = substr(act,1,1)
+ nsym = substr(act,2)
+ nmov = "="
+ }
+ tj_mk_A(stnam,sym, nsym,nmov,nsta)
+}
+
+##-----------------------------------------------------------------------------
+
+function tj_tupAdd1(act, tup, stBase ,num,nam,sym)
+{
+ sym = (tjTupCnt++ % nBasSym) # local toggle
+ if( sym == 0 ) { # first (input) symbol (of state)
+ num = ++nState # invent a new state (number & name)
+ # according to the base type convert the number into a state name
+ if( stBase == "1" ) {
+ nam = num
+ }else if( stBase == "0" ) {
+ nam = num-1
+ }else {
+ nam = substr("ABCDEFGHIJKLMNOPQRSTUVWXYZ",num,1)
+ }
+ stNam[num] = nam
+ stNum[nam] = num
+ }
+ if( tup == 5 ) {
+ tj_mkA5(stNam[nState], sym, act)
+ }else {
+ tj_mkA4(stNam[nState], sym, act)
+ }
+}
+
+function tj_tupAddStr(str, tup, stBase ,j)
+{
+ if( str == "" ) return
+ if( j = index(str, ",") ) { # has comma
+ if( j <= 2 ) { # del comma
+ tj_tupAddStr(substr(str,1,j-1) substr(str,j+1), tup, stBase)
+ }else { # split at comma
+ tj_tupAddStr(substr(str,1,j-1), tup, stBase)
+ tj_tupAddStr(substr(str, j+1), tup, stBase)
+ }
+ return
+ }
+ if( (tup == 4) && (length(str) == 4) ) { # double 4-tupel action
+ tj_tupAdd1(substr(str,1,2), tup, stBase)
+ tj_tupAdd1(substr(str,3,2), tup, stBase)
+ return
+ }
+ if( (tup == 5) && (length(str) == 6) ) { # double 5-tupel action
+ tj_tupAdd1(substr(str,1,3), tup, stBase)
+ tj_tupAdd1(substr(str,4,3), tup, stBase)
+ return
+ }
+ #FFS: triple action etc?
+ tj_tupAdd1(str, tup, stBase)
+}
+
+function tj_tupAdd1Arg(i, tup, stBase) # whether to stop (comment)
+{
+ if( tj_CmtRest(i) ) return 1 # yes, rest is comment
+ tj_tupAddStr($i, tup, stBase)
+ return 0 # no, go on
+}
+
+function tj_tupAddArgs(i, tup, stBase)
+{
+ for( ; i<=NF ; ++i ) {
+ if( tj_tupAdd1Arg(i, tup, stBase) ) break
+ }
+}
+
+##-----------------------------------------------------------------------------
+
+function tj_I_nAA( n,i)
+{
+ # a new state, explicitly named, with actions for all input symbols
+ n = ++nState ; tjTupCnt = 0
+ stNam[n] = $2
+ stNum[$2] = n
+ for( i=3 ; i<=NF ; ++i ) {
+ tj_mkA5($2, i-3, $i)
+ }
+}
+
+##-----------------------------------------------------------------------------
+## Macro machine construction
+## - tjMtype[1] first factor to push (0 == "bck")
+## - tjMtypeCnt so many factors remembered
+
+function tj_MacType(i ,j,k)
+{
+ for( j=i ; j<=NF ; ++j ) {
+ if( $j ~ "^#" ) break # allow comments
+ k = int($j)
+ if( k < 0 ) k = 0 # 0==bck
+ tjMtype[++tjMtypeCnt] = k
+ }
+}
+
+function tj_PushThem( j,k,mID) # --> 0 | ID of pushed top machine
+{
+ if( ! mIDtop ) {
+ pr("Pushing initial machine."); nl()
+ mID = pushInitial()
+ }else {
+ mID = mIDtop
+ }
+ for( j=1 ; j<=tjMtypeCnt ; ++j ) {
+ k = tjMtype[j]
+ if( k ) {
+ pr("Pushing macro factor " k "."); nl()
+ mID = pushKmac(k, mIDtop)
+ }else {
+ pr("Pushing BCK machine."); nl()
+ mID = pushKbck(1, mIDtop)
+ }
+ if( ! mID ) {
+ pr("*** Could not push machine."); nl()
+ break
+ }
+ }
+ return mID
+}
+
+function tj_RunType(mm)
+{
+ if( mm <= 1 ) {
+ doCH = 0
+ }else {
+ doCH = mm - 1
+ }
+}
+
+function tj_Set1ori(dir)
+{
+ tjIniOri = tj_mvAdd(dir)
+}
+
+function tj_1ori()
+{
+ if( tjIniOri ) return tjIniOri
+ return 1 ## default initial (state) orientation = R
+}
+
+##-----------------------------------------------------------------------------
+## Show vector
+## We want to show certain ranges with certain intervals, coded in an array:
+## [0] last value done
+## [1] next value to do: < 0: +inf (never)
+## [2] delta to use <=0: don't show
+## [3] up to this limit 0: unlimited
+## [4] [5] next delta, next limit ...
+##
+## - shBase[] global show vector
+## All numbers in this section are BIG numbers.
+
+function shv_Ini(shv, cur)
+{
+ #cur = 0+cur
+ shv[0] = bnDec(cur) # last
+ shv[1] = cur # next
+}
+
+function shv_Nxt(shv, cur ,r,b,j,stp,lim) # --> skip width (BIG)
+{ # did cur, set up next, ret skip width
+ #cur = 0+cur
+ r = bnSub(cur, shv[0])
+ #printf("cur - shv[0] = %s - %s = %s\n", cur,shv[0],r)
+ shv[0] = cur # ... and need new next
+ b = 0
+ for( j=2 ;; j += 2) {
+ stp = bn0(shv[j ])
+ lim = bn0(shv[j+1])
+ if( bnEQ0(lim) ) { # missing limit
+ if( bnLE0(stp) ) { # stp <= 0
+ shv[1] = -1 # unlimited nothing
+ }else {
+ #shv[1] = cur - ((cur-b) % stp) + stp
+ shv[1] = bnAdd(bnSub(cur, bnMod(bnSub(cur,b), stp)), stp)
+ }
+ #printf("1 stp/lim %s/%s --> nxt=%s\n",stp,lim,shv[1])
+ break;
+ }
+ if( bnLT(cur, lim) ) { # this limit to use
+ if( bnLE0(stp) ) { # no delta: next is limit
+ shv[1] = lim
+ }else {
+ #shv[1] = cur - ((cur-b) % stp) + stp
+ shv[1] = bnAdd(bnSub(cur, bnMod(bnSub(cur,b), stp)), stp)
+ }
+ #printf("2 stp/lim %s/%s --> nxt=%s\n",stp,lim,shv[1])
+ break;
+ }
+ b = lim # limit is next base
+ }
+ return r # skip width (BIG)
+}
+
+function shv_Frc(shv, cur) # --> 0 | skip width (BIG)
+{
+ #cur = 0+cur
+ if( bnEQ(cur, shv[0]) ) return 0 # did that last
+ return shv_Nxt(shv, cur)
+}
+
+function shv_Cur(shv, cur) # --> 0 | skip width (BIG)
+{
+ #cur = 0+cur
+ #printf("cur=%s shv[0/1]=%s/%s\n", cur,shv[0],shv[1])
+ if( bnLT(cur, shv[1]) ) return 0 # no, not yet
+ if( bnLT0(shv[1]) ) return 0 # no, unlimited
+ return shv_Nxt(shv, cur)
+}
+
+function shv_Set(j ,i,v)
+{
+ if( j <= NF ) {
+ for( i in shBase ) delete shBase[i]
+ for( i=j ; i <= NF ; ++i ) {
+ v = "" $i
+ if( bnIsOk(v) ) {
+ shBase[i] = v #FFS BIG
+ }else {
+ shBase[i] = 0 + v
+ }
+ #printf("shBase[%s] = %s (%s)\n", i, shBase[i], v)
+ }
+ }
+}
+
+##-----------------------------------------------------------------------------
+## - tjIL[] remember input (for debugging)
+## - tjILcnt
+
+## - supExpo whether exponents as superscripts
+
+function tj_Input()
+{
+ # The input line contains command input. Parse it.
+
+ tjIL[++tjILcnt] = $0 # just remember it
+ ht_Iline($0) # to document it in output
+
+ if( $1 == "C" ) { ht_CmtAdd(tj_ConcArgs(2)); return }
+ if( $1 == "#" ) { return }
+ if( tj_CmtRest(1) ) { return }
+
+ if( $1 == "nbs" ) { tj_SetNBS(0 + $2); return }
+ if( $1 == "A" ) { tj_I_nAA(); return }
+ if( $1 == "5T" ) { tj_tupAddArgs(2, 5, "A"); return }
+ if( $1 == "5t" ) { tj_tupAddArgs(2, 5, "1"); return }
+ if( $1 == "4T" ) { tj_tupAddArgs(2, 4, "1"); return }
+ if( $1 == "4t" ) { tj_tupAddArgs(2, 4, "0"); return }
+
+ if( $1 == "mtype" ) { tj_MacType(2); return }
+ if( $1 == "mmtyp" ) { tj_RunType(0+$2); return }
+ if( $1 == "iniori" ) { tj_Set1ori($2); return }
+ if( $1 == "gohalt" ) { gohalt = 0 + $2; return }
+
+ if( $1 == "H" ) { html = 0 + $2; return }
+ if( $1 == "E" ) { supExpo = 0 + $2; return }
+ if( $1 == "sympr" ) { symPrStyle = "" $2; return }
+ if( $1 == "T" ) { htTit = tj_ConcArgs(2); return }
+ if( $1 == "date" ) { htDate = tj_ConcArgs(2); return }
+ if( $1 == "edate" ) { htEdate = tj_ConcArgs(2); return }
+ if( $1 == "machv" ) { ht_AddVers($2, tj_ConcArgs(3)); return }
+ if( $1 == "iam" ) { htViam = $2; return }
+ if( $1 == "pref" ) { htPref = $2; return }
+ if( $1 == "HM" ) { HM = 0 + $2; return }
+ if( $1 == "mapLR" ) { tjMapLR = 0 + $2; return }
+
+ if( $1 == "M" ) { maxLines = 0 + $2; return }
+ if( $1 == "show" ) { shv_Set(2); return }
+ if( $1 == "short" ) { tjShortH = 0 + $2; return }
+
+ if( $1 == "bnspeed" ) {
+ bix_CheckAll()
+ if( (0+$2) >= 1 ) { bix_CheckNot("ioy"); }
+ if( (0+$2) >= 2 ) { bix_CheckNot("IO" ); }
+ return
+ }
+ ##FFS: reject unknown command
+}
+
+##=============================================================================
+## Print states
+
+function aug_for_st(txt, st, ashlt) # txt is quoted, already
+{
+ # "st" is an input state name
+ if( st in stNum ) { # normal (input) state
+ txt = htCCaug(txt, stNum[st])
+ }else { # halting state
+ if( ashlt == "" ) ashlt = 1
+ }
+ if( ashlt ) txt = htAaug(txt, hCAhalt)
+ return txt
+}
+
+function aug_st(st, ashlt)
+{
+ return aug_for_st(htQ(st), st, ashlt)
+}
+
+function aug_sta(sta, ashlt)
+{
+ # sta is an internal state number of "mmSim"
+ if( ashlt == "" ) ashlt = (staCla[sta] == "halting")
+ return aug_for_st(htQ(staNam[sta]), staNam[staBas[sta]], ashlt)
+ #FFS: ----------------^^^^^^^^^^^ may use [y] for bck's
+}
+
+##-----------------------------------------------------------------------------
+## Print basic transition table
+
+function tj_pr_basTT( i,j,st,tbtxt,ta,th,td,rs2,cs3)
+{
+ if( html ) {
+ tbtxt[-1] = "left"
+ tbtxt[ 0] = "stay"
+ tbtxt[ 1] = "right"
+ printf("</PRE><TABLE border>")
+ ta = "align=center"
+ th = "TH " ta
+ td = "TD " ta
+ rs2 = "rowspan=2"
+ cs3 = "colspan=3"
+ printf("<TR><%s %s>State</TH>\n", th,rs2)
+ for( j=0 ; j<nBasSym ; ++j ) {
+ printf("<%s %s>on<br>%d</TH>\n", th,rs2, j)
+ }
+ for( j=0 ; j<nBasSym ; ++j ) {
+ printf("<%s %s>on %d</TH>\n", th,cs3, j)
+ }
+ printf("</TR><TR>\n")
+ for( j=0 ; j<nBasSym ; ++j ) {
+ printf("<%s>Print</TH>\n", th)
+ printf("<%s>Move</TH>\n", th)
+ printf("<%s>Goto</TH>\n", th)
+ }
+ printf("</TR>\n")
+ for( i=1 ; i<=nState ; ++i ) {
+ st = stNam[i]
+ printf("<TR><%s>%s</TD>\n", td, aug_st(st,1))
+ for( j=0 ; j<nBasSym ; ++j ) {
+ printf("<%s>%s</TD>\n", td, htQ(aAct[st,j])) # FFS: H bold
+ }
+ for( j=0 ; j<nBasSym ; ++j ) {
+ printf("<%s>%s</TD>\n", td, htQ( aIsy[st,j] ))
+ printf("<%s>%s</TD>\n", td, htQ(tbtxt[aAdd[st,j]]))
+ printf("<%s>%s</TD>\n", td, aug_st( aSta[st,j],1))
+ }
+ printf("</TR>\n")
+ }
+ printf("<CAPTION align=bottom>Transition table</CAPTION>\n")
+ printf("</TABLE><PRE>\n")
+ }else {
+ printf("State")
+ for( j=0 ; j<nBasSym ; ++j ) printf(" sym=%d", j)
+ nl()
+ for( i=1 ; i<=nState ; ++i ) {
+ st = stNam[i]
+ printf("%5s", st)
+ for( j=0 ; j<nBasSym ; ++j ) {
+ printf(" %4s ", aAct[st,j])
+ }
+ nl()
+ }
+ nl()
+ }
+}
+
+##=============================================================================
+## Print config
+## - symPrStyle == "": normal, print what the symbol stands for (01 string)
+## == "i": internal, print the symbol numbers
+
+function tj_pr_cnf_1line(tape,style ,mID,dir,t,pa,pz,p,st,y,e)
+{
+ #printf("{{tj_pr_cnf_1line}}\n")
+ mID = tape["mID"]
+ if( style == "" ) {
+ if( supExpo <= 0 ) style = " "
+ else if( supExpo == 1 ) style = "^"
+ else style = "^b"
+ }
+ if( (style != " ") && !html ) style = " "
+ for( dir = -1 ; dir <= 1 ; ++dir ) { # left/head/right
+ if( dir == 0 ) { # head
+ #pr( " " staNam[tape["sta"]] )
+ # FFS: do not color the bck-symbol in the state?
+ t = tape["sta"]
+ if( staCla[t] != "funny" ) {
+ #FFS: use internal style for bck-sym(s) where from
+ # same symbol space
+ t = aug_sta(t)
+ t = htAaug(t, hCAhead)
+ printf(" %s", t)
+ }
+ }else {
+ t = tape["T",dir] # top index
+ if( t ) { # this side not empty
+ if( dir < 0 ) { # left side
+ pa = -1; pz = t
+ }else { # right side
+ pa = t; pz = 1
+ }
+ for( p=pa ; p>=pz ; --p ) {
+ # FFS: funny symbols (GLTs) have state in symbol?
+ y = tape[p]
+ if( y == "*" ) { # wildcard
+ pr(" [*]" )
+ }else if( y < 0 ) { # GLT-symbol
+ if( gltPos[y] > 1 ) {
+ pr(" " substr(gltTap[y],1,gltPos[y]-1))
+ }
+ st = staNam[gltSta[y]]
+ t = htQ(st ">")
+ t = aug_for_st(t, st, 1)
+ t = htAaug(t, hCAhead)
+ printf(" %s", t)
+ if( gltPos[y] <= length(gltTap[y]) ) {
+ pr(" " substr(gltTap[y],gltPos[y]))
+ }
+ }else { # normal symbol
+ if( symPrStyle == "" ) {
+ pr( " " y_str(mID, y) )
+ }else {
+ pr( " [" y "]" )
+ }
+ }
+ e = tape["E",p]
+ if( (e == "*") || ! bv_EQ(e, 1) ) {
+ if( length(e) > 1 ) {
+ e = bv_P(e)
+ }
+ if( style == " " ) {
+ pr("^" e)
+ }else {
+ prA(e, style)
+ }
+ }
+ }
+ }
+ }
+ }
+ if( tape["rst"] ) {
+ pr(" [" tape["rst"] "]")
+ }
+ # no NL
+}
+
+function tj_pr_symtab_warn()
+{
+ ## In case we have not printed the symbols meaning, but rather just
+ ## their number, we now have to explain the notation
+ if( symPrStyle == "" ) return
+ nl()
+ pr("Symbols are printed by their bracketed index into the symbol table,")
+ nl()
+ pr("which is printed below.")
+ nl()
+}
+
+function tj_pr_symtab(mID ,xlen,y,ycnt,yy)
+{
+ ## In case we have not printed the symbols meaning, but rather just
+ ## their number, we now have to list the meaning of these numbers.
+ if( symPrStyle == "" ) return
+ nl()
+ pr("Symbol Table:") ; nl()
+ ycnt = mYcnt[mID];
+ xlen = 1 ; yy = 10
+ while( (ycnt > yy) && (xlen < 6) ) {
+ xlen += 1 ; yy *= 10
+ }
+ for( y=0 ; y < ycnt ; ++y ) {
+ pr( str_bi(xlen+2, "[" y "]") " = " y_str(mID, y) )
+ nl()
+ }
+}
+
+function tj_PrConfHdr()
+{
+ pr(sprintf("%5s %*s %*s %s", "Steps", \
+ tjLstep,"BasSteps", \
+ tjLtpos,"BasTpos", "Tape contents"))
+ nl()
+}
+
+function tj_PrConf(tape,cinfo,style)
+{
+ pr( str_bi( 5, bv_P(cinfo["stpSoft"])))
+ pr(" " str_bi(tjLstep, bv_P(cinfo["stpBas" ])))
+ pr(" " str_bi(tjLtpos, bv_P(cinfo["tposBas"])))
+ pr(" ")
+ tj_pr_cnf_1line(tape, style)
+ nl()
+}
+
+##-----------------------------------------------------------------------------
+
+function cinfo_Init(cinfo)
+{
+ cinfo["stpBas" ] = 0
+ cinfo["stpSoft"] = 0 ## FFS: really needed? 1st Col in tj_PrConf()!
+ #cinfo["stpEff" ] = 0 ## FFS: unused?
+ cinfo["tposBas"] = 0 ## FFS: is expensive
+ delete cinfo[0] # has been used to arrify
+}
+
+function tj_pr_summary(tape,cinfo,lines,topsteps ,ones,b)
+{
+ nl()
+ ones = mt_Sum(tape)
+ pr( "Lines: " lines "\n")
+ pr( "Top steps: " topsteps "\n")
+ pr( "Macro steps: " cinfo["stpSoft"] "\n")
+ #pr("Macro steps: " cinfo["stpEff"] "\n")
+ pr( "Basic steps: " bv_P(cinfo["stpBas"]) "\n")
+ pr( "Tape index: " cinfo["tposBas"] "\n")
+ b = str_ba(9, ones_txt(":"))
+ pr( b " " bv_P(ones) "\n")
+ pr( sprintf("log10(%-8s): %8.3f\n", ones_txt(), bnLog10(ones)))
+ pr( sprintf("log10(%-8s): %8.3f\n", "steps", bnLog10(cinfo["stpBas"])))
+ if( ("rst" in tape) && ("" != tape["rst"]) ) {
+ pr(sprintf("Run state: %s\n", tape["rst"]))
+ }
+ b = bnDidShort()
+ if( "" != b ) {
+ nl()
+ pr("Some long numbers above are shortened: #(omitted digits) shown in ")
+ pr(sprintf("\"%s\".\n", b))
+ }
+ # FFS: MM-statistics:
+ # #(syms ) top/all
+ # #(states ) top/all
+ # #(transitions) top/all
+ # FFS: CTR-statistics
+}
+
+##=============================================================================
+## Config-Machines
+##
+## This is an experimental chapter. While the construction of good and fast
+## macro machines is well understood, the construction of "config machines"
+## is less well understood.
+##
+## A config machine (CM) builds a new sort of transition: config transitions
+## (CTRs). A CTR maps a complete configuration onto another one, not just
+## a single symbol under the head. Now, if it really were an exact complete
+## configuration, it would be pointless: when an exact configuration does
+## occur again, the TM is looping and does not halt.
+## (We will detect that, but it is not our main point)
+##
+## Of interest are similar configurations, which can be handled the same way.
+## Our main focus is on configurations, which differ only in some exponents
+## (tape repetition factors). We will replace some of the exponents by
+## named, symbolic variables, and operate with them for a fixed number of
+## normal macro steps. The resulting configuration still contains the
+## variables in the exponents.
+##
+## Such a step then is valid for all (non-negative) values of the variables.
+## Example:
+## a^(1+X) <S b^(0+Y) c after several steps
+## --> a^(0+X) <S b^(2+Y) c
+##
+## Those steps with variables within the exponents are valid only, as long
+## as they are independant of the exact value of the variables.
+## We hereby define, that all our variables always have non-negative integer
+## values, i.e. (Var >= 0) is true.
+## Therefore an exponent (1+X) is known to be >=1, and can be used for all
+## sorts of macro steps. Especially those macro steps walking through
+## a symbol without changing the exponent are helpful.
+##
+## The config mapping shown above is called a "pure additive" CTR (PA-CTR).
+## The resulting configuration is the same, except for some exponents
+## (typically two), which are changed by a constant amount. The variables
+## do not mix or change place.
+##
+## First, we will detect, remember and iterate PA-CTRs.
+##-----------------------------------------------------------------------------
+## Next, what happens between two PA-CTRs is tried to be catched
+## as one CTR, called a PPA-CTR: preparing a PA-CTR.
+##
+## As suffix appended to a PA-CTR we yield another state capable to do one,
+## and thus prepare a new type of machine with "states" consisting of
+## a complete tape except for some exponent variables, which are the "tape".
+
+##-----------------------------------------------------------------------------
+## There are many generalizations conceivable, but we do not yet do them.
+## Our CTRs are always made from a fixed amount of macro-TRs.
+##
+## The goal we try to achieve is: run the 5-state and even 6-state record
+## machines until they really halt.
+
+##=============================================================================
+## Config Transitions (CTRs)
+
+## - ctrxCnt #CTRs, last used index
+## - ctrxItape[ctrx,tx] input config to match (tape copy) (with wildcards)
+## - ctrxVcnt[ctrx] #(variable exponents) in input config, 1..n
+## - ctrxVtapx[ctrx,vx] tape index of var expo (in Itape, without the "E")
+## - ctrxVminE[ctrx,vx] minimum required expo value, variable offset
+## - ctrxOinfo[ctrx,cx] VarExpr for resulting cinfo[cx] (for repcnt=1)
+## - ctrxType[ctrx] "PA" | "PPA"
+##
+## ............................................................................
+## PA-CTRs have a output config nearly identical to its input.
+## + Symbols are identical.
+## + Variables do not move, and have a factor of one.
+## Hence we just have to add some constants to some exponents. (PA!)
+##
+## - ctrxVaddV[ctrx,vx] constant additive delta for PA variable
+##
+## ............................................................................
+## Other CTRs (PPAs) may have a completely different output tape.
+## It may contain wild cards, corresponding to the input tape values.
+## It may contain exponents with variables in them, at arbitrary places.
+## Such exponents may even contain multiple variables.
+## NB: it may be shorter or longer than the input tape.
+##
+## - ctrxOtape[ctrx,tx] output config that results (tape copy) (with wildcards)
+##
+## Since the tape indexes of the output tape are not easy to enumerate
+## from the above, we make a complete list of the tape indexes,
+## and what to do about them, when we execute such a step.
+## + Idel delete from Itape (does not occur in Otape)
+## + Icop identical in input and output
+## + Ocop copy from Otape (literally)
+## + Omfy take from Otape with Vars substituted
+## - ctrxAtyp[ctrx,tx] one of the above codes
+## - ctrxAtCnt[ctrx,atyp] #(tape index for atyp)
+## - ctrxAtLoc[ctrx,atyp,n] tape index of n'th atyp
+## That way we can enumerate the tape indexes for each type of action.
+
+function ctr__defAtyp(ctrx,tx,atyp ,x)
+{
+ ctrxAtyp[ctrx,tx] = atyp
+ x = ++ctrxAtCnt[ctrx,atyp]
+ ctrxAtLoc[ctrx,atyp,x] = tx
+ #printf("ctr__defAtyp(%s, %s, %s)\n", ctrx,tx,atyp)
+}
+
+function ctr_gNEW(itape,otape,cinfo,vinfo ,nv,ctrx,k,vx,tx)
+{
+ ## itape[] input tape, with WCs and vars
+ ## otape[] output tape, with WCs and vars
+ ## cinfo[] config info of otape[], based on itape[] zeros
+ ## General interface with vinfo[], containing (at least):
+ ## typ PA | PPA
+ ## nv #(vars)
+ ## vxtx[] maps var indexes to input tape indexes (without "E")
+ ## vxEmin[] maps var indexes to minimal exponent value (var offset)
+
+ nv = vinfo["nv"]
+
+ ctrx = ++ctrxCnt ## allocate a new CTR index
+ ctrxType[ctrx] = vinfo["typ"] ## and mark by type
+
+ ## Just copy input tape ...
+ for( k in itape ) {
+ ctrxItape[ctrx,k] = itape[k]
+ }
+ ## Just copy cinfo[] (may contain Vars) ...
+ for( k in cinfo ) {
+ ctrxOinfo[ctrx,k] = cinfo[k]
+ }
+
+ ## Build info about vars (in input tape) ...
+ ctrxVcnt[ctrx] = nv
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ tx = vinfo["vxtx",vx] ## without the "E"
+ ctrxVtapx[ctrx,vx] = tx ## without the "E"
+ ctrxVminE[ctrx,vx] = vinfo["vxEmin",vx]
+ }
+
+ if( "PA" == vinfo["typ"] ) { ## set up ctrxVaddV[]
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ tx = vinfo["vxtx",vx] ## without the "E"
+ k = bv_sub(otape["E",tx], itape["E",tx])
+ if( ! bnIsOk(k) ) {
+ pr("\n*** Internal error: ctr_NEW: " ctrx " delta " k); nl()
+ do_exit(1)
+ }
+ ctrxVaddV[ctrx,vx] = k
+ }
+ }
+
+ ## The following would not be strictly necessary for "PA".
+
+ ## Just copy output tape ...
+ for( k in otape ) {
+ ctrxOtape[ctrx,k] = otape[k]
+ }
+
+ for( k in itape ) {
+ if( ! (k in otape) ) {
+ ctr__defAtyp(ctrx, k, "Idel")
+ }
+ }
+ for( k in otape ) {
+ if( (k in itape) && (("" itape[k]) == ("" otape[k])) ) {
+ ctr__defAtyp(ctrx, k, "Icop")
+ }else { ## not in input, or changed in output
+ ## Variables must occur only in exponent positions.
+ ## FFS: how to detect reliably?
+ if( (substr(k,1,2) != ("E" SUBSEP)) || bv_IsConst(otape[k]) ) {
+ ctr__defAtyp(ctrx, k, "Ocop")
+ }else {
+ ctr__defAtyp(ctrx, k, "Omfy")
+ }
+ }
+ }
+
+ return ctrx
+}
+
+function ctr_NEW(typ,itape,otape,cinfo,nv,vxtx,vxEmin ,vx,vinfo)
+{
+ ## itape[] input tape, with WCs and vars
+ ## otape[] output tape, with WCs and vars
+ ## cinfo[] config info of otape[], based on itape[] zeros
+ ## nv #(vars)
+ ## vxtx[] maps var indexes to input tape indexes (without "E")
+ ## vxEmin[] maps var indexes to minimal exponent value (var offset)
+
+ vinfo["typ"] = typ
+ vinfo["nv" ] = nv
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ vinfo["vxtx" , vx] = vxtx[vx] ## without the "E"
+ vinfo["vxEmin", vx] = vxEmin[vx]
+ }
+ return ctr_gNEW(itape,otape,cinfo,vinfo)
+}
+
+function ctr_PA_check(ctrx,txadd ,nv,vx,tx)
+{
+ nv = ctrxVcnt[ctrx]
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ tx = ctrxVtapx[ctrx,vx] ## without the "E"
+ if( ctrxVaddV[ctrx,vx] != bn0(txadd[tx]) ) {
+ pr("\n*** Internal error: ctr_PA_check: " ctrx " vx=" vx)
+ pr(": " ctrxVaddV[ctrx,vx] " != " txadd[tx]) ; nl()
+ do_exit(1)
+ }
+ }
+}
+
+function ctr_PA_new(itape,otape,cinfo,nv,vxtx,vxE)
+{
+ return ctr_NEW( "PA", itape, otape, cinfo, nv, vxtx, vxE )
+}
+
+function ctr_PPA_new(itape,otape,cinfo,nv,vxtx,vxEmin)
+{
+ return ctr_NEW("PPA", itape, otape, cinfo, nv, vxtx, vxEmin)
+}
+
+function ctr_Match(tape,ctrx ,vx,vtx,k)
+{
+ ## Whether the tape[] matches the specified CTR's input pattern.
+ ## NB: we do not expect any wild cards in the tape[].
+ ## FFS: This function works for all types of CTRs (not just PAs).
+ ## FFS: PA-CTRs will never accept non-const values for decreasing vars.
+ ## But they may accept VarExprs for increasing vars.
+ ## PPAs may accept VarExprs for any Var (or WC).
+
+ #pr("/ match sta\n")
+ if( tape["sta"] != ctrxItape[ctrx,"sta"] ) return 0
+ #pr("/ match 0\n")
+ if( tape[ 0 ] != ctrxItape[ctrx, 0 ] ) return 0
+
+ ## Build a temporary array for the (full) tape indexes that may differ.
+ ## These are just the exponents replaced by Vars.
+ for( vx=1 ; vx <= ctrxVcnt[ctrx] ; ++vx ) {
+ vtx["E",ctrxVtapx[ctrx,vx]] = ctrxVminE[ctrx,vx]
+ }
+
+ for( k in tape ) {
+ if( k == "0" ) continue # see above (FFS)
+ #pr("/ match " k "\n")
+ if( k in vtx ) { # the input pattern has a var, here
+ #pr("/ match " tape[k] " >= " vtx[k] "\n")
+ if( ! bnGE(tape[k], vtx[k]) ) {
+ #pr("/ match FAIL\n")
+ return 0 # sorry
+ }
+ }else { # non-Var pos: WC or equal
+ ## FFS: check the ctrxItape[] elem to be there
+ ## NB: at non-var places wildcards may appear in the CTR
+ #pr("/ match " tape[k] " == " ctrxItape[ctrx,k] "\n")
+ if( "*" == ctrxItape[ctrx,k] ) continue
+ if( tape[k] != ctrxItape[ctrx,k] ) {
+ #pr("/ match FAIL\n")
+ return 0 # sorry
+ }
+ }
+ }
+ return 1 # no significant difference seen: matches
+}
+
+## When we expect to have many CTRs, we will better support searching
+## for a matching CTR. Currently we search linearly.
+
+function ctr_Search(tape ,ctrx,altctrx)
+{
+ ## NB: later we will pick the "best" (most productive)
+ altctrx = 0
+ for( ctrx=1 ; ctrx <= ctrxCnt ; ++ctrx ) {
+ #pr("/ search " ctrx " of " ctrxCnt "\n")
+ if( ctr_Match(tape, ctrx) ) {
+ if( "PA" == ctrxType[ctrx] ) {
+ return ctrx ## any PA is always wanted
+ }
+ altctrx = ctrx ## in case we cannot find PA
+ }
+ }
+ return altctrx ## may be a PPA
+}
+
+function ctr_do_1Step(ctrx,tape,cinfo ,vx,tx,vn,cx,d,vnval,k,atyp,i,n,ot)
+{
+ ## The ctrx has already been found to match.
+ ## FFS: we want to do this as part of another PPA definition,
+ ## i.e. with WCs and Vars in the tape!
+
+# if( "PA" == ctrxType[ctrx] ) {
+# pr("\n*** Internal error: ctr_do_1Step: " ctrx " is PA\n")
+# do_exit(1)
+# }
+ ## We shall execute it exactly once...
+ pr(sprintf("== Executing %3s-CTR %2s (once)", ctrxType[ctrx], ctrx))
+
+ ## First, we find the current values for the variables.
+ ## We subtract the minimum exponent value from the current exponent value.
+ ## We need them by name, not by index.
+ for( vx=1 ; vx <= ctrxVcnt[ctrx] ; ++vx ) {
+ tx = ctrxVtapx[ctrx,vx] ## Itape var position (without "E")
+ vn = bv_Vname(vx)
+ vnval[vn] = bnSub( tape["E",tx], ctrxVminE[ctrx,vx] )
+ pr(", " vn "=" bv_P(vnval[vn]))
+ }
+ nl()
+
+ ## Now we have to modify the tape. This is the hard part.
+ ## Coded as wild cards we may have a fixed outer part,
+ ## and a completely different inner part, including state, orientation
+ ## and top indexes of the HTs.
+ ## I.e., the result tape consists of
+ ## - some entries completely unchanged (WC in Itape & Otape) (copy Itape)
+ ## - some entries copied unchanged from the Otape (copy Otape)
+ ## - some entries constructed from the Otape by substituting (mfy Otape)
+ ## variables.
+ ## - some entries completely deleted (del Itape)
+
+ ## Handle "Idel" actions...
+ for( k in tape ) {
+ if( (ctrx SUBSEP k) in ctrxAtyp ) {
+ atyp = ctrxAtyp[ctrx, k]
+ if( "Idel" == atyp ) {
+ delete tape[k]
+ }
+ }else {
+ ## Non-expected index in tape. FFS: leave alone
+ pr("Unexpected: tape[" k "] = \"" tape[k] "\""); nl()
+ }
+ }
+ ## Actions "Icop" effectively are NOOPs...
+ ## Now, Idel and Icop are handled.
+
+ ## Handle "Ocop" actions...
+ atyp = "Ocop"
+ n = ctrxAtCnt[ctrx,atyp]
+ for( i=1 ; i<=n ; ++i ) {
+ k = ctrxAtLoc[ctrx,atyp,i]
+ tape[k] = ctrxOtape[ctrx,k]
+ }
+
+ ## Handle "Omfy" actions...
+ atyp = "Omfy"
+ n = ctrxAtCnt[ctrx,atyp]
+ for( i=1 ; i<=n ; ++i ) {
+ k = ctrxAtLoc[ctrx,atyp,i]
+ ot = ctrxOtape[ctrx,k]
+ ## FFS: vnval itself may contain vars, which shall be considered
+ ## constant with regard to this substitution!
+ d = bv_Subst(vnval, ot)
+ #pr("Subst " ot " --> " d " for [" k "]"); nl()
+ tape[k] = d
+ }
+
+ ## Now we update the cinfo[]: this needs the var values substituted
+ ## into the expressions, and evaluated.
+ for( cx in cinfo ) {
+ d = bv_Subst(vnval, ctrxOinfo[ctrx,cx])
+ #pr("Subst " ctrxOinfo[ctrx,cx] " --> " d "\n")
+ cinfo[cx] = bnAdd(cinfo[cx], d)
+ }
+}
+
+function ctr_StepN(ctrx,tape,cinfo ,nv,vx,vn,vnval,tx,cx,d,k,k1,yn,dneg,dpos,r)
+{
+ ## The ctrx has already been found to match, and to be PA.
+ ## Now perform the step as often as possible (at least once).
+
+ if( "PA" != ctrxType[ctrx] ) {
+ pr("\n*** Internal error: ctr_StepN: " ctrx " is not PA\n")
+ do_exit(1)
+ }
+ pr(sprintf("== Executing %3s-CTR %2s", "PA", ctrx))
+ ## First, we find the current values for the variables.
+ ## We need them by name, not by index.
+ nv = ctrxVcnt[ctrx]
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ tx = ctrxVtapx[ctrx,vx]
+ vn = bv_Vname(vx)
+ vnval[vn] = bnSub( tape["E",tx], ctrxVminE[ctrx,vx] )
+ pr(", " vn "=" bv_P(vnval[vn]))
+ }
+ ## Now, find the maximum allowed repetition count. It is limited by
+ ## each of the decreasing exponents, simultaneously...
+ dneg = ""
+ dpos = ""
+ k = "Inf" # Unlimited
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ d = ctrxVaddV[ctrx,vx]
+ if( bnLT0(d) ) {
+ ## The Var-value must not be negative before the step.
+ ## Thus we can do (V div d) + 1
+ vn = bv_Vname(vx)
+ k1 = bnInc( bnDiv(vnval[vn], bnNeg(d)) )
+ if( bnLT(k1,1) ) {
+ pr("\n*** Internal error: ctr_StepN: " k1 " steps for " vn "\n")
+ do_exit(1)
+ }
+ if( (k == "Inf") || bnLT(k1,k) ) { # currently most limiting
+ k = k1
+ dneg = bnNeg(d)
+ }
+ }else if( bnGT0(d) ) {
+ if( (dpos == "") || bnGT(d,dpos) ) {
+ dpos = d
+ }
+ }
+ }
+
+ ## Ok, we cannot really perform infinitely many steps.
+ ## That clearly is a run state...
+ if( k == "Inf" ) {
+ r = k " (doing just 1)"
+ k = 1 # just to show the effect
+ tape["rst"] = "loop/gtConfig" # greater config
+ }else {
+ r = bv_P(k)
+ }
+
+ pr(", repcount=" r)
+ if( (dpos != "") && (dneg != "") ) pr(", factor=" dpos "/" dneg)
+ nl()
+
+ ## Now update the changing exponents in the tape.
+ ## This is especially easy, not even the var values are needed,
+ ## we know the deltas by var/position
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ tx = ctrxVtapx[ctrx,vx]
+ d = ctrxVaddV[ctrx,vx]
+ d = bnMul(d, k) # k times add fix amount
+ tape["E",tx] = bnAdd(tape["E",tx], d)
+ #FFS: reducing to zero
+ }
+
+ ## Now we update the cinfo[]: this needs the var values substituted
+ ## into the expressions, and evaluated, and then summed over k steps.
+ ## Since all Vars are expected to change with each step, here we
+ ## have to be careful to do it right.
+ ## First, we express all var-values as a linear expression in
+ ## the same new variable, which then is summed from 0 to k-1.
+ ## V = Vinit + Vadd*Y 0 <= Y < k
+ yn = bv_Vname(nv+1) ## a new name
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ vn = bv_Vname(vx)
+ d = ctrxVaddV[ctrx,vx]
+ vnmap[vn] = bv_add(vnval[vn], bv_mul(d, "0+" yn))
+ }
+ for( cx in cinfo ) {
+ d = bv_Subst(vnmap, ctrxOinfo[ctrx,cx])
+ #pr("Subst " ctrxOinfo[ctrx,cx] " --> " d "\n")
+ d = bv_Sum0lt(yn, k, d)
+ #pr("Sum --> " d "\n")
+ cinfo[cx] = bnAdd(cinfo[cx], d)
+ }
+}
+
+##=============================================================================
+## Config History
+## - doCH 0=MM, 1=PA, 2=PPA
+## - chCnt last valid index (1 = first)
+## - chBad last invalid index (search block)
+## - chTap[chx,y] saved data from tape[y] for index chx
+## - chByx[chx] 0 | CTR index creating this config
+
+function ch_Save(tape,cinfo,byctrx ,chx,k) # --> new chx
+{
+ chx = ++chCnt
+ #printf("CH: save to %s\n", chx)
+ for( k in tape ) {
+ chTap[chx,k] = tape[k]
+ }
+ chByx[chx] = byctrx
+ return chx
+}
+
+function ch_SetBad(chx ,k)
+{
+ #printf("CH: setbad %s: Bad=%s Cnt=%s\n", chx,chBad,chCnt)
+ if( chCnt ) { # something in there
+ if( chBad < chx ) { # will increase
+ chBad = chx
+ ## Although not strictly necessary, we try to get rid
+ ## of all the remembered data, which we will not access, again.
+ ## It does not seem to speed us up, but it helps to
+ ## compute really large jobs by limiting the memory usage.
+ #printf("CH: ? %s >= %s\n", chBad, chCnt)
+ if( chBad >= chCnt ) { # nothing left, logically
+ #delete chTap # clear completely
+ ## deleting a complete array at once is not portable "awk".
+ for( k in chTap ) delete chTap[k]
+ #printf("CH: DELETED!\n")
+ }
+ }
+ }
+}
+
+function ch_SetFirst(chx)
+{
+ #printf("CH: setfirst %s\n", chx)
+ ch_SetBad(chx-1)
+}
+
+function ch_ToTape(chx,tape ,lo,hi,x)
+{
+ tape["mID"] = chTap[chx,"mID"]
+ tape["sta"] = chTap[chx,"sta"]
+ tape["rst"] = chTap[chx,"rst"]
+ tape[ 0 ] = chTap[chx, 0 ]
+ lo = chTap[chx,"T",-1]
+ hi = chTap[chx,"T", 1]
+ tape["T",-1] = lo
+ tape["T", 1] = hi
+ for( x=lo ; x <= hi ; ++x ) {
+ if( x == 0 ) continue
+ tape[ x] = chTap[chx, x]
+ tape["E",x] = chTap[chx,"E",x]
+ }
+}
+
+function ch_Pr(chx,txt ,tape) # for debugging
+{
+ if( txt ) printf("%s ", txt)
+ printf("chx %3d: ", chx)
+ tape[0] = 0 # arrify
+ ch_ToTape(chx, tape)
+ tj_pr_cnf_1line(tape)
+ nl()
+}
+
+function ch_Match(nchx,ochx ,lo,hi,x) # --> whether equal except exponents
+{
+ if( chTap[nchx,"sta"] != chTap[ochx,"sta"] ) return 0
+ if( chTap[nchx,"T",-1] != chTap[ochx,"T",-1] ) return 0
+ if( chTap[nchx,"T", 1] != chTap[ochx,"T", 1] ) return 0
+ if( chTap[nchx,"mID"] != chTap[ochx,"mID"] ) return 0
+ if( chTap[nchx,"rst"] != chTap[ochx,"rst"] ) return 0 #FFS
+ if( chTap[nchx, 0 ] != chTap[ochx, 0 ] ) return 0
+ lo = chTap[nchx,"T",-1]
+ hi = chTap[nchx,"T", 1]
+ for( x=lo ; x <= hi ; ++x ) {
+ if( x == 0 ) continue
+ if( chTap[nchx,x] != chTap[ochx,x] ) return 0
+ ## NB: we do NOT compare any exponents, here!
+ }
+ return 1
+}
+
+function ch_Tmatch(tape,ochx ,lo,hi,x) # --> whether equal except exponents
+{
+ ## The tape may contain wildcards!
+ if( tape["mID"] != chTap[ochx,"mID" ] ) return 0
+ if( tape["sta"] != chTap[ochx,"sta" ] ) return 0
+ if( tape["rst"] != chTap[ochx,"rst" ] ) return 0 #FFS
+ if( tape[ 0 ] != chTap[ochx, 0 ] ) return 0
+ if( tape["T",-1] != chTap[ochx,"T",-1] ) return 0
+ if( tape["T", 1] != chTap[ochx,"T", 1] ) return 0
+ lo = tape["T",-1]
+ hi = tape["T", 1]
+ for( x=lo ; x <= hi ; ++x ) {
+ if( x == 0 ) continue
+ if( tape[x] == "*" ) continue # wildcard
+ if( tape[x] != chTap[ochx,x] ) return 0
+ ## NB: we do NOT compare any exponents, here!
+ }
+ return 1
+}
+
+function ch_Find(nchx ,ochx)
+{
+ ## FFS: -2 ??
+ for( ochx = nchx-1 ; ochx > chBad ; --ochx ) {
+ if( ch_Match(nchx,ochx) ) {
+ #ch_Pr(ochx, "\tMatch")
+ return ochx # first match backwards
+ }
+ }
+ return 0 # not found
+}
+
+##-----------------------------------------------------------------------------
+
+function ch_Deltas(ochx,nchx,txadd,txser ,lo,hi,x)
+{
+ ## The configs did match, already
+ for( x in txadd ) {
+ if( x == 0 ) continue
+ delete txadd[x]
+ }
+ txadd[0] = 0
+ txser[0] = 0
+ lo = chTap[nchx,"T",-1]
+ hi = chTap[nchx,"T", 1]
+ for( x=lo ; x <= hi ; ++x ) {
+ if( x == 0 ) continue
+ if( chTap[nchx,"E",x] != chTap[ochx,"E",x] ) {
+ txadd[0] += 1
+ txadd[x] = bnSub( chTap[nchx,"E",x], chTap[ochx,"E",x] )
+ #printf("txadd[%d] = %s\n", x, txadd[x])
+ txser[ txser[0] += 1 ] = x
+ }
+ }
+ return txadd[0]
+}
+
+function ch_Wild_1R(chx ,t) # how many cells from right are unread
+{
+ if( chByx[chx] ) return 0 #FFS
+ t = chTap[chx,"T",1] # top index right side
+ if( t > 0 ) {
+ if( chTap[chx,0] > 0 ) --t # looks at that one: reduce
+ }
+ return t
+}
+
+function ch_Wild_1L(chx ,t) # how many cells from left are unread
+{
+ if( chByx[chx] ) return 0 #FFS
+ t = 0 - chTap[chx,"T",-1] # |top index left side|
+ if( t > 0 ) {
+ if( chTap[chx,0] < 0 ) --t # looks at that one: reduce
+ }
+ return t
+}
+
+function ch_WildR(ochx,nchx ,w,t) # how many cells from right are unread
+{
+ w = ch_Wild_1R(ochx)
+ while( (w > 0) && (++ochx <= nchx) ) {
+ t = ch_Wild_1R(ochx)
+ if( w > t ) w = t # collect minimum
+ }
+ return w
+}
+
+function ch_WildL(ochx,nchx ,w,t) # how many cells from left are unread
+{
+ w = ch_Wild_1L(ochx)
+ while( (w > 0) && (++ochx <= nchx) ) {
+ t = ch_Wild_1L(ochx)
+ if( w > t ) w = t # collect minimum
+ }
+ return w
+}
+
+function ch_WCpatch(tape,wcL,wcR ,t)
+{
+ #FFS: hat nix mit "ch" zu tun!
+ for( t=1 ; t <= wcR ; ++t ) {
+ tape["E",t] = "*"
+ tape[ t] = "*"
+ }
+ for( t=1 ; t <= wcL ; ++t ) {
+ tape["E",-t] = "*"
+ tape[ -t] = "*"
+ }
+}
+
+function ch_Vattach(tape,nv,vxtx ,vx,tx)
+{
+ #FFS: hat nix mit "ch" zu tun!
+ ## Attach Vars to the exponents at the variable positions
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ tx = vxtx[vx]
+ tape["E",tx] = tape["E",tx] "+" bv_Vname(vx)
+ #pr("Tag Var " vx " to TX " tx "\n")
+ }
+}
+
+function ch_Vrepl(tape,nv,vxtx,vxE ,vx,tx)
+{
+ #FFS: hat nix mit "ch" zu tun!
+ ## Replace the exponents at the variable positions
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ tx = vxtx[vx]
+ tape["E",tx] = vxE[vx]
+ }
+}
+
+function ch_Vrun(ochx,nchx,tape,cinfo,show) # --> succ
+{
+ ## The tape has variables in it, and maybe some wildcards.
+ ## Try to run it, along the saved config history ochx..nchx
+ ## This function does not only work for PA!
+
+ while( ochx < nchx ) {
+ ## FFS: chByx[ochx+1] --> do a PPA
+ #pr("ch_Vrun: chx " ochx);nl()
+ if( ! mt_TRX(tape) ) {
+ if( show>0 ) pr("<< Cannot determine macro-TR, here!\n")
+ return 0
+ }
+ if( ! mt_DoStep(tape, cinfo) ) {
+ if( show>0 ) pr("<< Cannot execute macro-TR, here!\n")
+ return 0
+ }
+ ++ochx # just did a step
+ if( show>0 ) tj_PrConf(tape, cinfo)
+
+ if( tape["rst"] || (staCla[tape["sta"]] != "normal") ) {
+ if( show>0 ) pr("<< Will not accept this, here!\n")
+ return 0
+ }
+ # Check, whether we look like the saved one:
+ if( ! ch_Tmatch(tape, ochx) ) {
+ if( show>0 ) pr("<< That does not match the saved config!\n")
+ return 0
+ }
+ }
+ return 1 # yes, did run it completely with some expos
+}
+
+##-----------------------------------------------------------------------------
+## Check termination condition after some run with vars along a saved history.
+##
+## Since we ran along a saved history, we know that ch_Tmatch() succeeded,
+## all the way along, i.e.:
+## - all non-WC symbols in the tape matched literally
+## - all exponents are completely unchecked, so far,
+## except that they were sufficient to decide for the intermediate
+## transitions.
+##
+## Now, the exponents of the resulting config need to be checked, at least.
+##
+## Where the PA had expected an exponent difference, we now accept
+## any constant delta... which must be exactly the expected value.
+## If we have stuck exponents elsewhere, we need to assert, that the
+## output still is literally equal, as otherwise we cannot use our
+## PA-CTR repetition machinery.
+##
+## Hence, a PA-CTR has to distinguish those Vars with expected deltas
+## from those which are just there to make it match more generally.
+## At the expected deltas... exactly those deltas must result,
+## and elsewhere all exponent deltas must be exactly zero.
+##
+## For PPA-CTRs it is a bit more difficult.
+##
+## Since we must work towards a certain PA-CTR (which we fixed, already),
+## we must match the resulting tape against the input condition of that
+## PA-CTR (the PAIC).
+## - Where the PAIC has a WC, anything is ok.
+## - Where the PAIC has no var, we must match literally, we must neither
+## provide a WC nor something with a var in it.
+## - Where the PAIC has a var (in an exponent), we may provide something
+## with vars in it (those vars are from a different universe),
+## as long as the minimum value requirement is guaranteed.
+## - Where the PAIC has a var (in an exponent), we may not provide a WC
+## (FFS: that may change in future).
+##
+## When the PAIC does match, we are done. Since this planned PPA-CTR
+## will not be iterated, any var expr is ok, elsewhere.
+##
+## Unifying these different demands, we have:
+## - PA: all exponents must have exactly constant delta 0, except for
+## some explicitly named places, where another exact constant is needed.
+## - PPA: the tape must match exactly the input condition of a certain CTR.
+##
+## Hence we need for checking correct output condition:
+## - typ
+## - itape[] with WCs and vars (setup before running)
+## - otape[] with WCs and vars (generated by running)
+## - nv #vars
+## - vxtx[] where the vars are located in itape[] exponents
+## - txadd[] or similar (e.g. vxadd[]) for PA
+## - pactrx for PPA
+##
+## In order to embed the complete setup for running with vars and WCs,
+## we also need:
+## - ochx first history index (inclusive)
+## - nchx last history index (inclusive)
+## - wcL #(left WCs)
+## - wcR #(right WCs)
+## - show show level
+##
+## In order to embed the definition of the new CTR, we also need the
+## more arguments for ctr_NEW():
+## - cinfo[] (generated by running)
+## - vxEmin[] could be deduced from itape[] ? FFS
+## - def whether we shall define a new CTR (in case of success)
+##
+## That completes what we need to do a complete Vtry(), for PA or PPA.
+## We stuff (most of) the data into an array vinfo[], and pass that along,
+## thus avoiding excessively long argument lists.
+## What is in a vinfo[]:
+## - typ
+## - ochx first history index (inclusive)
+## - nchx last history index (inclusive)
+## - nv #vars
+## - vxtx[] where the vars are located in itape[] exponents
+## - txadd[] or similar (e.g. vxadd[]) for PA
+## - pactrx for PPA
+## - wcL #(left WCs)
+## - wcR #(right WCs)
+##
+## What is NOT in a vinfo[]:
+## - show show level
+## - def whether we shall define a new CTR (in case of success)
+## - itape[] with WCs and vars (setup before running)
+## - otape[] with WCs and vars (generated by running)
+## - cinfo[]
+##
+## What is needed for Vbest?
+## - where to place vars forcibly
+## - what deltas are needed at forced vars (PA)
+## - min permissable var offsets
+## - where are WCs meaningful?
+## - a strategy to follow
+##
+## These are the elements of Vbest optimization:
+## - never omit necessary vars
+## - adding WCs from both sides, as much as is possible
+## - minimize a single vars offset
+## - try adding more vars
+
+function ctr_tape_is_IC(tape,ctrx ,vx,k,vtx)
+{
+ ## Whether the tape[] matches the specified CTR's input pattern.
+ ## The tape may contain WCs, and also Vars in exponents.
+ ## (such vars are independant of those in the Itape condition)
+ ## FFS: like ctr_Match(tape,ctrx), but with Vars in tape,
+ ## (which are independant of the Vars in the CTR).
+
+ if( ! (ctrx in ctrxVcnt) ) {
+ pr("\n*** Internal error: ctr_tape_is_IC: ctrx=" ctrx); nl()
+ do_exit(1)
+ return 0
+ }
+
+ #pr("/ match sta\n")
+ if( tape["sta" ] != ctrxItape[ctrx,"sta" ] ) return 0
+ #pr("/ match 0\n")
+ if( tape[ 0 ] != ctrxItape[ctrx, 0 ] ) return 0
+ #pr("/ match T 1\n")
+ if( tape["T", 1] != ctrxItape[ctrx,"T", 1] ) return 0
+ #pr("/ match T-1\n")
+ if( tape["T",-1] != ctrxItape[ctrx,"T",-1] ) return 0
+
+ ## Build a temporary array for the (full) tape indexes that may differ.
+ ## These are just the exponents replaced by Vars.
+ for( vx=1 ; vx <= ctrxVcnt[ctrx] ; ++vx ) {
+ vtx["E",ctrxVtapx[ctrx,vx]] = ctrxVminE[ctrx,vx]
+ }
+
+ for( k in tape ) {
+ if( k == "0" ) continue # see above (FFS)
+ #pr("/ match " k "\n")
+ if( k in vtx ) { # the input pattern has a var, here
+ #pr("/ match " tape[k] " >= " vtx[k] "\n")
+ if( tape[k] == "*" ) { # known to be at least 0 (FFS: 1?)
+ ##FFS: symbol will not match, anyhow?
+ if( bnGT0(vtx[k]) ) { # CTR demands val > 0
+ #pr("/ match FAIL\n")
+ return 0 # sorry
+ }
+ }else if( ! bv_GEconst(tape[k], vtx[k]) ) {
+ #pr("/ match FAIL\n")
+ return 0 # sorry
+ }
+ }else { # non-Var pos: WC or equal
+ ## FFS: check the ctrxItape[] elem to be there
+ ## NB: at non-var places wildcards may appear in the CTR
+ #pr("/ match " tape[k] " == " ctrxItape[ctrx,k] " @ " k "\n")
+ if( "*" == ctrxItape[ctrx,k] ) {
+ ## FFS: Even a wildcard exponent currently must be > 0
+ if( "*" == tape[k] ) continue # "*" matches "*" ok
+ if( substr(k,1,2) == ("E" SUBSEP) ) {
+ if( ! bv_GT0(tape[k]) ) {
+ #pr("/ match FAIL\n")
+ return 0 # sorry
+ }
+ }
+ continue
+ }
+ if( tape[k] != ctrxItape[ctrx,k] ) {
+ #pr("/ match FAIL\n")
+ return 0 # sorry
+ }
+ }
+ }
+ return 1 # no significant difference seen: matches
+}
+
+function ctr_termOK(vinfo,itape,otape ,nx,vx,tx,e1,e2,d,w,txok,lo,hi)
+{
+ if( "PA" == vinfo["typ"] ) {
+ ## Check all exponents. At the Var places they must have delta
+ ## indicated in vinfo["txadd",tx], and otherwise they must
+ ## be literally equal.
+ nv = vinfo["nv"]
+ txok[0] = 1 ## arrify & prepare skipping
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ tx = vinfo["vxtx",vx]
+ e1 = itape["E",tx]
+ e2 = otape["E",tx]
+ d = bv_sub(e2, e1) ## real delta in tape
+ if( ! bv_IsConst(d) ) {
+ return 0 ## some Var moved!
+ }
+ w = bn0(vinfo["txadd",tx]) ## wanted (needed) delta
+ if( ("" d) != w ) { ## FFS: should not happen
+ pr("\n*** Warning: unexpected delta " d " vs " w); nl()
+ do_exit(1) ## FFS
+ return 0 ## not expected delta
+ }
+ txok[tx] = 1
+ }
+ lo = itape["T",-1]
+ hi = itape["T", 1]
+ for( tx=lo ; tx <= hi ; ++tx ) {
+ if( tx in txok ) continue
+ e1 = "" itape["E",tx]
+ e2 = "" otape["E",tx]
+ ## At this non-var place they must be strictly equal.
+ if( e1 != e2 ) {
+ return 0 ## non-var expo changed
+ }
+ }
+ return 1 ## hurra
+ }else if( "PPA" == vinfo["typ"] ) {
+ return ctr_tape_is_IC(otape, vinfo["pactrx"])
+ }else {
+ pr("\n*** Internal error: ctr_termOK: [typ] " vinfo["typ"]); nl()
+ do_exit(1)
+ return 0
+ }
+}
+
+function ctr_gVsetvars(vinfo,tape ,nv,vx,tx)
+{
+ nv = vinfo["nv"]
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ tx = vinfo["vxtx",vx]
+ tape["E",tx] = vinfo["vxEmin",vx] "+" bv_Vname(vx)
+ }
+}
+
+function ch_gV_try(vinfo,show,def ,k,ochx,nchx,typ,nv,ctrx,itape,otape,cinfo)
+{ # --> succ
+ ## generalized ch_Vtry
+ ## Here, we mainly
+ ## - set up a fresh (local) tape[] and cinfo[],
+ ## - call ch_Vrun()
+ ## - check that the result is ok
+ ## - eventually define a new CTR
+
+ if( 1 ) {
+ cinfo[0] = 0 # arrify (portability work around)
+ itape[0] = 0 # arrify (portability work around)
+ otape[0] = 0 # arrify (portability work around)
+ for( k in cinfo ) delete cinfo[k]
+ for( k in itape ) delete itape[k]
+ for( k in otape ) delete otape[k]
+ }
+
+ ochx = vinfo["ochx"]
+ nchx = vinfo["nchx"]
+ typ = vinfo["typ" ]
+ nv = vinfo["nv" ]
+
+ cinfo[0] = 0 ## arrify
+ cinfo_Init(cinfo)
+
+ itape[0] = 0 ## arrify
+ ch_ToTape(ochx, itape)
+
+ ch_WCpatch(itape, vinfo["wcL"], vinfo["wcR"])
+ ctr_gVsetvars(vinfo, itape)
+
+ for( k in itape ) otape[k] = itape[k] ## otape[] = itape[]
+
+ if( show>0 ) {
+ pr(sprintf(">> Try to prove a %s-CTR with %s Vars...\n", typ, nv))
+ tj_PrConf(otape, cinfo)
+ }
+
+ if( ! ch_Vrun(ochx, nchx, otape, cinfo, show) ) {
+ return 0
+ }
+
+ ## Check result (mostly exponents) ...
+ if( ! ctr_termOK(vinfo, itape, otape) ) {
+ if( show>0 ) pr(sprintf("<< not a %s step, sorry\n", typ))
+ return 0
+ }
+
+ if( show>0 ) {
+ if( !def ) pr("<< Success!\n")
+ if( show>1 ) {
+ for( k in cinfo ) {
+ pr(sprintf("cinfo[%8s] = ", k) cinfo[k] "\n")
+ }
+ }
+ }
+
+ if( def ) {
+ ctrx = ctr_gNEW(itape, otape, cinfo, vinfo)
+ #ctr_PA_check(ctrx,txadd)
+ if( show>0 ) {
+ pr(sprintf("<< Success! ==> defined new CTR %s (%s)\n",ctrx,typ))
+ }
+ ch_SetFirst(nchx) # FFS
+ }
+ return 1 # yes, could completely run that CTR
+}
+
+function ch_gVtry(vinfo,show)
+{
+ return ch_gV_try(vinfo,show,0)
+}
+
+function ch_gVdef(vinfo,show)
+{
+ return ch_gV_try(vinfo,show,1)
+}
+
+function ctr_gV_init(vinfo,ochx,nchx) ## vinfo-CTOR
+{
+ #vinfo["typ" ] was used to arrify
+ vinfo["ochx" ] = ochx
+ vinfo["nchx" ] = nchx
+ vinfo["nv" ] = 0
+ vinfo["pactrx"] = 0
+ vinfo["wcL" ] = 0
+ vinfo["wcR" ] = 0
+}
+
+function ctr_gVset_newvar(vinfo,tx ,ochx,vx)
+{ ## --> new vx
+ ## Invent a new var at index tx.
+ ## The offset is taken from the the saved history entry "ochx".
+ ochx = vinfo["ochx"]
+
+ vx = ++vinfo["nv"] ## allocate 1 more variable
+ vinfo["vxtx", vx] = tx
+ vinfo["vxEmin",vx] = chTap[ochx,"E",tx]
+ vinfo["vxMin" ,vx] = 1 ## absolute minimum FFS: 0-reduc
+ return vx
+}
+
+function ctr_gVset_seradd(vinfo,txser,txadd ,i,tx,d,m,vx)
+{
+ ## This is typically for a PA-CTR.
+ ## We use txser[] (forcing a serialisation) and txadd[] to build up
+ ## an initial set of vars.
+ ## We also want a minimum of 1 PA-step, which may directly imply
+ ## a minimum var offset, when the txadd[] is negative.
+
+ for( i=1 ; i<=txser[0] ; ++i ) {
+ tx = txser[i]
+ vx = ctr_gVset_newvar(vinfo, tx)
+ d = txadd[tx]
+ m = (bnLT0(d) ? bnInc(bnNeg(d)) : 1) ## FFS: 0-reduc
+ vinfo["vxMin" ,vx] = m ## must not become smaller
+ vinfo["txadd" ,tx] = d
+ }
+}
+
+function ctr_gVhasVat(vinfo,tx ,nv,vx) ## --> 0 | vx @ tx
+{
+ nv = vinfo["nv"]
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ if( tx == vinfo["vxtx",vx] ) {
+ return vx ## yes, found here
+ }
+ }
+ return 0 ## no such var found
+}
+
+function ctr_gVhasWCat(vinfo,tx)
+{
+ if( tx > 0 ) return tx <= vinfo["wcR"]
+ if( tx < 0 ) return tx >= 0 - vinfo["wcL"]
+ return 0
+}
+
+##-----------------------------------------------------------------------------
+
+function ctr_gVopt_WC(vinfo,show ,ochx,tL,tR,t,x)
+{
+ ## The current values of "wcL" and "wcR" are considered accepted.
+ ## We try to increase them both, independantly.
+ ochx = vinfo["ochx"]
+ tL = 0 - chTap[ochx,"T",-1] # |top index left side|
+ tR = chTap[ochx,"T", 1] # |top index right side|
+
+ if( show>0 ) { pr("try R WC up to " tR); nl() }
+ t = vinfo["wcR"]
+ while( t < tR ) {
+ x = t+1 # the critical tape index
+ if( ctr_gVhasVat(vinfo, x) ) break
+ vinfo["wcR"] = t+1 # try that one
+ if( ! ch_gVtry(vinfo,show) ) {
+ vinfo["wcR"] = t # sorry... set back to trusted value
+ break
+ }
+ t += 1 # accept!
+ }
+
+ if( show>0 ) { pr("try L WC up to " tL); nl() }
+ t = vinfo["wcL"]
+ while( t < tL ) {
+ x = 0 - (t+1) # the critical tape index
+ if( ctr_gVhasVat(vinfo, x) ) break
+ vinfo["wcL"] = t+1 # try that one
+ if( ! ch_gVtry(vinfo,show) ) {
+ vinfo["wcL"] = t # sorry... set back to trusted value
+ break
+ }
+ t += 1 # accept!
+ }
+ if( show>0 ) pr("WC done\n")
+}
+
+function ctr_gVopt_VX(vinfo,vx,show ,lo,hi,d,mid)
+{
+ ## We consider this var "vx" to have been verified in the current setup.
+ ## I.e. its current "vxEmin" is ok.
+ ## We try to reduce it down to the permissable minimum.
+ ## FFS: wer genau reduziert "show"?
+
+ hi = vinfo["vxEmin",vx] ## accepted
+ lo = vinfo["vxMin" ,vx] ## may force a smallest value
+ if( "" == lo ) { ## not yet any forced Min
+ lo = 1 ## FFS: 0-reduc
+ }
+ if( ! bnLT(lo,hi) ) {
+ return ## no room to do anything
+ }
+
+ vinfo["vxEmin",vx] = lo ## try the minimum, optimistically
+ if( ch_gVtry(vinfo, show) ) { ## worked! just accept that
+ return
+ }
+
+ lo = bnInc(lo) ## smallest left candidate, >0
+ while( bnLT(lo,hi) ) { ## still room to choose
+ ## Value hi is accepted already, lo is not.
+ ## Basically we bisect from below.
+ ## But occasionally the range is huge, while we expect the
+ ## result to most often be rather tiny.
+ d = length(hi) - length(lo)
+ if( d >= 3 ) { ## 10*lo < hi/10
+ mid = lo "0" ## try 10*lo
+ }else {
+ d = bnSub(hi, lo) ## so many choices
+ #d = bnDiv(d, 2) ## down rounded half of choices
+ d = bnDiv(d, 3) ## down rounded, lower part preferred
+ mid = bnAdd(lo, d) ## candidate to test
+ }
+ vinfo["vxEmin",vx] = mid ## try candidate
+ if( ch_gVtry(vinfo,show) ) { ## worked!
+ hi = mid ## so we just accept it
+ }else { ## mid did not work
+ lo = bnInc(mid) ## so we must grow larger
+ }
+ }
+ vinfo["vxEmin",vx] = hi ## restore (accept) what worked
+}
+
+function ctr_gVopt_vars(vinfo,show ,nv,vx)
+{
+ ## Just optimize the existing vars.
+ nv = vinfo["nv"]
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ ctr_gVopt_VX(vinfo, vx, show)
+ }
+}
+
+function ctr_gVopt_tryVat(vinfo,tx,show ,vx,k,nvinfo)
+{ ## --> 0 | new vx
+ ## We shall try an optional var at tx.
+ ## If it is allowed, and does run, we also optimize it.
+
+ if( ctr_gVhasVat( vinfo,tx) ) return 0
+ if( ctr_gVhasWCat(vinfo,tx) ) return 0
+
+ for( k in vinfo ) nvinfo[k] = vinfo[k] ## temp copy
+
+ vx = ctr_gVset_newvar(nvinfo,tx)
+ if( ! ch_gVtry(nvinfo, show) ) {
+ return 0 ## leave vinfo[] unchanged
+ }
+
+ ## Worked! We accept and copy the new vinfo[] back...
+ for( k in nvinfo ) vinfo[k] = nvinfo[k]
+
+ ctr_gVopt_VX(vinfo,vx,show)
+ return vx ## that one invented
+}
+
+function ctr_gVopt_moreV(vinfo,show ,ochx,lo,hi,tx)
+{
+ ## We try to add more vars... as many as possible
+
+ ochx = vinfo["ochx"]
+ lo = chTap[ochx,"T",-1]
+ hi = chTap[ochx,"T", 1]
+ for( tx=lo ; tx<=hi ; ++tx ) {
+ if( 0 == tx ) continue
+ ctr_gVopt_tryVat(vinfo,tx,show)
+ }
+}
+
+##-----------------------------------------------------------------------------
+
+function ch_CTR_rest(vinfo,show,subshow) ## --> succ
+{
+ ctr_gVopt_vars( vinfo, subshow)
+ ctr_gVopt_WC( vinfo, subshow)
+ ctr_gVopt_moreV(vinfo, subshow)
+
+ ## Now we finally go to define it (and show it along that try)...
+
+ if( ! ch_gVdef(vinfo,show) ) {
+ pr("\n*** ch_CTR_rest: internal error\n")
+ do_exit(1)
+ return 0
+ }
+ return 1
+}
+
+function ch_PA_best(ochx,nchx,show ,txadd,txser,vinfo,subshow)
+{ ## --> succ
+ subshow = show-1
+
+ txadd[0] = 0 ## arrify
+ txser[0] = 0 ## arrify
+ ch_Deltas(ochx,nchx, txadd, txser)
+
+ vinfo["typ"] = "PA" ## arrify & type it
+ ctr_gV_init(vinfo, ochx, nchx)
+ ctr_gVset_seradd(vinfo, txser, txadd)
+
+ if( ! ch_gVtry(vinfo,subshow) ) {
+ return 0 ## sorry
+ }
+ ## We already have a proof, and will succeed!
+
+ return ch_CTR_rest(vinfo,show,subshow)
+}
+
+function ch_PPA_best(ochx,nchx,pactrx,show ,vinfo) ## --> succ
+{
+ subshow = show-1
+
+ vinfo["typ"] = "PPA" ## arrify & type it
+ ctr_gV_init(vinfo, ochx, nchx)
+ vinfo["pactrx"] = pactrx
+
+ if( ! ch_gVtry(vinfo,subshow) ) { ## huh?
+ pr("\n*** ch_PA_best: internal error: cannot redo plain\n")
+ ch_gVtry(vinfo, show+2)
+ do_exit(1)
+ return 0 ## sorry
+ }
+ ## We already have a proof, and will succeed!
+
+ return ch_CTR_rest(vinfo,show,subshow)
+}
+
+##-----------------------------------------------------------------------------
+
+function ch_isPAstep(itape,tape,txadd ,lo,hi,x,d,e1,e2)
+{
+ #FFS: hat nix mit "ch" zu tun!
+ ## Both tapes are with Vars, and constructed along a saved history.
+ ## We know already, that the non-exponent part does match exactly.
+ ## We walk through the tape, and compute the exponent differences.
+ ## Where we had exponents, it may be constant != 0,
+ ## elsewhere it must be strictly zero.
+ ## Both tapes may contain wildcards at non-var positions.
+ lo = tape["T",-1]
+ hi = tape["T", 1]
+ for( x=lo ; x <= hi ; ++x ) {
+ if( x == 0 ) continue
+ e1 = "" itape["E",x]
+ e2 = "" tape["E",x]
+ #pr("[" x "] " e1 " -> " e2 "\n")
+ if( ! (x in txadd) ) { # non-var place
+ ## When the exponents are strictly equal, that is always ok,
+ ## even when we have wildcards.
+ ## When they are not strictly equal, the difference cannot
+ ## be zero, and we fail. FFS
+ if( e1 != e2 ) {
+ return 0 # non-var expo changed
+ }
+ }else { # var place: we expect a delta
+ ## Now we must not see any wildcards.
+ d = bv_sub(e2, e1)
+ #pr("[" x "] " e1 " -> " e2 ": d = " d "\n")
+ if( ! bv_IsConst(d) ) {
+ return 0 # var moved
+ }
+ ## Check against the suspected delta in txadd:
+ ## This should never fail!
+ if( ! bv_EQ(d, txadd[x]) ) {
+ pr( "\n*** Internal error: ch_isPAstep: [" x "] " \
+ d " != " txadd[x] "\n" )
+ do_exit(1)
+ return 0 # huh?
+ }
+ }
+ }
+ return 1 # hurra
+}
+
+function ch_Vtry(ochx,nchx,txadd,nv,vxtx,vxE,wcL,wcR,def,show \
+ ,vtape,tape,cinfo,tx,vx,k,ctrx)
+{
+ ## Re-Run the indicated part of the saved history,
+ ## but at the variable places replace the exponents,
+ ## and stick vars to them. For PA, only!
+ ##
+ ## ochx..nchx History index of starting & ending configuration
+ ## Both are inclusive, and for a PA nearly identical
+ ## txadd[] maps tape index to the additive delta of the PA
+ ## nv #(variables) for the CTR (1..nv = var index)
+ ## vxtx[] maps var index to tape index (without the "E")
+ ## xvE[] maps var index to minimum exponent (offset)
+ ## wcL so many left wild cards
+ ## wcR so many right wild cards
+ ## def whether to define as a new PA-CTR
+ ## show verbose level
+
+ cinfo[0] = 0 # arrify (portability work around)
+ vtape[0] = 0 # arrify (portability work around)
+ for( k in cinfo ) delete cinfo[k]
+ for( k in vtape ) delete vtape[k]
+
+ cinfo[0] = 0 # arrify
+ cinfo_Init(cinfo)
+
+ tape[0] = 0 # arrify
+ ch_ToTape(ochx,tape)
+
+ ch_WCpatch(tape, wcL, wcR)
+ ch_Vrepl( tape, nv, vxtx, vxE)
+ ch_Vattach(tape, nv, vxtx)
+ for( k in tape ) vtape[k] = tape[k] # input, with vars attached
+
+ if( show>0 ) {
+ pr(">> Try to prove a PA-CTR with " nv " Vars...\n")
+ tj_PrConf(tape, cinfo)
+ }
+ if( ! ch_Vrun(ochx, nchx, tape, cinfo, show) ) {
+ return 0
+ }
+ # Check exponents...
+ if( ! ch_isPAstep(vtape, tape, txadd) ) {
+ if( show>0 ) pr("<< not a PA step, sorry\n")
+ return 0
+ }
+ if( show>0 ) {
+ if( !def ) pr("<< Success!\n")
+ if( show>1 ) {
+ for( k in cinfo ) {
+ pr(sprintf("cinfo[%8s] = ", k) cinfo[k] "\n")
+ }
+ }
+ }
+ if( def ) {
+ ctrx = ctr_PA_new(vtape, tape, cinfo, nv, vxtx, vxE)
+ ctr_PA_check(ctrx,txadd)
+ if( show>0 ) pr("<< Success! ==> defined new CTR " ctrx "\n")
+ ch_SetFirst(nchx) # FFS
+ }
+ return 1 # yes, could completely run a PA-CTR
+}
+
+function ch_Vbest(ochx,nchx,txadd,txser,show \
+ ,nv,tx,vx,vxtx,vxE,d,vxMin,lo,hi,mid,wcL,wcR,wcok)
+{
+ ## Use the suspected deltas to build a variable map,
+ ## and extract the original exponents at the variable places
+ ## txser[] is just use to serialize the tx in txadd[].
+ ## (we want the output to be deterministic)
+ ## - nv #(vars)
+ ## - vxtx[vx] tape index where V(vx) is to be attached
+ ## - vxE[vx] exponents that should work
+ nv = 0
+ for( d=1 ; d<=txser[0] ; ++d ) {
+ tx = txser[d]
+ vx = ++nv
+ vxtx[vx] = tx
+ vxE[ vx] = chTap[ochx,"E",tx]
+ }
+ if( ! ch_Vtry(ochx,nchx,txadd,nv,vxtx,vxE,0,0,0,show-1) ) {
+ return 0
+ }
+ ## Ha, we already have done a proof!
+ ## We are going to optimize the vxE[] vector, and finally use
+ ## the best one.
+ ## Determine a minimum value that is allowed/possible
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ d = txadd[vxtx[vx]]
+ if( bnLT0(d) ) {
+ vxMin[vx] = bnInc(bnNeg(d)) #FFS: 0-reduc
+ }else {
+ vxMin[vx] = 1 #FFS: 0-reduc
+ }
+ }
+ ## Now we are going to independantly reduce the var-offsets:
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ hi = vxE[ vx] # accepted
+ lo = vxMin[vx] # smallest acceptable
+ if( ! bnLT(lo,hi) ) {
+ continue # no room to do anything
+ }
+ vxE[vx] = lo # try the minimum, optimistically
+ if( ch_Vtry(ochx,nchx,txadd,nv,vxtx,vxE,0,0,0,show-1) ) {
+ # worked... ready with minimizing vx
+ continue
+ }
+ lo = bnInc(lo) # smallest left candidate
+ while( bnLT(lo,hi) ) { # still room to choose
+ ## hi is checked already, lo not
+ ## We bisect from below...
+ d = bnSub(hi, lo) # so many choices
+ d = bnDiv(d, 2) # down rounded half of choices
+ mid = bnAdd(lo, d) # candidate to test
+ vxE[vx] = mid # try candidate
+ if( ch_Vtry(ochx,nchx,txadd,nv,vxtx,vxE,0,0,0,show-1) ) { # ok
+ hi = mid # so we just accept it
+ }else { # mid did not work
+ lo = bnInc(mid) # so we must grow larger
+ }
+ }
+ vxE[vx] = hi # restore what worked
+ }
+ ## Ok, we know our minimum exponents. That must not fail, later.
+ ## Now we try to detect wildcards, i.e. parts of the config,
+ ## which are completely untouched and therefore irrelevant.
+ ## We could find that by experiments, but we also can try to compute
+ ## it over the config history:
+ wcL = ch_WildL(ochx,nchx)
+ wcR = ch_WildR(ochx,nchx)
+ #printf("wcL/wcR = %s/%s\n", wcL, wcR)
+ ## Verify that it works... we can have erred by one,
+ ## since a cell may well change by uniting it with the neighbour.
+ ## Also, a variable place must not be a wildcard.
+ if( wcL || wcR ) {
+ for( vx=1 ; vx <= nv ; ++vx ) {
+ tx = vxtx[vx]
+ if( tx > 0 ) {
+ if( tx <= wcR ) wcR = tx - 1
+ }else if( tx < 0 ) {
+ if( (-tx) <= wcL ) wcL = (-tx) - 1
+ }
+ }
+ }
+ wcok = ( !wcL && !wcR )
+ if( ! wcok ) {
+ wcok = ch_Vtry(ochx,nchx,txadd,nv,vxtx,vxE,wcL,wcR,0,show-1)
+ }
+ if( !wcok && wcL ) {
+ if( wcok = ch_Vtry(ochx,nchx,txadd,nv,vxtx,vxE,wcL-1,wcR,0,show-1) ) {
+ --wcL
+ }
+ }
+ if( !wcok && wcR ) {
+ if( wcok = ch_Vtry(ochx,nchx,txadd,nv,vxtx,vxE,wcL,wcR-1,0,show-1) ) {
+ --wcR
+ }
+ }
+ if( !wcok && wcR && wcL ) {
+ if( wcok = ch_Vtry(ochx,nchx,txadd,nv,vxtx,vxE,wcL-1,wcR-1,0,show-1) ) {
+ --wcL ; --wcR
+ }
+ }
+ if( wcL || wcR ) {
+ if( ! ch_Vtry(ochx,nchx,txadd,nv,vxtx,vxE,wcL,wcR,0,show-1) ) {
+ pr("\n*** ch_Vbest: internal error: WC\n")
+ ch_Vtry(ochx,nchx,txadd,nv,vxtx,vxE,wcL,wcR,0,1+show)
+ do_exit(1)
+ wcL = 0
+ wcR = 0
+ }
+ }
+ #FFS
+ ## Going to define and show it:
+ if( ! ch_Vtry(ochx,nchx,txadd,nv,vxtx,vxE,wcL,wcR,1,show) ) {
+ pr("\n*** ch_Vbest: internal error\n")
+ do_exit(1)
+ }
+ return 1
+}
+
+function ch_TryDefNewPA(chx ,ochx,txadd,txser) # --> succ
+{
+ ochx = ch_Find(chx)
+ if( ochx ) {
+ #ch_Pr(ochx, "\tMatch")
+ bv_SHpush(0) ## no shortening during proofs
+ if( 1 ) {
+ if( ch_PA_best(ochx,chx,1) ) {
+ bv_SHpop()
+ return 1
+ }
+ }else {
+ txadd[0] = 0 ## arrify
+ txser[0] = 0 ## arrify
+ ch_Deltas(ochx, chx, txadd, txser)
+ if( ch_Vbest(ochx, chx, txadd, txser, 1) ) {
+ bv_SHpop()
+ return 1
+ }
+ }
+ bv_SHpop()
+ }
+ return 0 # no, sorry, nothing happened
+}
+
+##=============================================================================
+## Create PPA-CTRs
+
+## For a PPA to work, we generally accept most everything, except that
+## the intended PA-CTR must still match the result:
+## - we must not generate a WC where the PA-CTR does not expect one.
+## - we must not generate a too small Var where the PA-CTR has needs.
+## - We must not generate a Var, where the PA-CTR does not at all expect one.
+## ==>
+## - Except for the var places (in exponents) of the upcoming PA-CTR
+## we must not allow any vars in the result.
+## - Except for the WC places of the PA-CTR we must not generate WCs
+## - Even at the var places, all Vars must not have negative coeff,
+## and the base value must be large enough for the minimum of the PA-CTR.
+##
+## FFS: Our PA-CTRs should have more vars!
+
+function ch_TryDefNewPPA(chx,pactrx ,x,ochx)
+{ ## --> succ (defined new PPA)
+ ## We are currently at "chx" and about to execute a PA-CTR, "pactrx".
+ ## Let's try to prepare that better, the next time.
+ if( !chx || !pactrx ) {
+ return 0
+ }
+ ochx = 0
+ for( x=chx ; x > chBad ; --x ) {
+ ochx = x ## accept it
+ ctrx = 0 + chByx[x]
+ if( ctrx && ("PA" == ctrxType[ctrx]) ) {
+ break
+ }
+ }
+ if( ochx && ((ochx+1) < chx) ) {
+ ## Ok, there are at least 2 steps to unify.
+ ## FFS: Sollten wir mit einem Step zufrieden sein,
+ ## um es zu vereinheitlichen? --> das sehen wir erst spaeter.
+ #ch_ppaVbest(pactrx,ochx,chx)
+ bv_SHpush(0) ## no shortening during proofs
+ x = ch_PPA_best(ochx, chx, pactrx, 1)
+ bv_SHpop()
+ return x
+ }
+ return 0 #FFS sorry
+}
+
+##=============================================================================
+## Run top level machine
+## - push wanted machine
+## - Count/Collect lines/bassteps/baspos
+
+function tj_Run( mID,sta,tape,cinfo,lines,topsteps,brk,sk,skc,chx,ctrx, \
+ nocnf,newdef)
+{
+ tj_pr_symtab_warn()
+ nl()
+ if( doCH ) {
+ tjLstep = 20
+ tjLtpos = 20
+ }
+ mID = mIDtop
+ symPrID = mID # for that one we may print y-table
+ sta = sta_make_empty(mSspc[mID], 1, tj_1ori())
+ tape[0] = 0 # arrify
+ mt_Init(tape, mID, sta)
+ cinfo[0] = 0 # arrify
+ cinfo_Init(cinfo)
+ lines = 0
+ topsteps = 0
+
+ shv_Ini(shBase, cinfo["stpBas"])
+ bv_SHpush(tjShortH)
+
+ tj_PrConfHdr()
+ tj_PrConf(tape, cinfo)
+ ++lines
+ sk = shv_Cur(shBase, cinfo["stpBas"])
+ #printf("sk = %s (cur = %s)\n", sk, cinfo["stpBas"])
+
+ skc = 0
+ brk = 0
+ while( (lines < maxLines) && !brk ) {
+ if( tape["rst"] ) break
+ newdef = 0
+ ctrx = ctr_Search(tape)
+ if( !ctrx && doCH ) {
+ if( ch_TryDefNewPA(chx) ) {
+ newdef = 1
+ nocnf = 1 # showed proof
+ #pr("// research CTR\n")
+ ctrx = ctr_Search(tape)
+ #pr("// research CTR --> " ctrx "\n")
+ }
+ }
+ if( ctrx ) {
+ if( nocnf ) {
+ tj_PrConf(tape, cinfo); nocnf = 0 # showed config
+ }
+ if( "PA" == ctrxType[ctrx] ) { # could exec PA
+ if( chx && !newdef && (doCH >= 2) ) {
+ ## Try DEF new PPA here, except just defined a PA
+ #pr("Try DEF new PPA " chBad+1 ".." chx " " chCnt); nl()
+ if( ch_TryDefNewPPA(chx, ctrx) ) {
+ tj_PrConf(tape, cinfo); nocnf=0 # showed config
+ }
+ }
+ ctr_StepN(ctrx,tape,cinfo)
+ ch_SetBad(chx) # FFS, still the last config
+ ++topsteps
+ }else { # could exec PPA
+ ## FFS: only do it, when immediately following a PA!
+ ## That way we _connect_ PA-CTRs!
+ if( chx && chByx[chx] && ("PA" == ctrxType[chByx[chx]]) ) {
+ ctr_do_1Step(ctrx,tape,cinfo)
+ ch_SetBad(chx)
+ ++topsteps
+ }else {
+ #pr("HINT: ignoring CTR " ctrx); nl() #FFS
+ ctrx = 0 ## forget about the CTR (ignore it)
+ }
+ }
+ }
+ if( ! ctrx ) { ## still not done a CTR, here...
+ if( ! mt_DoStep(tape, cinfo) ) {
+ brk = 1
+ }else {
+ ++topsteps
+ }
+ }
+
+ #printf("shv = %s,%s, %s,%s\n", shBase[0],shBase[1],shBase[2],shBase[3])
+ sk = shv_Cur(shBase, cinfo["stpBas"])
+ #printf("sk = %s (cur = %s)\n", sk, cinfo["stpBas"])
+ #printf("shv = %s,%s, %s,%s\n", shBase[0],shBase[1],shBase[2],shBase[3])
+ if( (sk != "0") || tape["rst"] ) {
+ tj_PrConf(tape, cinfo); nocnf = 0 # showed config
+ ++lines
+ skc = 0
+ }else if( ! brk ) {
+ if( ! skc ) {
+ pr("...\n")
+ skc = 1
+ }
+ }
+
+ if( doCH && ! brk ) {
+ chx = ch_Save(tape, cinfo, ctrx)
+ #if( chx ) ch_Pr(chx, "\tSaved")
+ }
+ }
+ bv_SHpop() #FFS
+ sk = shv_Frc(shBase, cinfo["stpBas"])
+ if( sk != "0" ) {
+ tj_PrConf(tape, cinfo)
+ ++lines
+ }
+ tj_pr_symtab(mID)
+ tj_pr_summary(tape,cinfo,lines,topsteps)
+ symPrID = ""
+ #dump_ALL("at END")
+}
+
+function tj_Doit( i)
+{
+ bignum_app_signature()
+ ht_BEG()
+ #for( i=1 ; i<=tjILcnt ; ++i ) { pr("IN: " tjIL[i]); nl() }
+ tj_pr_basTT()
+ if( ht_VersPr(0) ) nl()
+ if( tj_PushThem() ) {
+ tj_Run()
+ }
+ ht_END()
+}
+
+##=============================================================================
+## Module initialization
+
+function tjSetup()
+{
+ gohalt = 1
+ maxLines = 201
+ shBase[2] = 1 # default: show everything
+ doCH = 0
+ tjLstep = 8
+ tjLtpos = 7
+ tj_SetNBS(2) # nBasSym = 2
+
+ basics_init() # we use "basics"
+ htSupp_init() # we use "htSupp"
+ mmSim_init() # we use "mmSim"
+ bignum_init() # we use "bignum"
+ #FFS
+}
+
+function tmJob_init()
+{
+ if( tmJob_inited_ ) return # did it, already
+ tmJob_inited_ = 1 # committed to do so
+ # register us globally:
+ g_Id[++g_IdCnt] = "$Id: tmJob.awk,v 1.34 2010/05/06 18:26:17 heiner Exp $"
+
+ tjSetup()
+}
+
+##-----------------------------------------------------------------------------
+## End of module "tmJob"
+##=============================================================================
+BEGIN {
+ tmJob_init()
+}
+ {
+ tj_Input()
+ }
+function do_exit(xc)
+{
+ if( xc ) {
+ g_xc = xc
+ }else {
+ g_xc = -1
+ }
+ exit(g_xc)
+}
+END {
+ if( g_xc ) exit(g_xc)
+ tj_Doit()
+}