[Teem-commit] teem/src/ell cubicEll.c,1.8,1.9 eigen.c,1.19,1.20

[Gordon K. <kin...@us...>]

Update of /cvsroot/teem/teem/src/ell
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv1306

Modified Files:
	cubicEll.c eigen.c 
Log Message:
based on recent observations, the task of sorting the cubic roots has been seen to be pretty trivial, so now the sorting is done in ell_cubic

Index: cubicEll.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ell/cubicEll.c,v
retrieving revision 1.8
retrieving revision 1.9
diff -C2 -d -r1.8 -r1.9
*** cubicEll.c	25 Mar 2004 03:12:52 -0000	1.8
--- cubicEll.c	1 Apr 2004 15:02:47 -0000	1.9
***************
*** 35,41 ****
  **   ellCubicRootTriple: root[0] == root[1] == root[2]
  **   ellCubicRootSingleDouble: single root[0]; double root[1] == root[2]
  **   ellCubicRootThree: root[0], root[1], root[2]
  **
! ** The values stored in root[] are NOT SORTED in anyway.
  **
  ** This does NOT use biff
--- 35,43 ----
  **   ellCubicRootTriple: root[0] == root[1] == root[2]
  **   ellCubicRootSingleDouble: single root[0]; double root[1] == root[2]
+ **                          or double root[0] == root[1], single root[2]
  **   ellCubicRootThree: root[0], root[1], root[2]
  **
! ** The values stored in root[] are, in a change from the past, sorted
! ** in descending order!  No need to sort them any more!
  **
  ** This does NOT use biff
***************
*** 58,64 ****
      theta = acos(R/sqrt(QQQ))/3.0;
      t = 2*sqrt(Q);
      root[0] = t*cos(theta) - sub;
!     root[1] = t*cos(theta + 2*AIR_PI/3.0) - sub;
!     root[2] = t*cos(theta - 2*AIR_PI/3.0) - sub;
      /*
      if (!AIR_EXISTS(root[0])) {
--- 60,68 ----
      theta = acos(R/sqrt(QQQ))/3.0;
      t = 2*sqrt(Q);
+     /* yes, these are sorted, because the C definition of acos says
+        that it returns values in in [0, pi] */
      root[0] = t*cos(theta) - sub;
!     root[1] = t*cos(theta - 2*AIR_PI/3.0) - sub;
!     root[2] = t*cos(theta + 2*AIR_PI/3.0) - sub;
      /*
      if (!AIR_EXISTS(root[0])) {
***************
*** 107,113 ****
  	fprintf(stderr, "%s: rescued double root:% 20.15f\n", me, nr);
        } 
!       root[0] = x;
!       root[1] = nr;
!       root[2] = nr;
        return ell_cubic_root_single_double;
      }
--- 111,123 ----
  	fprintf(stderr, "%s: rescued double root:% 20.15f\n", me, nr);
        } 
!       if (x > nr) {
! 	root[0] = x;
! 	root[1] = nr;
! 	root[2] = nr;
!       } else {
! 	root[0] = nr;
! 	root[1] = nr;
! 	root[2] = x;
!       }
        return ell_cubic_root_single_double;
      }
***************
*** 118,124 ****
        /* one double root and one single root */
        u = airCbrt(R);
!       root[0] = 2*u - sub;
!       root[1] = -u - sub;
!       root[2] = -u - sub;
        return ell_cubic_root_single_double;
      } 
--- 128,140 ----
        /* one double root and one single root */
        u = airCbrt(R);
!       if (u > 0) {
! 	root[0] = 2*u - sub;
! 	root[1] = -u - sub;
! 	root[2] = -u - sub;
!       } else {
! 	root[0] = -u - sub;
! 	root[1] = -u - sub;
! 	root[2] = 2*u - sub;
!       }
        return ell_cubic_root_single_double;
      } 
***************
*** 136,141 ****
  
  
- 
- 
- 
- 
--- 152,153 ----

Index: eigen.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ell/eigen.c,v
retrieving revision 1.19
retrieving revision 1.20
diff -C2 -d -r1.19 -r1.20
*** eigen.c	1 Mar 2004 21:58:55 -0000	1.19
--- eigen.c	1 Apr 2004 15:02:47 -0000	1.20
***************
*** 194,201 ****
      + (m[3]*m[1] - m[0]*m[4])*m[8];
    roots = ell_cubic(eval, A, B, C, newton);
!   if (ell_cubic_root_three == roots 
!       || ell_cubic_root_single_double == roots) {
!     ELL_SORT3(eval[0], eval[1], eval[2], tmp);
!   }
    ELL_3V_SCALE(_eval, 1.0/scale, eval);
    return roots;
--- 194,198 ----
      + (m[3]*m[1] - m[0]*m[4])*m[8];
    roots = ell_cubic(eval, A, B, C, newton);
!   /* no longer need to sort here */
    ELL_3V_SCALE(_eval, 1.0/scale, eval);
    return roots;
***************
*** 245,249 ****
    switch (roots) {
    case ell_cubic_root_three:
-     ELL_SORT3(e0, e1, e2, t);
      /* if (ell_debug) {
        printf("ell_3m_eigensolve: evals: %20.15f %20.15f %20.15f\n", 
--- 242,245 ----