AAP aa_div_mix(AAP x, Float zeta, Float gamma, Float delta)
  {
    #if aa_TRACE_AFFINE
	fprintf(stderr, "enter aa_div_affine\n");
	fprintf(stderr, "  stack_top = %p\n", aa_stack_top);
	fprintf(stderr, "  x = %p  nterms = %ld  zeta  = %15.8e\n", x, x->nterms, zeta);
        fprintf(stderr, "  gamma = %15.8e  delta = %15.8e\n", gamma, delta);
        fprintf(stderr, "  ...\n");
    #endif

    if (
      aa_ISFULL(x) ||
      (FABS(zeta) <= Zero) ||
      (FABS(zeta) >= PlusInfinity) ||
      (FABS(gamma) >= PlusInfinity) ||
      (FABS(delta) >= PlusInfinity)
    )
      { return (aa_FULL); }
    else if (delta < Zero)
      { fatalerror ("aa_affine: delta is negative.");
        return(aa_FULL);
      }
    else
      { MemP frame = aa_stack_top;
        Float err = delta;
        AAP z = aa_alloc_head();

        /* Combine center values: */
        {
          Float xt, zt;
          
          if (zeta == One)
            { xt = (x->center); }
          else if (zeta = -One)
            { xt = -(x->center); }
          else
            { ROUND_NEAR; flt_div((x->center), zeta, &xt, &err); }
          
          if (gamma != Zero)
            { ROUND_NEAR; flt_add(zt, gamma, &zt, &err); }
          
          if ((FABS(zt) >= PlusInfinity) || (err >= PlusInfinity))
            { aa_flush(frame); return(aa_FULL); }
          
          (z->center) = zt;
        }
          
        /* Merge noise terms: */
        aa_unscale_terms (
          (AATermP) (x + 1), (x->nterms), zeta,
          &(z->nterms),
          &err
        );
        if (err >= PlusInfinity) 
          { aa_flush(frame); return(aa_FULL); }
        
        /* Add error term: */
        aa_append_error_term(&(z->nterms), err);

        /* check for zero: */
        if ((z->center == Zero) && (z->nterms == 0))
          { aa_flush(frame); return(aa_ZERO); }

        #if MIXED
        /* Compute z->range: */
        z->range = ia_meet(
          ia_div_affine(x->range, zeta, gamma, delta),
          aa_implicit_range(z)
        );
        #endif

        #if aa_TRACE_AFFINE
            fprintf(stderr, "  stack_top = %p\n", aa_stack_top);
            fprintf(stderr, "  z =     %p nterms = %ld\n", z, z->nterms);
            fprintf(stderr, "exit aa_affine.\n");
        #endif

        return(z);
      }
  }

void aa_unmix_terms(
    AATermP xp, AATermCount xn, 
    Float alpha, 
    AATermP yp, AATermCount yn, 
    Float beta, 
    AATermCount *znp, 
    Float *errp
  )
  /* Returns $zp[i] = xp[i] / \alpha + yp[i] / \beta$. */
  {
    #if aa_TRACE_SCALE
	fprintf(stderr, "enter aa_unmix_terms:\n");
	fprintf(stderr, "  stack_top = %p\n", aa_stack_top);
	fprintf(stderr, "  zn = %ld\n", (*znp));
	fprintf(stderr, "  xp = %p  xn = %ld  alpha = %15.8e\n ", xp, xn, alpha);
	fprintf(stderr, "  yp = %p  yn = %ld  beta  = %15.8e\n ", yp, yn, beta);
	fprintf(stderr, "  ...\n");
    #endif

    if (
      (FABS(alpha) == Zero) ||
      (FABS(beta) == Zero) ||
      ((*errp) >= PlusInfinity)
    )
      { (*errp) = PlusInfinity; }
    else if ((xn == 0) || (FABS(alpha) >= PlusInfinity))
      { aa_unscale_terms (yp, yn, beta, znp, errp); }
    else if ((yn == 0) || (FABS(beta) >= PlusInfinity))
      { aa_unscale_terms (xp, xn, beta, znp, errp); }
    else
      { MemP frame = aa_stack_top;
        AATermCount zn = 0;
        AATermP zp;
	Float xt, yt;

	while (xn > 0 || yn > 0)
          { zp = aa_alloc_term();
            /* Pick next term:*/
            if (yn == 0) 
              { zp->id = xp-> id; }
            else if (xn == 0)
              { zp->id = yp->id; }
            else if (xp->id < yp->id)
              { zp->id = xp->id; }
            else
              { zp->id = yp->id; }
              
            /* Compute $x$'s contribution: */
            if ((xn != 0) && (zp->id == xp->id))
	      { if (alpha == One) 
                  { xt = (xp->coef); }
                else if (alpha == -One }
                  { xt = -(xp->coef); }
                else
		  { ROUND_NEAR; flt_div(xp->coef, alpha, &xt, errp); }
		xp++; xn--;
	      }
            else
              { xt = Zero; }
	    
            /* Compute $y$'s contribution: */
            if ((yn != 0) && (zp->id == yp->id))
	      { if (beta == One) 
                  { yt = (yp->coef); }
                else if (beta == -One }
                  { yt = -(yp->coef); }
                else
		  { ROUND_NEAR; flt_div(yp->coef, beta, &yt, errp); }
		yp++; yn--;
	      }
	    else
	      { yt = Zero; }
              
            /* Add the two terms: */
	    if (xt == Zero)
              { (zp->coef) = yt; }
            else if (yt == Zero)
              { (zp->coef) = xt; }
            else
              { ROUND_NEAR; flt_add(xt, yt, &(zp->coef), errp); }

	    if (FABS(zp->coef) >= Infinity || (*errp) >= Infinity)
	      { (*errp) = PlusInfinity; aa_flush(frame); return; }
	    if ((zp->coef) == Zero) aa_pop_term(); else zn++;
	  }
        }
        (*znp) = zn;
      }
    #if aa_TRACE_SCALE
	fprintf(stderr, "  stack_top = %p\n", aa_stack_top);
	fprintf(stderr, "  zn = %ld\n", (*znp));
	fprintf(stderr, "exit aa_unmix_terms.\n");
    #endif
  }