#define ROOT_DOUBLE 1
#define ROOT_TWO    2

/*
** Eigensolver for symmetric 2x2 matrix:
**
**  M00  M01
**  M01  M11
**
**
** Return value indicates something about the eigenvalue solution to
** the quadratic characteristic equation; see ROOT_ #defines above
** 
** HEY: the numerical precision issues here merit some more scrutiny.
*/
int Diderot_evals2x2 (
    float2 * eval,
    const Diderot_real_t _M00, const Diderot_real_t _M01,
    const Diderot_real_t _M11)
{
    Diderot_real_t mean, Q, D, M00, M01, M11;
    int roots;

  /* copy the given matrix elements */
    M00 = _M00;
    M01 = _M01;
    M11 = _M11;

  /*
  ** subtract out the eigenvalue mean (will add back to evals later);
  ** helps with numerical stability
  */
    mean = (M00 + M11)/2.0;
    M00 -= mean;
    M11 -= mean;
  
    Q = M00 - M11;
    D = 4.0*M01*M01 + Q*Q;
    if (D > EPSILON) {
      /* two distinct roots */
        Diderot_real_t vv;
        vv = sqrt(D)/2.0;
        eval->s0 = vv;
        eval->s1= -vv;
        roots = ROOT_TWO;
    }
    else {
      /* double root */
      eval->s0 = eval->s1 = 0.0;
      roots = ROOT_DOUBLE;
    }

  /* add back in the eigenvalue mean */
    eval->s0 += mean;
    eval->s1 += mean;

    return roots;
}

int Diderot_evecs2x2 (
    float2 * eval, float2 * evec,
    const Diderot_real_t _M00, const Diderot_real_t _M01,
    const Diderot_real_t _M11)
{
    Diderot_real_t mean, Q, D, M00, M01, M11;
    int roots;

  /* copy the given matrix elements */
    M00 = _M00;
    M01 = _M01;
    M11 = _M11;

  /*
  ** subtract out the eigenvalue mean (will add back to evals later);
  ** helps with numerical stability
  */
    mean = (M00 + M11)/2.0;
    M00 -= mean;
    M11 -= mean;
  
    Q = M00 - M11;
    D = 4.0*M01*M01 + Q*Q;
    if (D > EPSILON) {
      /* two distinct roots */
        Diderot_real_t vv;
        Diderot_vec2_t r1, r2;
        vv = sqrt(D)/2.0;
        eval->s0 = vv;
        eval->s1 = -vv;
      /* null space of T = M - evec[0]*I == 
         [M00 - vv      M01  ]
         [  M01      M11 - vv]
         is evec[0], but we know evec[0] and evec[1] are orthogonal,
         so row span of T is evec[1]
      */
        r1 = VEC2(M00 - vv, M01);
        r2 = VEC2(M01, M11 - vv);
        if (dot(r1,r2) > 0.0) {
            evec[1] = VEC2(r1.s0+r2.s0, r1.s1+r2.s1);
        }
        else {
            evec[1] = VEC2(r1.s0-r2.s0, r1.s1-r2.s1);
        }
        evec[1] = normalize(evec[1]);
        evec[0] = VEC2(evec[1].s1, -evec[1].s0);
        evec[0] = normalize(evec[0]);
        roots = ROOT_TWO;
    }
    else {
      /* double root */
        eval->s0 = eval->s1 = 0.0;
      /* use any basis for eigenvectors */
        evec[0] = VEC2(1.0, 0.0);
        evec[1] = VEC2(0.0, 1.0);
        roots = ROOT_DOUBLE;
    }

    /* add back in the eigenvalue mean */
    eval->s0 += mean;
    eval->s1 += mean;

    return roots;
}

<script type="text/javascript">//<![CDATA[
		function toggle_ffErrors() {
			var errorsblock = document.getElementById("ffErrorsBlock");
			var errorsgroup = document.getElementById("ffErrors");
			if (errorsblock.style.display == "none") {
				errorsblock.style.display = "block";
				errorsgroup.style.right = "10px";
			} else {
				errorsblock.style.display = "none";
				errorsgroup.style.right = "300px";
			}
		}
//]]></script>
<div id="ffErrors">
<a href="javascript:toggle_ffErrors();">Click to toggle</a>
<div id="ffErrorsBlock">
	<div class="error">does not end with &lt;/html&gt; tag</div>
	<div class="error">does not end with &lt;/body&gt; tag</div>
	<div class="info">The output has ended thus:  add back in the eigenvalue mean */
    eval-&gt;s0 += mean;
    eval-&gt;s1 += mean;

    return roots;
}</div>
</div></div>
</body></html>