Initial commit with Phodd's libraries and a README

Signed-off-by: Gavin Howard <gavin@yzena.com>

README.md

@@ -0,0 +1,20 @@
+# `bc` Libs
+
+This is a series of libraries for the [`bc` and `dc`][1] that I (Gavin D.
+Howard) built.
+
+These libraries are based on the [Phodd's excellent `bc` libraries][2], but they
+have been changed to use the extensions in my `bc`. I have also added some
+useful and not-so-useful `dc` scripts.
+
+## Why?
+
+Why have I built this repository?
+
+To build the [oral tradition of software engineering][3]. `bc` and `dc` are
+important historical tools, and I feel like adding libraries for my versions
+will push that forward.
+
+[1]: https://git.yzena.com/gavin/bc
+[2]: http://phodd.net/gnu-bc/
+[3]: https://www.youtube.com/watch?v=4PaWFYm0kEw

array.bc

@@ -0,0 +1,805 @@
+#!/usr/local/bin/bc -l
+
+### Array.BC - tools for managing arrays
+###            uses the undocumented pass-by-reference for arrays
+
+max_array_ = 4^8-1
+
+## Conventions used in this library
+/*
+ * All function names begin with the letter 'a'
+ * All array parameters have a double underscore at the end of the name
+   since some versions of bc are scope-confused by the undocumented pbr
+ * The last letter of function names indicates the type of array parameter(s)
+   expected. These are:
+   * r  = range; A range of array indices will be provided by the caller and
+          the array type is therefore moot.
+        * A number before the r suggests different ranges from different
+          array parameters.
+        * Range-specifying parameters always follow the array parameter.
+   * 0  = zero terminated; An element of zero indicates the element after
+          the last in the array. First element is always index zero
+          Advantage: Functions have the shortest possible parameter lists.
+          Disadvantage: Zero can't be used within this type of array.
+   * b  = before; All elements before a given sentinel value are treated as
+          the array. First element is index zero as for zero terminated.
+   * u  = upto; only used for the aprintu function; deliberately displays the
+          terminator for arrays that are otherwise zero/sentinel terminated
+   * l  = length provided; Array index zero contains the length of the rest
+          of the array, putting the first element at index one and the last
+          at the index specified by the length.
+   * (no last letter) = Array begins at index zero and the length will be
+          provided by the function caller
+*/
+##
+
+## Internal functions ##
+
+# For those functions which use start and end as parameters for a range,
+# use POSIX scope to sanitise those variables
+define asanerange_() {
+  auto t;
+  if(start>end){t=start;start=end;end=t}
+  if(start<0)start=0;
+  if(end>max_array_)end=max_array_;
+}
+define asanerange2_() {
+  auto start,end;
+  start=astart;end=aend;.=asanerange_();astart=start;aend=end
+  start=bstart;end=bend;.=asanerange_();bstart=start;bend=end
+}
+
+## Sorting ##
+
+# Sort a subsection of an array
+# . from index 'start' to index 'end' inclusive
+define asortr(*a__[],start,end) { # non-recursive run-finding mergesort
+  auto r[],b__[],i,j,k,ri,p,p2,ri_p2,m,mn,mid,subend;
+  .=asanerange_();
+  r[ri=0]=start;
+  i=start;while(i<end){
+    while(i<end&&a__[j=i]<=a__[++j])i=j
+    if(i>r[ri])r[++ri]=++i;
+    while(i<end&&a__[j=i]>a__[++j])i=j
+    if(i>(k=r[ri])){ # reverse backward runs
+      j=i;while(j>k){t=a__[j];a__[j--]=a__[k];a__[k++]=t}
+      r[++ri]=++i;
+    }
+  }
+  r[++ri]=++end;
+  for(p=1;p<ri;p=p2){
+    ri_p2=ri-(p2=p+p)
+    for(m=0;m<=ri_p2;m=mn){
+      # merge a__[r[m]..r[m+p]-1] with a__[r[m+p]..a[r[m+2p]-1]
+      i=r[m];mid=j=r[m+p];subend=r[mn=m+p2];k=0
+      while(i<mid&&j<subend)
+         if(a__[i]<a__[j]){ b__[k++]=a__[i++]
+                     }else{ b__[k++]=a__[j++]}
+      while(i<mid          )b__[k++]=a__[i++]
+      while(k)a__[--j]=b__[--k]
+    }
+    if(m<ri){
+      # roll in any outlier
+      # merge a[r[m-p-p]..a[r[m]-1] with a[r[m]..r[ri]-1]
+      i=r[m-p2];mid=j=r[m];k=0
+      while(i<mid&&j<end)
+         if(a__[i]<a__[j]){b__[k++]=a__[i++]
+                     }else{b__[k++]=a__[j++]}
+      while(i<mid       )  b__[k++]=a__[i++]
+      while(k)a__[--j]=b__[--k]
+      r[ri=m]=end
+    }
+  }
+  return 0;
+}
+
+define asortr_old(*a__[],start,end) { # plain recursive mergesort
+  auto os,i,j,k,t,mid,b__[];
+  if(end==start)return 0;
+  .=asanerange_();
+  if(end-start==1){
+    if(a__[start]>a__[end]){t=a__[start];a__[start]=a__[end];a__[end]=t}
+    return 0;
+  }
+  os=scale;scale=0;mid=(start+end)/2;scale=os
+  i=start;j=mid+1
+  if(i<mid).=asortr_old(a__[],start,mid)
+  if(j<end).=asortr_old(a__[],mid+1,end)
+  k=0
+  while(i<=mid&&j<=end)
+     if(a__[i]<a__[j]){b__[k++]=a__[i++]
+                 }else{b__[k++]=a__[j++]}
+  while(i<=mid        )b__[k++]=a__[i++]
+ #while(        j<=end)b__[k++]=a__[j++]
+  while(k>0)a__[--j]=b__[--k]
+  return 0;
+}
+
+# Sort all elements of array before zero terminator
+define asort0(*a__[]){
+  auto i;
+  for(i=0;a__[i];i++){}
+  .=asortr(a__[],0,i-1)
+}
+
+# Sort all elements of array before given terminator, x
+define asortb(*a__[],x){
+  auto i;
+  for(i=0;a__[i]!=x;i++){}
+  .=asortr(a__[],0,i-1);
+}
+
+# Sort all elements of array with length in [0]
+define asortl(*a__[]){
+  .=asortr(a__[],1,a__[0]);
+}
+
+# Sort the first n elements of an array
+define asort(*a__[],n) { .=asortr(a__[],0,n-1) }
+
+## Unique values / Run length finding ##
+
+# Store values from a in v & how many times they occur together in a run in r
+# e.g. if a=={7,8,8,9,9,6,9}, v <- {7,8,9,6,9} and r <- {1,2,2,1,1}
+define arunlengthr(*v__[],*r__[], a__[],start,end) {
+  auto vri,i,prev;
+  .=asanerange_();
+  if(end==start)return 0;
+  vri=0;prev=v__[vri]=a__[start];r__[vri]=1
+  for(i=start+1;i<=end;i++)if(a__[i]==prev){
+     .=r__[vri]++
+    } else {
+      r__[++vri]=1; prev=v__[vri]=a__[i]
+    }
+  return ++vri
+}
+define arunlength0(*v__[],*r__[], a__[]) {
+  auto vri,i,prev;
+  if(!a__[0])return 0;
+  vri=0;prev=v__[vri]=a__[0];r__[vri]=1
+  for(i=1;i<=max_array_&&a__[i];i++)if(a__[i]==prev){
+     .=r__[vri]++
+    } else {
+      r__[++vri]=1; prev=v__[vri]=a__[i]
+    }
+  .=vri++
+  r__[vri]=v__[vri]=0
+  return vri
+}
+define arunlengthb(*v__[],*r__[], a__[],x) {
+  auto vri,i,prev;
+  if(a__[0]==x)return 0;
+  vri=0;prev=v__[vri]=a__[0];r__[vri]=1
+  for(i=1;i<=max_array_&&a__[i]!=x;i++)if(a__[i]==prev){
+     .=r__[vri]++
+    } else {
+      r__[++vri]=1; prev=v__[vri]=a__[i]
+    }
+  .=vri++
+  r__[vri]=v__[vri]=x
+  return vri
+}
+define arunlengthl(*v__[],*r__[], a__[]) {
+  auto vri,i,prev;
+  if(!a__[0])return 0;
+  vri=1;prev=v__[vri]=a__[1];r__[vri]=1
+  for(i=2;i<=max_array_&&a__[i]!=x;i++)if(a__[i]==prev){
+     .=r__[vri]++
+    } else {
+      r__[++vri]=1; prev=v__[vri]=a__[i]
+    }
+  return r__[0]=v__[0]=vri
+}
+define arunlength(*v__[],*r__[], a__[], n) {
+  return arunlengthr(v__[],r__[], a__[],0,n-1)
+}
+
+# Fill v with the unique values from a. Works best on a sorted array.
+define auniqr(*v__[], a__[],start,end) {
+  auto r[];
+  return arunlengthr(v__[],r[], a__[],start,end)
+}
+define auniq0(*v__[], a__[]) {
+  auto r[];
+  return arunlength0(v__[],r[], a__[])
+}
+define auniqb(*v__[], a__[],x) {
+  auto r[];
+  return arunlengthb(v__[],r[], a__[],x)
+}
+define auniql(*v__[], a__[]) {
+  auto r[];
+  return arunlengthl(v__[],r[], a__[])
+}
+define auniq(*v__[], a__[],n) {
+  auto r[];
+  return arunlengthr(v__[],r[], a__[],0,n-1)
+}
+
+## Order reversal ##
+
+# Reverse a subsection of an array
+define areverser(*a__[],start,end) {
+  auto t;
+  .=asanerange_();
+  if(end==start)return 0
+  while(start<end){t=a__[start];a__[start++]=a__[end];a__[end--]=t}
+  return 0;
+}
+
+# Reverse order of all elements of array before given terminator, x
+define areverseb(*a__[],x){
+  auto i;
+  for(i=0;a__[i]!=x;i++){}
+  .=areverser(a__[],0,i-1);
+}
+
+# Reverse order of all elements of array before zero terminator
+define areverse0(*a__[]){
+  auto i;
+  for(i=0;a__[i];i++){}
+  .=areverser(a__[],0,i-1)
+}
+
+# Reverse order of all elements of array 1 .. a[0]
+define areversel(*a__[]) {
+  .=areverser(a__[],1,a__[0]);
+}
+
+# Reverse the first n elements of an array
+define areverse(*a__[],n) { .=areverser(a__[],0,n-1) }
+
+## Copying ##
+
+# Make a__[astart..aend] = b__[bstart..bend]
+# . Does not expand a__[] to make room if b__[bstart..bend] is too large!
+define acopy2r(*a__[],astart,aend,b__[],bstart,bend){
+  auto i,j;
+  .=asanerange2_();
+  i=astart;j=bstart;while(i<=aend&&j<=bend)a__[i++]=b__[j++]
+}
+
+# Make a__[start..end] = b__[start..end]
+define acopyr(*a__[],b__[],start,end){
+  auto i;
+  .=asanerange_();
+  for(i=start;i<=end;i++)a__[i]=b__[i]
+}
+
+# copy a zero-terminated array; a = b
+define acopy0(*a__[],b__[]) {
+  auto e,i;
+  for(i=0;i<=max_array_&&e=b__[i];i++)a__[i]=e
+  if(i<=max_array_)a__[i]=0
+}
+
+# copy an x terminated array; a = b
+define acopyb(*a__[],b__[],x) {
+  auto e,i;
+  for(i=0;i<=max_array_&&(e=b__[i])!=x;i++)a__[i]=e
+  if(i<=max_array_)a__[i]=x
+}
+
+# Copy array whose length is in element [0]
+define acopyl(*a__[],b__[]){
+  .=acopyr(a__[],b__[],0,a__[0])
+}
+
+# Copy first 'count' elements of a from b
+define acopy(*a__[],b__[],count){
+  auto i;
+  if(count<0)count=0;
+  if(count>max_array_)count=max_array_+1;
+  for(i=0;i<count;i++)a__[i]=b__[i];return 0
+}
+
+## Convert between array types ##
+
+define aconv0fromr(*a__[],  b__[],start,end){
+  auto i,j;
+  .=asanerange_();
+  i=0;for(j=start;j<=end;j++)a__[i++]=b__[j]
+  a__[i]=0
+}
+define aconvbfromr(*a__[],x,b__[],start,end){
+  auto i,j;
+  .=asanerange_();
+  i=0;for(j=start;j<=end;j++)a__[i++]=b__[j]
+  a__[i]=x
+}
+define aconvlfromr(*a__[],  b__[],start,end){
+  auto i,j;
+  .=asanerange_();
+  i=1;for(j=start;j<=end;j++)a__[i++]=b__[j]
+  a__[0]=end-start+1
+}
+
+define aconvrfrom0(*a__[],start,end,b__[]){
+  auto i,j;
+  .=asanerange_();
+  j=0;for(i=start;i<=end&&b__[j];i++)a__[i]=b__[j++]
+}
+define aconvbfrom0(*a__[],x,        b__[]){
+  auto i;
+  for(i=0;a__[i]=b__[i];i++){}
+  a__[i]=x
+}
+define aconvlfrom0(*a__[],          b__[]){
+  auto i,j;
+  i=1;j=0;while(a__[i++]=b__[j++]){}
+  a__[0]=j
+}
+
+define aconv0fromb(*a__[],          b__[],x){
+  auto i;
+  for(i=0;(a__[i]=b__[i])!=x;i++){}
+  a__[i]=0
+}
+define aconvrfromb(*a__[],start,end,b__[],x){
+  auto i,j;
+  .=asanerange_();
+  j=0;for(i=start;i<=end&&b__[j]!=x;i++)a__[i]=b__[j++]
+}
+define aconvlfromb(*a__[],          b__[],x){
+  auto i,j;
+  i=1;j=0;while((a__[i++]=b__[j++])!=x){}
+  a__[0]=j
+}
+
+define aconvrfroml(*a__[],start,end,b__[]){
+  auto i,j;
+  .=asanerange_();
+  j=1;for(i=start;i<=end&&j<=b__[0];i++)a__[i]=b__[j++]
+}
+define aconv0froml(*a__[],          b__[]){
+  auto i,j;
+  i=0;j=1;while(j<=b__[0])a__[i++]=b__[j++];
+  a__[i]=0;
+}
+define aconvbfroml(*a__[],x,        b__[]){
+  auto i,j;
+  i=0;j=1;while(j<=b__[0])a__[i++]=b__[j++];
+  a__[i]=x;
+}
+
+## Copy in a small number of elements ##
+
+define aset8(*a__[],start,a,b,c,d,e,f,g,h){
+  auto b[];
+  b[0]=a;b[1]=b;b[2]=c;b[3]=d
+  b[4]=e;b[5]=f;b[6]=g;b[7]=h
+  .=acopy2r(a__[],start,start+7,b[],0,7)
+}
+
+define aset16(*a__[],start,a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p){
+  auto b[];
+  b[0]=a;b[1]=b;b[2]=c;b[3]=d
+  b[4]=e;b[5]=f;b[6]=g;b[7]=h
+  b[8]=i;b[9]=j;b[A]=k;b[B]=l
+  b[C]=m;b[D]=n;b[E]=o;b[F]=p
+  .=acopy2r(a__[],start,start+F,b[],0,F)
+}
+
+## Append one array to another ##
+
+define aappendr(*a__[],aend,b__[],bstart,bend) {
+  auto i,j;
+  if(aend>max_array_)aend=max_array_;
+  if(bstart>bend){j=bstart;bstart=bend;bend=j}
+  if(bstart<0)bstart=0;
+  if(bend>max_array_)bend=max_array_;
+  j=bstart;
+  for(i=aend+1;i<=max_array_&&j<=bend;i++)a__[i]=b__[j++]
+}
+
+define aappend0(*a__[],b__[]) {
+  auto i,j;
+  for(i=0;i<=max_array_&&a__[i];i++){}
+  if(i>max_array_)return 0;
+  for(j=0;i<=max_array_&&b__[j];j++)a__[i++]=b__[j]
+  if(i>max_array_)return 0;
+  a__[i]=0;
+}
+
+define aappendb(*a__[],b__[],x) {
+  auto i,j;
+  for(i=0;i<=max_array_&&a__[i]!=x;i++){}
+  if(i>max_array_)return 0;
+  for(j=0;i<=max_array_&&b__[j]!=x;j++)a__[i++]=b__[j]
+  if(i>max_array_)return 0;
+  a__[i]=x;
+}
+
+define aappendl(*a__[],b__[]) {
+  .=aappendr(a__[],a__[0],b__[],1,b__[0])
+  a__[0]+=b__[0]
+}
+
+define aappend(*a__[],aend,b__[],count) {
+  return aappendr(a__[],aend,b__[],0,count-1)
+}
+
+define aappendelem0(*a__[],e) {
+  auto i;
+  for(i=0;i<=max_array_&&a__[i];i++){}
+  if(i>max_array_)return 0;
+  a__[i++]=e;
+  if(i>max_array_)return 0;
+  a__[i]=0;
+}
+
+define aappendelemb(*a__[],x,e) {
+  auto i;
+  for(i=0;i<=max_array_&&a__[i]!=x;i++){}
+  if(i>max_array_)return 0;
+  a__[i++]=e;
+  if(i>max_array_)return 0;
+  a__[i]=x;
+}
+
+define aappendelem(*a__[],aend,e) {
+  if(++aend>max_array_)return 0;
+  a__[aend]=e;return 0
+}
+
+define aappendeleml(*a__[],e) {
+  if(++a__[0]>max_array_){a__[0]=max_array_;return 0}
+  a__[a__[0]]=e
+}
+
+## part a, part b, part c
+
+define aparts3r(*a__[],astart,aend,b__[],bstart,bend,c__[],cstart,cend) {
+  auto i,j;
+  .=asanerange2_();
+  if(cstart>cend){i=cstart;cstart=cend;cend=i}
+  if(cstart<0)cstart=0;
+  if(cend>max_array_)cend=max_array_;
+  i=0;
+  if(astart>0)for(j=astart;j<aend;j++)if(i++<=max_array_)a__[i]=a__[j];
+  for(j=bstart;j<bend;j++)if(i++<=max_array_)a__[i]=b__[j];
+  for(j=cstart;j<cend;j++)if(i++<=max_array_)a__[i]=c__[j];
+}
+
+## Insertion of elements and other arrays ##
+
+define ainsertatr(*a__[],astart,aend,pos,b__[],bstart,bend){
+  auto i;
+  .=asanerange2_();
+  if(pos>aend-astart)return aappendr(a__[],aend,b__[],bstart,bend);
+  if(pos<=0){
+    .=aappendr(b__[],bend,a__[],astart,aend);
+    .=acopyr(a__[],b__[],bstart,aend-astart+bend)
+    return 0;
+  }
+  .=aparts3r(a__[],astart,astart+pos-1,b__[],bstart,bend,a__[],astart+pos,aend);
+}
+
+define ainsertat0(*a__[],pos,b__[]){
+  auto i,j,alen,blen;
+  for(i=0;i<=max_array_&&b__[i];i++){}
+  blen=i
+  for(i=0;i<=max_array_&&a__[i];i++){}
+  alen=i
+  if(pos>=alen)return aappend0(a__[],b__[]);
+  if(pos<0)pos=0;
+  for(--i;i>=pos;i--){
+    if((j=blen+i)>max_array_)continue;
+    a__[j]=a__[i]
+  }
+  for(i=0;i<blen;i++){
+    if((j=pos+i)>max_array_)continue;
+    a__[j]=b__[i]
+  }
+  a__[alen+blen]=0
+}
+define ainsertatb(*a__[],pos,b__[],x){
+  auto i,j,alen,blen;
+  for(i=0;i<=max_array_&&b__[i]!=x;i++){}
+  blen=i
+  for(i=0;i<=max_array_&&a__[i]!=x;i++){}
+  alen=i
+  if(pos>=alen)return aappend0(a__[],b__[]);
+  if(pos<0)pos=0;
+  for(--i;i>=pos;i--){
+    if((j=blen+i)>max_array_)continue;
+    a__[j]=a__[i]
+  }
+  for(i=0;i<blen;i++){
+    if((j=pos+i)>max_array_)continue;
+    a__[j]=b__[i]
+  }
+  a__[alen+blen]=x
+}
+define ainsertatl(*a__[],pos,b__[]){
+  auto i;
+  if(pos>a__[0])return aappendl(a__[],b__[]);
+  if(pos<1)pos=1;
+  for(i=a__[0];i>=pos;i--){
+    if((j=b__[0]+i)>max_array_)continue;
+    a__[j]=a__[i]
+  }
+  for(i=1;i<=b__[0];i++){
+    if((j=pos+i-1)>max_array_)continue;
+    a__[j]=b__[i]
+  }
+  a__[0]+=b__[0]
+}
+define ainsertat(*a__[],acount,pos,b__[],bcount){
+  return ainstertatr(a__[],0,acount-1,pos,b__[],0,bcount-1);
+}
+
+define ainsertelematr(*a__[],start,end,pos,e){
+  auto b[];b[0]=e
+  return ainsertatr(a__[],start,end,pos,b[],0,0);
+}
+define ainsertelemat0(*a__[],pos,e){
+  auto b[];b[0]=e;b[1]=0
+  return ainsertat0(a__[],pos,b[]);
+}
+define ainsertelematb(*a__[],x,pos,e){
+  auto b[];b[0]=e;b[1]=x
+  return ainsertatb(a__[],pos,b[],x);
+}
+define ainsertelematl(*a__[],pos,e){
+  auto b[];b[0]=1;b[1]=e
+  return ainsertatl(a__[],pos,b[]);
+}
+define ainsertelemat(*a__[],count,pos,e){
+  auto b[];b[0]=e
+  return ainstertatr(a__[],0,count-1,pos,b[],0,0);
+}
+
+## Comparison -> Lexical ##
+
+# Compare a__[astart..aend] and b__[bstart..bend]
+define acompare2r(a__[],astart,aend,b__[],bstart,bend){
+  auto a,b,i,j;
+  .=asanerange2_();
+  i=astart;j=bstart;
+  while(i<=aend&&j<=bend){
+    a=a__[i++];b=b__[j++]
+    if(a<b)return -1;
+    if(a>b)return  1;
+  }
+  if(i>aend&&j>bend)return 0;#sub arrays are equivalent
+  if(i>aend)return -1;#left array is shorter
+  return 1;#right array is shorter
+}
+
+# Compare a__[start..end] and b__[start..end]
+define acomparer(a__[],b__[],start,end){
+  auto i;
+  .=asanerange_();
+  for(i=start;i<=end;i++){
+    if(a__[i]<b__[i])return -1;
+    if(a__[i]>b__[i])return  1;
+  }
+  return 0;
+}
+
+# compare zero-terminated arrays; a <=> b
+define acompare0(a__[],b__[]) {
+  auto i;
+  for(i=0;i<=max_array_&&a__[i]&&b__[i];i++){
+    if(a__[i]<b__[i])return -1;
+    if(a__[i]>b__[i])return  1;
+  }
+  if(i>max_array_||(a__[i]==0&&b__[i]==0))return 0;
+  if(a__[i]==0)return -1;#left array is shorter
+  return 1;
+}
+
+# compare x terminated arrays; a <=> b
+define acompareb(a__[],b__[],x) {
+  auto i;
+  for(i=0;i<=max_array_&&a__[i]!=x&&b__[i]!=x;i++){
+    if(a__[i]<b__[i])return -1;
+    if(a__[i]>b__[i])return  1;
+  }
+  if(i>max_array_||(a__[i]==x&&b__[i]==x))return 0;
+  if(a__[i]==x)return -1;#left array is shorter
+  return 1;
+}
+
+# compare length-specified arrays; a <=> b
+define acomparel(a__[],b__[]) {
+  auto i;
+  for(i=1;i<=max_array_&&i<=a__[0]&&i<=b__[0];i++){
+    if(a__[i]<b__[i])return -1;
+    if(a__[i]>b__[i])return  1;
+  }
+  if(i>max_array_||(i>a__[0]&&i>b__[0]))return 0;
+  if(i>a__[0])return -1;#left array is shorter
+  return 1;
+}
+
+# Compare first 'count' elements of a and b
+define acompare(a__[],b__[],count){
+  auto i;
+  if(count<0)count=0;
+  if(count>max_array_)count=max_array_+1;
+  for(i=0;i<count;i++){
+    if(a__[i]<b__[i])return -1;
+    if(a__[i]>b__[i])return  1;
+  }
+  return 0;
+}
+
+## Comparison -> Equality ##
+
+# Check equality of a__[astart..aend] and b__[bstart..bend]
+define   aequals2r(a__[],astart,aend,b__[],bstart,bend){
+return !acompare2r(a__[],astart,aend,b__[],bstart,bend)}
+
+# Check equality of a__[start..end] and b__[start..end]
+define   aequalsr(a__[],b__[],start,end){
+return !acomparer(a__[],b__[],start,end)}
+
+# Check equality of zero-terminated arrays; a == b
+define   aequals0(a__[],b__[]){
+return !acompare0(a__[],b__[])}
+
+# Check equality of x terminated arrays; a == b
+define   aequalsb(a__[],b__[],x){
+return !acompareb(a__[],b__[],x)}
+
+# Check equality of length specified arrays; a == b
+define   aequalsl(a__[],b__[]){
+return !acomparel(a__[],b__[])}
+
+# Check equality of first 'count' elements of a and b
+define   aequals(a__[],b__[],count){
+return !acompare(a__[],b__[],count)}
+
+## Subarray matching ##
+
+# Returns position in large that first occurrence of small is found
+# . or -1 if no match is found
+define amatcharray2r(small__[],astart,aend,large__[],bstart,bend){
+  auto i,j;
+  .=asanerange2_();
+  if(aend-astart>bend-bstart)return -1;
+  i=astart;j=bstart;
+  while(j<=bend){
+    while(         j<=bend&&small__[i]!=large__[j]){.=j++}
+    if(j>bend)return -1;
+    .=i++   ;.=j++
+    while(i<=aend&&j<=bend&&small__[i]==large__[j]){.=i++;.=j++}
+    if(i>aend)return j-aend+astart-1;
+    i=astart;.=j++
+  }
+  return -1;
+}
+
+# As above but large and small are zero terminated arrays
+define amatcharray0(small__[],large__[]){
+  auto i,j;
+  for(i=0;i<=max_array_&&small__[i];i++){}
+  for(j=0;j<=max_array_&&large__[j];j++){}
+  if(i>j)return -1;
+  return amatcharray2r(small__[],0,i-1,large__[],0,j-1)
+}
+
+# As above but large and small are arrays terminated by x
+define amatcharrayb(small__[],large__[],x){
+  auto i,j;
+  for(i=0;i<=max_array_&&small__[i]!=x;i++){}
+  for(j=0;j<=max_array_&&large__[j]!=x;j++){}
+  if(i>j)return -1;
+  return amatcharray2r(small__[],0,i-1,large__[],0,j-1)  
+}
+
+# As above but large and small are arrays terminated by x
+define amatcharrayl(small__[],large__[]){
+  if(small__[0]>large__[0])return -1;
+  return amatcharray2r(small__[],1,small__[0],large__[],1,large__[0])  
+}
+
+define amatcharray(small__[],scount,large__[],lcount){
+  if(scount>lcount) return -1;
+  return amatcharray2r(small__[],0,scount-1,large__[],0,lcount-1);
+}
+
+## Single element finding ##
+
+define amatchelementr(e,a__[],start,end){
+  auto i;
+  .=asanerange_();
+  for(i=start;i<=end;i++)if(a__[i]==e)return i;
+  return -1;
+}
+
+define amatchelement0(e,a__[]){
+  auto i;
+  for(i=0;i<=max_array_&&a__[i];i++)if(a__[i]==e)return i;
+  if(i<=max_array_&&e==0)return i; 
+  return -1;
+}
+
+define amatchelementb(e,a__[],x){
+  auto i;
+  for(i=0;i<=max_array_&&a__[i]!=x;i++)if(a__[i]==e)return i;
+  if(i<=max_array_&&e==x)return i;
+  return -1;
+}
+
+define amatchelementl(e,a__[]){
+  auto i;
+  for(i=1;i<=max_array_&&i<=a__[0];i++)if(a__[i]==e)return i;
+  return -1;
+}
+
+define amatchelement(e,a__[],count){
+  return amatchelementr(e,a__[],0,count-1);
+}
+
+## Output ##
+
+aprint_wide_=1
+
+# Prints a__[start..end]
+define aprintr(*a__[],start,end) {
+  auto i;
+  if(start>end){print"{}";return 0}
+  .=asanerange_();
+  print "{";for(i=start;i<end;i++){
+    print a__[i],", ";if(!aprint_wide_)print"\n "
+  };print a__[i],"}\n"
+}
+
+# Treats a__[] as a zero terminated array; prints all elements prior to 0
+define aprint0(*a__[]) {
+  auto e,f,i;
+  print "{";
+  if(!(e=a__[i=0])){print "}\n";return 0}
+  while(++i<=max_array_&&f=a__[i]){
+    print e,", "
+    if(!aprint_wide_)print"\n "
+    e=f;
+  }
+  print e,"}\n"
+}
+
+# Treats a__[] as an 'x' terminated array; prints all elements [b]efore 'x'
+define aprintb(*a__[],x) {
+  auto e,f,i;
+  print "{";
+  if(x==(e=a__[i=0])){print "}\n";return 0}
+  while(++i<=max_array_&&(f=a__[i])!=x){
+    print e,", "
+    if(!aprint_wide_)print"\n "
+    e=f;
+  }
+  print e,"}\n"
+}
+
+# Treats a__[] as an 'x' terminated array 
+# . but prints all elements [u]pto and including 'x'
+define aprintu(*a__[],x) {
+  auto i;
+  print "{";
+  for(i=0;a__[i]!=x&&i<max_array_;i++) {
+    print a__[i],", "
+    if(!aprint_wide_)print"\n "
+  }
+  if(i==max_array_){print a__[i]}else{print x}
+  print "}\n"
+}
+
+# Treats a__[] as a length specified array and prints it
+define aprintl(*a__[]) {
+  auto e,f,i;
+  if(a__[0]==0){print"{}";return 0}
+  print "{";
+  for(i=1;i<a__[0];i++) {
+    print a__[i],", "
+    if(!aprint_wide_)print"\n "
+  }
+  print a__[i],"}\n"
+}
+
+# Prints the first 'count' elements of a__[]
+define aprint(*a__[],count) {
+ .=aprintr(a__[],0,count-1)
+}

cf.bc

@@ -0,0 +1,314 @@
+#!/usr/local/bin/bc -l
+
+### CF.BC - Continued fraction experimentation using array pass by reference
+
+# Best usage for the output functions in this library:
+#   .=output_function(params)+newline()
+# This suppresses the output and also appends a newline which is not added
+# by default.
+
+## Initialisation / Tools / Workhorses ##
+
+# Borrowed from funcs.bc
+define int(x) { auto os;os=scale;scale=0;x/=1;scale=os;return(x) }
+define abs(x) { if(x<0)return(-x)else return(x) }
+
+# Borrowed from output_formatting.bc
+define newline() { print "\n" }
+
+# sanity check for the below
+define check_cf_max_() {
+  auto maxarray;
+  maxarray = 4^8-4;cf_max=int(cf_max)
+  if(1>cf_max||cf_max>maxarray)cf_max=maxarray
+  return 0;
+}
+
+# global var; set to halt output and calculations at a certain point
+cf_max=-1; # -1 to flag to other libs that cf.bc exists before value is set
+
+# Workhorse function to prepare for and tidy up after CF generation
+define cf_tidy_(){
+  # POSIX scope; expects vars cf__[], max, p and i
+  # Tidy up the end of the CF
+  #  assumes the last element of large CFs to be invalid due to rounding
+  #    and deletes that term.
+  #  for apparently infinite (rather than just long) CFs, uses a bc trick
+  #    to signify special behaviour to the print functions.
+  #  for apparently finite CFs, checks whether the CF ends ..., n, 1]
+  #  which can be simplified further to , n+1]
+  auto j,diff,rl,maxrl,bestdiff
+  cf__[i]=0
+  if(p>max){
+    .=i--;
+    if(p<max*A^5){ # assume infinite
+      cf__[i]=0.0; # bc can tell between 0 and 0.0
+      # Identify repeating CFs and store extra information
+      #  that is used by the print functions.
+      max=i;mrl=0;
+      p=max-1;if((i=p-1)>0)for(i=i;i;i--){
+        if(cf__[i]==cf__[p]){
+          diff=p-i;rl=0
+          for(j=p-1;j>diff;j--){if(cf__[j-diff]!=cf__[j])break;.=rl++}
+          if(j<=i)rl+=diff
+          if(rl>mrl){mrl=rl;bestdiff=diff}
+        }
+      }
+      if(3*mrl>=p){cf__[++max]=p-mrl-1;cf__[++max]=bestdiff;cf__[0]+=0.0}
+      return 0;
+    }
+    cf__[i]=0; # bc can differentiate between 0 and 0.0
+  } 
+  if(i<2)return 0;
+  .=i--;if(abs(cf__[i])==1){cf__[i-1]+=cf__[i];cf__[i]=0}
+  return 0;
+}
+
+# Workhorse function for cf_new and cfn_new
+define cf_new_(near) {
+  # POSIX scope; expects array *cf__[] and var x
+  auto os, i, p, max, h, n, temp;
+  .=check_cf_max_()
+  os=scale;if(scale<scale(x))scale=scale(x)
+  if(scale<6+near){
+    print "cf";if(near)print"n";print"_new: scale is ",scale,". Poor results likely.\n"
+  }
+  max=A^int(scale/2);p=1
+  if(near)h=1/2;
+  for(i=0;i<cf_max&&p<max;i++) {
+    if(near){
+      n=int(temp=x+h);if(n>temp).=n-- # n=floor(x+.5)
+    } else {
+      n=int(x)
+    }
+    p*=1+abs(cf__[i]=n);
+    x-=cf__[i]
+    if(x==0||p==0){.=i++;break}
+    x=1/x
+  }
+  scale=os
+  return cf_tidy_();
+}
+
+## Making/unmaking continued fractions and ordinary fractions ##
+
+# Create a continued fraction representation of x in pbr array cf__[]
+define cf_new(*cf__[],x) {
+  return x+cf_new_(0);
+}
+
+# Create a continued fraction representation of x in pbr array cf__[]
+# using signed terms to guarantee largest magnitude following term
+define cfn_new(*cf__[],x) { 
+  return x+cf_new_(1);
+}
+
+# Copy a continued fraction into pbr array cfnew__[] from pbv array cf__[]
+define cf_copy(*cfnew__[],cf__[]){
+  auto e,i;
+  for(i=0;i<=cf_max&&e=cf__[i];i++)cfnew__[i]=e
+  if(i<=cf_max)cfnew__[i]=cf__[i] # might be 0.0
+  e=cf__[0];if(int(e)==e&&scale(e)){
+    # copy extra info
+    if(++i<=cf_max)cfnew__[i]=cf__[i]
+    if(++i<=cf_max)cfnew__[i]=cf__[i]
+  }
+  return (i>cf_max);
+}
+
+# Convert pbv array cf__[] into a number
+define cf_value(cf__[]) {
+  auto n, d, temp, i;
+  .=check_cf_max_();
+  if(cf__[1]==0)return cf__[0];
+  for(i=1;i<cf_max&&cf__[i];i++){}
+  n=cf__[--i];d=1
+  for(i--;i>=0;i--){temp=d;d=n;n=temp+cf__[i]*d}
+  return(n/d);
+}
+
+# Convert pbv array cf__[] into a new cf of the other type
+# . with respect to nearness.
+define cf_toggle(*cf__[],cf2__[]){
+  auto os,p,i,x,zero,sign,near;
+  sign=1;if(cf2__[0]<0)sign=-1
+  near=0
+  for(i=1;x=cf2__[i];i++){
+    p*=1+abs(x)
+    if(x*sign<0)near=1
+  }
+  zero=x;
+  os=scale;scale=0
+   for(i=2;p;i++)p/=A
+   if(i<os)i=os
+  scale=i
+   x=cf_value(cf2__[])
+   .=cf_new_(1-near)
+  scale=os
+  for(i=0;cf__[i];i++){}
+  cf__[i]=zero;
+  return 0
+}
+
+# Return the value of the specified convergent of a CF
+define cf_get_convergent(cf__[],c){
+  auto ocm,v;
+  if(c==0)return cf_value(cf__[]);
+  .=check_cf_max_();ocm=cf_max
+  cf_max=c;v=cf_value(cf[])
+  cf_max=ocm;return v;
+}
+
+# Return the value of the specified convergent of the CF of x
+define get_convergent(x,c){
+  auto cf[];
+  .=cf_new(cf[],x)
+  return cf_get_convergent(cf[],c)
+}
+
+# Create denominator, numerator and optional intpart in
+# . first three elements of pbr array f__[]
+# . from CF in pbv array cf[]()
+# NB: returned elements are in reverse of the expected order!
+define frac_from_cf(*f__[],cf__[],improper) {
+  auto n,d,i,temp;
+  .=check_cf_max_();
+  improper=!!improper;
+  if(cf__[0]==(i=int(cf__[0])))cf__[0]=i
+  if(cf__[1]==0){f__[0]=1;f__[1]=0;return f__[2]=cf__[0]}
+  for(i=1;i<cf_max&&cf__[i];i++){}
+  n=cf__[--i];d=1
+  for(i--;i>=!improper;i--){
+    temp=n;n=d;d=temp # reciprocal = swap numerator and denominator
+    n+=cf__[i]*d
+  }
+  temp=0
+  if(!improper){temp=n;n=d;d=temp}#correct for having stopped early
+  if(d<0){n=-n;d=-d} # denominator always +ve
+  if(!improper&&cf__[0]!=0){
+    temp=cf__[0]
+  }
+  f__[0]=d;f__[1]=n;f__[2]=temp
+  return temp+n/d;
+}
+
+# Upscale an allegedly rational number to the current scale
+define upscale_rational(x) {
+  auto os,f[],cf[];
+  if(scale<=scale(x))return x
+  os=scale;scale=scale(x)
+  .=cf_new(cf[],x);          # Sneaky trick to upscale (literally)
+  .=frac_from_cf(f[],cf[],1) # any rational value of x
+  scale=os
+  x=f[1]/f[0] # x is now (hopefully) double accuracy!
+  return x;
+}
+
+## Output ##
+
+# Set to 1 to truncate quadratic surd output before repeat occurs
+cf_shortsurd_=0
+
+# Output pbv array cf__[] formatted as a CF
+define cf_print(cf__[]) {
+  auto i,sign,surd,sli,sri;
+  .=check_cf_max_();
+  sign=1;if(cf__[0]<0||(cf__[0]==0&&cf__[1]<0)){sign=-1;print "-"}
+  # Check for surd flag and find lead in and repeat period
+  surd=0;sri=int(i=cf__[0]);if(scale(i)&&sri==i){
+    for(i=1;cf__[i];i++){}
+    sli=cf__[++i]
+    surd=cf__[++i]
+  }
+  print "[", sign*sri;
+  if(cf__[1]==0){print "]";return 0}
+  print "; ";
+  if(surd)if(sli){sri=sli-1}else{sri=surd-1}
+  for(i=2;i<=cf_max&&cf__[i];i++){
+    print sign*cf__[i-1]
+    if(!surd||sri){
+      print ", ";if(surd).=sri--
+    } else {
+      if(sli){
+        sli=0
+      }else if(cf_shortsurd_){print ";...]";return 0}
+      print "; ";sri=surd-1
+    }
+  }
+  print sign*cf__[i-1];
+  # detect a 0.0 which signifies that the cf has been determined
+  # to be infinite (not 100% accurate, but good)
+  if(scale(cf__[i])){
+    if(!surd||sri){print ","}else{print ";"}
+    print "..."
+  }
+  print "]";
+  return (i>cf_max);
+}
+
+# Print a number as a continued fraction
+define print_as_cf(x) {
+  auto cf[];
+  .=cf_new(cf[],x)+cf_print(cf[])
+  return x;
+}
+
+# Print a number as a signed continued fraction
+define print_as_cfn(x) {
+  auto cf[];
+  .=cfn_new(cf[],x)+cf_print(cf[])
+  return x;
+}
+
+# Output pbv array cf__[] as a fraction or combination of int and fraction
+# . the 'improper' parameter makes the choice.
+# . . set to non-zero for top-heavy / improper fractions
+define cf_print_frac(cf__[],improper) {
+  auto f[],v;
+  v=frac_from_cf(f[],cf__[],improper)
+  if(f[1]==0){print f[2];return f[2]}
+  if(!improper&&cf__[0]!=0){
+    print cf__[0]
+    if((f[1]<0)==(f[0]<0)){print"+"} # if n and d have same sign...
+  }
+  print f[1],"/",f[0] # n/d
+  return v;
+}
+
+# Print x as a fraction or combination of int and fraction
+# . the 'improper' parameter makes the choice.
+# . . set to non-zero for top-heavy / improper fractions
+define print_frac(x,improper) {
+  auto cf[]
+  .=cf_new(cf[],x)+cf_print_frac(cf[],improper)
+  return x;
+}
+
+# Print x as a fraction or combination of int and fraction
+# . the 'improper' parameter makes the choice.
+# . . set to non-zero for top-heavy / improper fractions
+# . This alternative function rounds to the nearest integer
+# . . for the integer part and then shows the addition
+# . . or _subtraction_ of the fraction.
+# . . e.g. 1+2/3 is shown as 2-1/3
+define print_frac2(x,improper) {
+  auto cf[];
+  .=cfn_new(cf[],x)+cf_print_frac(cf[],improper)
+  return x;
+}
+
+define cf_print_convergent(cf__[],c){
+  auto ocm,n,v;
+  if(c==0)return cf_print_frac(cf__[],0);
+  if(c<0){n=1;c=-c}else{n=c}
+  .=check_cf_max_();ocm=cf_max;
+  for(n=n;n<=c;n++){cf_max=n;v=cf_print_frac(cf__[],1);if(n!=c)print "\n"}
+  cf_max=ocm;return v;
+}
+
+define print_convergent(x,c){
+  auto cf[],v;
+  v=cf_new(cf[],x)
+  v=cf_print_convergent(cf[],c)
+  return v;
+}

cf_engel.bc

@@ -0,0 +1,107 @@
+#!/usr/local/bin/bc -l funcs.bc cf.bc
+
+### CF_ENGEL.BC - Engel expansion experimentation using array pass by reference
+
+# Workhorse: Create an Engel expansion of x in pbr array en__[]
+define engel_new_(mode) {
+  # expects *en__[] and x to be declared elsewhere
+  auto i,fc,max,q,p
+  
+  # error checking
+  .=check_cf_max_()
+  p=(mode<1&&abs(x)>=1);
+  q=(scale<6);
+  if(p||q){print "engel";if(mode==-1)print "fall";if(mode==0)print "alt";print "_new: "}
+  if(p){print "Can't work with integer part.\n"; x-=int(x)}
+  if(q){print "scale is ",scale,". Poor results likely.\n"}
+  
+  max=A^scale
+  i=fc=0;p=1
+  for(i=0;x&&p<max&&i<cf_max;i++){
+    q=1/abs(x)
+    if(mode==1){
+      q=ceil(q) # proper engel expansion
+    } else if(mode==-1){
+      q=floor(q) # secondary engel expansion
+    } else if(fc=!fc){q=floor(q)}else{q=ceil(q)} # tertiary engel expansion
+    p*=q
+    en__[i]=q*=sgn(x)
+    x=x*q-1
+  }
+  if(p>=max){
+    if(!--i)i=1
+    max=A^int(scale/2)
+    if(abs(en__[i])>max){
+      en__[i]=0;return 0
+    }else{
+      en__[i]=0.0
+      if(i-=2<3)return 0;
+      q=en__[i]
+      if(en__[--i]==q){
+        while(en__[i--]==q){}
+        en__[i+3]=0
+        .=en__[i+2]--
+      }
+    }
+  }else{
+    en__[i]=0
+  }
+  return 0;
+}
+
+# Create Engel expansion of x in pbr array en__[]
+# . all terms same sign as x
+define engel_new(*en__[],x) {
+  return x+engel_new_(1);
+}
+
+# Create secondary Engel expansion of x in pbr array en__[]
+# . first term same sign as x, all other terms opposite sign
+# . terms in implied Egyptian fraction sum have alternating signs
+define engelfall_new(*en__[],x) {
+  return x+engel_new_(-1);
+}
+
+# Create tertiary Engel expansion of x in pbr array en__[]
+# . first term same sign as x, following terms alternate in sign
+# . terms in implied Egyptian fraction alternate pairwise e.g. +--++--++...
+define engelalt_new(*en__[],x) {
+  return x+engel_new_(0);
+}
+
+# Output pbv array en__[] formatted as an Engel expansion
+define engel_print(en__[]){
+  auto i;
+  .=check_cf_max_();
+  print "{";
+  if(en__[1]==0){print en__[0],"}";return 0}
+  for(i=1;en__[i];i++)print en__[i-1],", ";
+  print en__[i-1];
+  if(scale(en__[i]))print ",...";
+  print "}";return 0;
+}
+
+# Print a number as an Engel expansion
+define print_as_engel(x){
+  auto en[];
+  .=engel_new(en[],x)+engel_print(en[])
+  return x;
+}
+define print_as_engelfall(x){
+  auto en[];
+  .=engelfall_new(en[],x)+engel_print(en[])
+  return x;
+}
+define print_as_engelalt(x){
+  auto en[];
+  .=engelalt_new(en[],x)+engel_print(en[])
+  return x;
+}
+
+# Turn the Engel expansion in pbv array en__[] into its value
+define engel_value(en__[]) {
+  auto i,p,v;
+  .=check_cf_max_();
+  p=1;v=0;for(i=0;i<cf_max&&p*=en__[i];i++)v+=1/p
+  return v;
+}

cf_misc.bc

@@ -0,0 +1,127 @@
+#!/usr/local/bin/bc -l funcs.bc cf.bc
+
+### CF_MISC.BC - Miscellaneous functions separated from CF.BC
+
+## CF Alteration ##
+
+# Take the absolute value of all terms in pbv array cf__[]
+# . WARNING: This irrevocably changes the array!
+define cf_abs_terms(*cf__[]) {
+  auto i,t,changed;
+  .=check_cf_max_();
+  changed=0;
+  for(i=0;i<cf_max&&cf__[i];i++){
+    t=cf__[i];if(t<0){t=-t;changed=1}
+    cf__[i]=t
+  }
+  return changed;
+}
+
+# Take the absolute value, less 1, of all terms in pbv array cf__[]
+# . WARNING: This irrevocably changes the array!
+define cf_abs1_terms(*cf__[]) {
+  auto i,err;
+  .=check_cf_max_();
+  for(i=0;i<cf_max&&cf__[i];i++)
+    if(0==(cf__[i]=abs(cf__[i])-1)){err=1;break}
+  if(err){
+    print "abs1_cf error: invalid cf for this transformation. truncation occurred\n";
+  }
+  return err;
+}
+
+# Return the value of a CF as if all terms are positive
+define cf_value_abs(cf__[]) {
+  auto cp[];
+  .=cf_copy(cp[],cf__[])+cf_abs_terms(cp[])
+  return cf_value(cp[])
+}
+
+# Return the value of a CF as if all terms are positive and reduced by 1
+define cf_value_abs1(cf__[]) {
+  auto cp[];
+  .=cf_copy(cp[],cf__[])+cf_abs1_terms(cp[])
+  return cf_value(cp[])
+}
+
+# Convert x through the cfn and abs transformations
+# . and return the value of the resultant CF
+define cfn_flip_abs(x) {
+  auto cf[];
+  .=cfn_new(cf[],x)+cf_abs_terms(cf[])
+  return cf_value(cf[])
+}
+
+# Convert x through the cfn and abs1 transformations
+# . and return the value of the resultant CF
+define cfn_flip_abs1(x) {
+  auto cf[];
+  .=cfn_new(cf[],x)+cf_abs1_terms(cf[])
+  return cf_value(cf[])
+}
+
+## Binary RLE <--> CF conversion ##
+
+# Minkowski Question Mark function - Inverse of Conway Box
+# . Treat the fractional part of x as a CF and transform it into a
+# . representation of alternating groups of bits in a binary number
+define cf_question_mark(cf__[]) {
+  auto os,n,i,b,x,t,tmax,sign,c;
+  .=check_cf_max_();
+  tmax=A^scale
+  sign=1;if(cf__[0]<0)sign=-1
+  b=0;t=1
+  for(i=1;i<cf_max&&cf__[i]&&t<tmax;i++){
+    if((c=cf__[i]*sign)<0){
+      print "question_mark_cf: terms not absolute values. aborting\n"
+      return 0;
+    }
+    for(j=c;j&&t<tmax;j--){x+=b/t;t+=t}
+    b=!b
+  }
+  os=scale
+  if(t<tmax){
+    c=0;while(t>1){.=c++;t/=A}
+    scale=c
+  }
+  x=(cf__[0]+sign*x)/1;
+  scale=os
+  return upscale_rational(x);
+}
+
+# As above but only generates a CF as intermediary
+define question_mark(x) { # returns ?(x)
+  auto cf[];
+  .=cf_new(cf[],x)
+  return cf_question_mark(cf[])
+}
+
+# Conway Box function - Inverse of Minkowski Question Mark
+# . Transform the fractional part of x by making a CF from a run-length
+# . encoding of the binary digits, and using that as the new fractional part
+define cf_conway_box(*cf__[],x) { # cf__[] = [[_x_]]
+  auto os,f[],max,p,i,b,bb,n,j,ix,sign,which0;
+  os=scale;scale+=scale
+  x=upscale_rational(x)
+  max=A^os;p=1
+  sign=1;if(x<0){sign=-1;x=-x}
+  x-=(b=int(x));cf__[0]=sign*b
+  b=0
+  n=1;j=1
+  for(i=0;i<=cf_max&&i<scale;i++){
+    x+=x;bb=int(x)
+    if(bb==b){.=n++}else{cf__[j++]=sign*n;p*=1+n;n=1}
+    b=bb;x-=b
+  }
+  if(n){cf__[j++]=sign*n;p*=1+n}
+  i=j;while(!i).=i--;.=i++
+  scale=os
+  return cf_tidy_();
+}
+
+# As above but only generates a CF as intermediary
+define conway_box(x) {
+  auto os,f[],cf[];
+  .=cf_conway_box(cf[],x)
+  return cf_value(cf[])
+}

cf_sylvester.bc

@@ -0,0 +1,92 @@
+#!/usr/local/bin/bc -l funcs.bc cf.bc
+
+### CF_SYLVESTER.BC - Sylvester expansion experimentation using array pass by reference
+
+# Workhorse: Create an Sylvester expansion of x in pbr array sy__[]
+define sylvester_new_(mode) {
+  # expects *sy__[] and x to be declared elsewhere
+  auto i,max,q,iq,p,h
+  
+  # error checking
+  .=check_cf_max_()
+  #p=(abs(x)>=1);
+  q=(scale<6);
+  #if(p||q){print "sylvester";if(mode)print "2";print "_new: "}
+  #if(p){print "Can't work with integer part.\n"; x-=int(x)}
+  if(q){print "scale is ",scale,". Poor results likely.\n"}
+  
+  max=A^int(scale/2)
+  i=0;p=1;h=3/2
+  for(i=0;x&&p<max&&i<cf_max;i++){
+    q=1/abs(x)
+    if(!mode){
+      q=ceil(q) # proper sylvester expansion
+    } else {
+      q=floor(q+h) # secondary sylvester expansion
+    }
+    p*=q
+    sy__[i]=q*=sgn(x)
+    x-=1/q
+  }
+  if(p>=max){
+    if(!mode)if(!--i)i=1
+    if(abs(sy__[i])>max){
+      sy__[i]=0;.=i--
+      while(iq=int(q=(sy__[i]*sy__[i-1])/(sy__[i]+sy__[i-1]))==q){
+        sy__[i]=0;sy__[--i]=iq
+      }
+    }else{
+      sy__[i]=0.0
+    }
+  }else{
+    sy__[i]=0
+  }
+  return 0;
+}
+#echo 'scale=250;sylvester_new(a[],7/15);.=sylvester_print(a[])+newline();sylvester_value(a[])' | bc -l funcs.bc *.bc
+
+
+# Create Sylvester expansion of x in pbr array sy__[]
+# . all terms same sign as x
+define sylvester_new(*sy__[],x) {
+  return x+sylvester_new_(0);
+}
+
+# Create secondary Sylvester expansion of x in pbr array sy__[]
+# . first term same sign as x, all other terms opposite sign
+# . terms in implied Egyptian fraction sum have alternating signs
+define sylvester2_new(*sy__[],x) {
+  return x+sylvester_new_(1);
+}
+
+# Output pbv array sy__[] formatted as an Sylvester expansion
+define sylvester_print(sy__[]){
+  auto i;
+  .=check_cf_max_();
+  print "{";
+  if(sy__[1]==0){print sy__[0],"}";return 0}
+  for(i=1;sy__[i];i++)print sy__[i-1],", ";
+  print sy__[i-1];
+  if(scale(sy__[i]))print ",...";
+  print "}";return 0;
+}
+
+# Print a number as an Sylvester expansion
+define print_as_sylvester(x){
+  auto sy[];
+  .=sylvester_new(sy[],x)+sylvester_print(sy[])
+  return x;
+}
+define print_as_sylvester2(x){
+  auto sy[];
+  .=sylvester2_new(sy[],x)+sylvester_print(sy[])
+  return x;
+}
+
+# Turn the Sylvester expansion in pbv array sy__[] into its value
+define sylvester_value(sy__[]) {
+  auto i,p,n,d;
+  .=check_cf_max_();
+  n=0;d=1;for(i=0;i<cf_max&&p=sy__[i];i++){n=n*p+d;d*=p}
+  return n/d;
+}

collatz.bc

@@ -0,0 +1,300 @@
+#!/usr/local/bin/bc -l
+
+### Collatz.BC - The 3x+1 or hailstones problem
+
+# Global variable
+# The original Collatz iteration has rules:
+#   odd  x -> 3x+1
+#   even x -> x/2
+# The condensed Collatz iteration has rules:
+#   odd  x -> (3x+1)/2
+#   even x -> x/2
+# ...since the usual odd step always produces an even value
+# The odd-only Collatz iteration has rules:
+#   odd  x -> odd part of 3x+1
+#   even x -> odd part of x
+# This var sets the mode of the functions in this library
+#   0 =>  odd-only Collatz
+#   1 =>  original Collatz - note that these two entries ...
+#   2 => condensed Collatz - ... match the divisor on the odd step
+collatz_mode_=1
+
+# sanity check
+define check_collatz_mode_() {
+  auto os;
+  if(collatz_mode_==0||collatz_mode_==1||collatz_mode_==2)return collatz_mode_
+  if(collatz_mode_<0||collatz_mode_>2)collatz_mode_=1
+  if(scale(collatz_mode_)){os=scale;scale=0;collatz_mode_/=1;scale=os}
+  return collatz_mode_
+}
+
+## Step forwards and back
+
+
+# Generate the next hailstone
+define collatz_next_(x) {
+  auto os,t;
+  os=scale;scale=0;x/=1
+  t=x/2;if(x!=t+t)t=3*x+1
+  if(collatz_mode_){
+    if(collatz_mode_==2&&t>x){x=t/2}else{x=t}
+  } else {
+    while(x==t+t||t>x){x=t;t/=2}
+  }
+  scale=os;return x
+}
+
+define collatz_next(x) {
+  .=check_collatz_mode_()
+  return collatz_next_(x)
+}
+
+# Take a guess at the previous hailstone - since in some cases there are
+# two choices, this function always chooses the option of lowest magnitude
+define collatz_prev(x) {
+  auto os,a,b,c;
+  os=scale;scale=0;x/=1
+  if(check_collatz_mode_()){
+    a=collatz_mode_*x-1;b=a/3
+    x+=x
+    if(3*b!=a||b==1||b==-1){scale=os;return x}
+    if((b>0)==(b<x))x=b
+  } else {
+    # oddonly mode shouldn't really return an even number
+    #  but when x is even or divisible by three, there _is_
+    #  no previous odd hailstone, so an even number must suffice.
+    if(!x%2||!x%3){scale=os;return x+x}
+    for(a=1;1;a+=a){
+      b=a*x-1;c=b/3
+      if(3*c==b){b=c/2;if(c!=b+b){scale=os;return c}}
+    }
+  }
+  scale=os;return x
+}
+
+## Chain examination
+
+max_array_ = 4^8
+
+# Determine whether an integer, x, reaches 1 under the Collatz iteration
+# . defined for both positive and negative x, so will
+# . return 0 under some circumstances!
+define is_collatz(x) {
+  auto os,t,i,tape[],tapetop
+  os=scale;scale=0;x/=1
+  if(x==0){scale=os;return 0}
+  .=check_collatz_mode_()
+  tapetop=-1
+  while(x!=1&&x!=-1){
+    t = collatz_next_(x)
+    # Search backwards for previous occurrence of t (which is more
+    #   likely to be near end of tape since chains lead to loops)
+    for(i=tapetop;i>0;i--)if(tape[i]==t){scale=os;return 0}
+    if(tapetop++>max_array_){
+      print "is_collatz: can't calculate; chain too long. assuming true.\n"
+      scale=os;return 1
+    }
+    tape[tapetop]=x=t
+  }
+  return x
+}
+
+# Print the chain of iterations of x until a loop or 1
+# . was cz_chain
+define collatz_print(x) {
+  auto os,t,i,tape[],tapetop
+  os=scale;scale=0;x/=1
+  x;if(x==0){scale=os;return 0}
+  .=check_collatz_mode_()
+  tapetop=-1
+  while(x!=1&&x!=-1){
+    t = collatz_next_(x)
+    # Search backwards for previous occurrence of t (which is more
+    #   likely to be near end of tape since chains lead to loops)
+    for(i=tapetop;i>0;i--)if(tape[i]==t){scale=os;"looping ";return t}
+    if(tapetop++>max_array_){
+      print "collatz_print: can't calculate; chain too long.\n"
+      scale=os;return t
+    }
+    tape[tapetop]=x=t;t
+  }
+}
+
+# Find the number of smallest magnitude under the Collatz iteration of x
+# . assuming the conjecture is true, this returns 1 for all positive x
+define collatz_root(x) {
+  auto os,t,i,tape[],tapetop
+  os=scale;scale=0;x/=1
+  if(x==0){scale=os;return 0}
+  .=check_collatz_mode_()
+  tapetop=-1
+  while(x!=1&&x!=-1){
+    t = collatz_next_(x)
+    # Search backwards for previous occurrence of t (which is more
+    #   likely to be near end of tape since chains lead to loops)
+    for(i=tapetop;i>0;i--)if(tape[i]==t){
+      #go back the other way looking for the lowest absolute value
+      while(++i<=tapetop)if((tape[i]>0)==(tape[i]<t))t=tape[i]
+      scale=os;return t
+    }
+    if(tapetop++>max_array_){
+      print "collatz_print: can't calculate; chain too long.\n"
+      scale=os;return (x>0)-(x<0)
+    }
+    tape[tapetop]=x=t
+  }
+  return x
+}
+
+# Returns the loopsize should the iteration become stuck in a loop
+# . assuming the conjecture is true, this returns 3 for the
+# . 4,2,1,4,etc. loop for all positive x.
+define collatz_loopsize(x) {
+  auto os,t,i,tape[],tapetop
+  os=scale;scale=0;x/=1
+  if(x==0){scale=os;return 1}
+  .=check_collatz_mode_()
+  tapetop=-1
+  while(x!=1&&x!=-1){
+    t = collatz_next_(x)
+    # Search backwards for previous occurrence of t (which is more
+    #   likely to be near end of tape since chains lead to loops)
+    for(i=tapetop;i>0;i--)if(tape[i]==t){scale=os;return tapetop-i+1}
+    if(tapetop++>max_array_){
+      print "collatz_loopsize: can't calculate; chain too long.\n"
+      scale=os;return 0
+    }
+    tape[tapetop]=x=t
+  }
+  if(collatz_mode_==0)return 1
+  if(collatz_mode_==1)return 3
+  if(collatz_mode_==2)return 2
+}
+
+# How many iterations to 1 (or loop)?
+define collatz_chainlength(x) {
+  auto os,t,i,c,tape[],tapetop
+  os=scale;scale=0;x/=1
+  if(x==0){scale=os;return 0}
+  .=check_collatz_mode_()
+  tapetop=-1
+  while(x!=1&&x!=-1){
+    .=c++
+    t = collatz_next_(x)
+    # Search backwards for previous occurrence of t (which is more
+    #   likely to be near end of tape since chains lead to loops)
+    for(i=tapetop;i>0;i--)if(tape[i]==t){scale=os;return 2-c }# infinity
+    if(tapetop++>max_array_){
+      print "collatz_chainlength: can't calculate; chain too long.\n"
+      scale=os;return -c
+    }
+    tape[tapetop]=x=t
+  }
+  return c
+}
+
+# Highest point on way to 1 or before being stuck in a loop
+define collatz_magnitude(x) {
+  auto os,t,i,m,tape[],tapetop
+  os=scale;scale=0;x/=1
+  if(x==0){scale=os;return 0}
+  .=check_collatz_mode_()
+  tapetop=-1
+  m=x