@x
\datethis
@y
\datethis
\let\maybe=\iffalse
@z
@x
the unary operators \<unop>, the binary operators \<binop>, and the
ternary operators \<ternop> are explained below. One short example
@y
the unary operators \<unop>, the binary operators \<binop>, and the
ternary operators \<ternop> are explained below. This extension
implements seven operations motivated by Olivier Coudert's
1997 conference paper about graph optimization with ZDDs.
There are three unary operations:
\.{\tth f0} (the maximal elements of $f_0$);
\.{\char`\_f0} (the minimal elements of $f_0$);
\.{\char`\#f0} (the minimal elements that meet all elements of $f_0$).
And four binary operations:
\.{f0;f1} (elements of $f_0$ that aren't supersets of elements of $f_1$);
\.{f0,f1} (elements of $f_0$ that aren't subsets of elements of $f_1$);
\.{f0:f1} (maximal elements of $f_0\sqcup f_1$);
\.{f0\char`\$f1} (maximal elements of $f_0\mathbin{;}f_1$).
One short example
@z
@x
delta~(11).
@y
delta~(11),
maxdot~(12),
maxnonsup~(13),
nonsup~(14),
nonsub~(15).
Unary operators `maxmal', `minmal', and `sharp' are also cached
with the binary opcode~3, with second operand respectively
|botsink|, |botsink+1|, and |botsink+2|.
@z
@x
char *binopname[]=
@y
char *unopname[]={"?","^","_","#"};
char *binopname[]=
@z
@x
@*Ternary operations. All operations can be reduced to binary operations,
@y
@ Here are some unary operations that interact with binary operations.
The first one extracts the maximal members of each family~$f$.
I'm (temporarily?)\ calling denoting it by~$f^{\uparrow}$.

@<Sub...@>=
node* maxmal_rec(node*f) {
  node *r,*r0,*r1;
  var *vf;
  if (f<=topsink) return oo,f->xref++,f;
    /* $\emptyset^{\uparrow}=\emptyset$, $\epsilon^{\uparrow}=\epsilon$ */
  oo,r=cache_lookup(f,botsink,node_(3));
  if (r) return r;
  @<Find $f^{\uparrow}$ recursively@>;
}

@ @<Find $f^{\uparrow}$ recursively@>=
rmems++;
vf=thevar(f); /* already fetched from memory when we did the cache lookup */
o,r=maxmal_rec(node_(f->lo));
if (!r) return NULL;
r1=maxmal_rec(node_(f->hi));
if (!r1) {
  deref(r);@+return NULL;
}
r0=nsub_rec(r,r1);
deref(r);
if (!r0) {
  deref(r1);@+return NULL;
}
r=unique_find(vf,r0,r1);
if (r) {
  if ((verbose&128)&&(vf<tvar))
    printf("   %x=^%x (level %d)\n",id(r),id(f),vf-varhead);
  cache_insert(f,botsink,node_(3),r);
}
return r;  

@ @<Templates for subroutines@>=
node* nsub_rec(node *f,node *g);
node* nsup_rec(node *f,node *g);

@ The |haseps| subroutine tests whether a given family contains
the empty set. The test is essentially `$f\cap\epsilon=\epsilon$?';
but I don't use |and_rec| for this special case, because |and_rec|
doesn't cache intermediate results at this stage of the computation.

@<Sub...@>=
int haseps(node *f) {
  node *r;
  int b;
  if (f<=topsink) return f-botsink;
  r=cache_lookup(f,topsink,node_(1));
  if (r) return oo,r->xref--,r-botsink;
  o,b=haseps(node_(f->lo));
  cache_insert(f,topsink,node_(1),botsink+b);
  return b;
}

@ The minimal elements of $f$, denoted $f^{\downarrow}$,
are found by an algorithm dual to the maximal one.

I'm not sure if the |haseps| call is really helpful here.
Hopefully I'll have time to experiment some day. If it actually
slows things down, it can safely be removed.

@<Sub...@>=
node* minmal_rec(node*f) {
  node *r,*r0,*r1;
  var *vf;
  if (f<=topsink) return oo,f->xref++,f;
    /* $\emptyset^{\uparrow}=\emptyset$, $\epsilon^{\uparrow}=\epsilon$ */
  if (haseps(f)) return oo,topsink->xref++,topsink;
  oo,r=cache_lookup(f,topsink,node_(3));
  if (r) return r;
  @<Find $f^{\downarrow}$ recursively@>;
}

@ @<Find $f^{\downarrow}$ recursively@>=
rmems++;
vf=thevar(f); /* already fetched from memory when we did the cache lookup */
o,r=minmal_rec(node_(f->hi));
if (!r) return NULL;
r0=minmal_rec(node_(f->lo));
if (!r0) {
  deref(r);@+return NULL;
}
r1=nsup_rec(r,r0);
deref(r);
if (!r1) {
  deref(r0);@+return NULL;
}
r=unique_find(vf,r0,r1);
if (r) {
  if ((verbose&128)&&(vf<tvar))
    printf("   %x=_%x (level %d)\n",id(r),id(f),vf-varhead);
  cache_insert(f,topsink,node_(3),r);
}
return r;  

@ The third unary operator computes $f^\sharp$, defined to be the
minimal elements of the family of all sets that have a nonzero
intersection with every member of~$f$. (In particular, if $f$
consists of the prime implicants of a monotone function,
$f^\sharp$ gives the prime clauses, and conversely.)

Again I am unsure about the value of |haseps| here.

@<Sub...@>=
node* sharp_rec(node*f) {
  node *r,*r0,*r1;
  var *vf;
  if (f<=topsink || haseps(f)) {
    if (f==botsink) return oo,topsink->xref++,topsink;
    else return oo,botsink->xref++,botsink;
  }
  oo,r=cache_lookup(f,botsink+2,node_(3));
  if (r) return r;
  @<Find $f^\sharp$ recursively@>;
}

@ @<Find $f^\sharp$ recursively@>=
rmems++;
vf=thevar(f); /* already fetched from memory when we did the cache lookup */
o,r=or_rec(node_(f->lo),node_(f->hi));
if (!r) return NULL;
r0=sharp_rec(r);
deref(r);
if (!r0) return NULL;
r=sharp_rec(node_(f->lo));
if (!r) {
  deref(r0);@+return NULL;
}
r1=nsup_rec(r,r0);
deref(r);
if (!r1) {
  deref(r0);@+return NULL;
}
r=unique_find(vf,r0,r1);
if (r) {
  if ((verbose&128)&&(vf<tvar))
    printf("   %x=#%x (level %d)\n",id(r),id(f),vf-varhead);
  cache_insert(f,botsink+2,node_(3),r);
}
return r;  

@ Now, the binary operations. The family |nsup(f,g)| consists of all
members of~$f$ that don't contain any members of~$g$.
In this documentation I call it $f\searrow g$.
(It can also be written $f\setminus(f\sqcup g)$; but this
implementation is intended to streamline that calculation.)

@<Sub...@>=
node* nsup_rec(node *f,node *g) {
  var *v,*vf,*vg;
  node *r,*r0,*r1;
  if (g==botsink) return oo,f->xref++,f;
  if (f==botsink || f==g || haseps(g)) return oo,botsink->xref++,botsink;
  oo,vf=thevar(f),vg=thevar(g);
  while (vf>vg) {
    oo,g=node_(g->lo),vg=thevar(g);
    if (g==botsink) return oo,f->xref++,f;
    if (f==g) return oo,botsink->xref++,botsink;
  }
  r=cache_lookup(f,g,node_(14));
  if (r) return r;
  @<Find $f\searrow g$ recursively@>;
}

@ @<Find $f\searrow g$ recursively@>=
rmems++; /* track recursion overhead */
v=vf;
if (v<vg) {
  o,r0=nsup_rec(node_(f->lo),g);
  if (!r0) return NULL;
  r1=nsup_rec(node_(f->hi),g);
  if (!r1) {
    deref(r0);@+return NULL;
  }
}@+else {
  oo,r0=nsup_rec(node_(f->hi),node_(g->hi));
  if (!r0) return NULL;
  r=nsup_rec(node_(f->hi),node_(g->lo));
  if (!r) {
    deref(r0);@+return NULL;
  }
  r1=and_rec(r,r0);
  deref(r);@+deref(r0);
  if (!r1) return NULL;
  r0=nsup_rec(node_(f->lo),node_(g->lo));
  if (!r0) {
    deref(r1);@+return NULL;
  }
}
r=unique_find(v,r0,r1);
if (r) {
  if ((verbose&128)&&(v<tvar))
    printf("   %x=%x;%x (level %d)\n",id(r),id(f),id(g),v-varhead);
  cache_insert(f,g,node_(14),r);
}
return r;  

@ Dually, |nsub_rec(f,g)| computes the ZDD for all
members of~$f$ that aren't contained in any members of~$g$.
In this documentation I call that family $f\nearrow g=
f\setminus(f\sqcap g)$.

@<Sub...@>=
node* nsub_rec(node *f,node *g) {
  var *v,*vf,*vg;
  node *r,*r0,*r1;
  if (g==botsink) return oo,f->xref++,f;
  if (f<=topsink || f==g) return oo,botsink->xref++,botsink;
  oo,r=cache_lookup(f,g,node_(15));
  if (r) return r;
  @<Find $f\nearrow g$ recursively@>;
}

@ @<Find $f\nearrow g$ recursively@>=
rmems++; /* track recursion overhead */
vf=thevar(f),v=vg=thevar(g); /* already fetched */
if (vf<v) {
  v=vf;
  o,r0=nsub_rec(node_(f->lo),g);
  if (!r0) return NULL;
  r1=node_(f->hi);
  r1->xref++;
}@+else {
  o,r1=nsub_rec((v==vf? o,node_(f->lo):f),node_(g->lo));
  if (!r1) return NULL;
  r=nsub_rec((v==vf? node_(f->lo):f),node_(g->hi));
  if (!r) {
    deref(r1);@+return NULL;
  }
  r0=and_rec(r,r1);
  deref(r);@+deref(r1);
  if (!r0) return NULL;
  r1=nsub_rec((v==vf? node_(f->hi):botsink),node_(g->hi));
  if (!r1) {
    deref(r0);@+return NULL;
  }
}
r=unique_find(v,r0,r1);
if (r) {
  if ((verbose&128)&&(v<tvar))
    printf("   %x=%x,%x (level %d)\n",id(r),id(f),id(g),v-varhead);
  cache_insert(f,g,node_(15),r);
}
return r;  

@ Coudert's |mdot_rec| routine computes the maximal elements of |prod_rec|.

@<Sub...@>=
node* mdot_rec(node*f, node*g) {
  var *v,*vf,*vg;
  node *r,*r0,*r1,*r01,*r10,*r11;
  if (f>g) r=f, f=g, g=r; /* wow */
  if (f<=topsink) {
    if (f==botsink) return oo,f->xref++,f; /* $0\odot g=0$ */
    else return maxmal_rec(g);
  }
  o,v=vf=thevar(f);
  o,vg=thevar(g);
  if (vf>vg) r=f, f=g, g=r, v=vg;
  r=cache_lookup(f,g,node_(12));
  if (r) return r;
  @<Find $f\odot g$ recursively@>;
}

@ @<Find $f\odot g$ recursively@>=
rmems++; /* track recursion overhead */
if (vf!=vg) {
  o,r0=mdot_rec(node_(f->lo),g);
  if (!r0) return NULL;
  r1=mdot_rec(node_(f->hi),g);
  if (!r1) {
    deref(r0); /* too bad, but we have to abort in midstream */
    return NULL;
  }
}@+else {
  o,r10=or_rec(node_(g->lo),node_(g->hi));
  if (!r10) return NULL;
  o,r1=mdot_rec(node_(f->hi),r10);
  deref(r10);
  if (!r1) return NULL;
  r01=mdot_rec(node_(f->lo),node_(g->hi));
  if (!r01) {
    deref(r1);@+return NULL;
  }
  r=or_rec(r01,r1);
  deref(r01);@+deref(r1);
  if (!r) return NULL;
  r1=maxmal_rec(r);
  deref(r);
  r0=mdot_rec(node_(f->lo),node_(g->lo));
  if (!r0) {
    deref(r1);@+return NULL;
  }
}
r=nsub_rec(r0,r1);
deref(r0);
if (!r) return NULL;
r=unique_find(v,r,r1);
if (r) {
  if ((verbose&128)&&(v<tvar))
    printf("   %x=%x:%x (level %d)\n",id(r),id(f),id(g),v-varhead);
  cache_insert(f,g,node_(12),r);
}
return r;  

@ Finally, |mnsup_rec| computes the maximals of |nsup_rec|. Coudert
claims this is more efficient than doing |nsup_rec| first, and I
tend to believe him; this routine will vindicate his judgment (or not).

@<Sub...@>=
node* mnsup_rec(node *f,node *g) {
  var *v,*vf,*vg;
  node *r,*r0,*r1;
  if (g==botsink) return maxmal_rec(f);
  if (f==botsink || f==g || haseps(g)) return oo,botsink->xref++,botsink;
  oo,vf=thevar(f),vg=thevar(g);
  while (vf>vg) {
    oo,g=node_(g->lo),vg=thevar(g);
    if (g==botsink) return maxmal_rec(f);
    if (f==g) return oo,botsink->xref++,botsink;
  }
  r=cache_lookup(f,g,node_(13));
  if (r) return r;
  @<Find $(f\searrow g)^\uparrow$ recursively@>;
}

@ @<Find $(f\searrow g)^\uparrow$ recursively@>=
rmems++; /* track recursion overhead */
v=vf;
if (v<vg) {
  o,r0=mnsup_rec(node_(f->lo),g);
  if (!r0) return NULL;
  r1=mnsup_rec(node_(f->hi),g);
  if (!r1) {
    deref(r0);@+return NULL;
  }
}@+else {
  o,r=or_rec(node_(g->lo),node_(g->hi));
  if (!r) return NULL;
  o,r1=mnsup_rec(node_(f->hi),r);
  deref(r);
  if (!r1) return NULL;
  r0=mnsup_rec(node_(f->lo),node_(g->lo));
  if (!r0) {
    deref(r1);@+return NULL;
  }
}
r=nsub_rec(r0,r1);
deref(r0);
if (!r) return NULL;
r=unique_find(v,r,r1);
if (r) {
  if ((verbose&128)&&(v<tvar))
    printf("   %x=%x$%x (level %d)\n",id(r),id(f),id(g),v-varhead);
  cache_insert(f,g,node_(13),r);
}
return r;  

@*Ternary operations. All operations can be reduced to binary operations,
@z
@x
case 11: r=delta_rec(f,g);@+break; /* $f\bindel g$ */
@y
case 11: r=delta_rec(f,g);@+break; /* $f\bindel g$ */
case 12: r=mdot_rec(f,g);@+break; /* $f\odot g$ */
case 13: r=mnsup_rec(f,g);@+break; /* $(f\searrow g)^\uparrow$ */
case 14: r=nsup_rec(f,g);@+break; /* $f\searrow g$ */
case 15: r=nsub_rec(f,g);@+break; /* $f\nearrow g$ */
@z
@x
@*Parsing the commands.
@y
@ @<Sub...@>=
node *unary_top(int curop,node *f) {
  node *r;
  unsigned long long oldmems=mems, oldrmems=rmems, oldzmems=zmems;
  if (verbose&2)
     printf("beginning to compute %s %x:\n",unopname[curop],id(f));
  cacheinserts=0;
  while (1) {
    switch (curop) {
case 1: r=maxmal_rec(f);@+break; /* $f^\uparrow$ */
case 2: r=minmal_rec(f);@+break; /* $f^\downarrow$ */
case 3: r=sharp_rec(f);@+break; /* $f^\sharp$ */
default: fprintf(stderr,"This can't happen!\n");@+exit(-69);
}
    if (r) break;
    attempt_repairs(); /* try to carry on */
  }
  if (verbose&(1+2))
    printf(" %x=%s%x (%llu mems, %llu rmems, %llu zmems, %.4g)\n",@|
                id(r),unopname[curop],id(f),
                mems-oldmems,rmems-oldrmems,zmems-oldzmems,@|
        mems-oldmems+rfactor*(rmems-oldrmems)+zfactor*(zmems-oldzmems));
  return r;
}

@*Parsing the commands.
@z
@x
case 'c': p=getconst(*c++);@+if (!p) reporterror;@+break;
@y
case 'c': p=getconst(*c++);@+if (!p) reporterror;@+break;
case '^': p=botsink;@+curop=-1;@+goto second; /* maxmal */
case '_': p=botsink;@+curop=-2;@+goto second; /* minmal */
case '#': p=botsink;@+curop=-3;@+goto second; /* sharp */
@z
@x
case '_': curop=11;@+break; /* delta */
@y
case '_': curop=11;@+break; /* delta */
case ':': curop=12;@+break; /* mdot */
case '$': curop=13;@+break; /* mnsup */
case ';': curop=14;@+break; /* nsup */
case ',': curop=15;@+break; /* nsub */
@z
@x
if (curop<=maxbinop) r=binary_top(curop,p,q);
@y
if (curop<0) r=unary_top(-curop,q);
else if (curop<=maxbinop) r=binary_top(curop,p,q);
@z