add missing copyrights

1 files changed, 107 insertions, 61 deletions

diff --git a/bench/eig33.cpp b/bench/eig33.cpp
index 2016c2c01..df07ad79d 100644
--- a/bench/eig33.cpp
+++ b/bench/eig33.cpp
@@ -1,3 +1,56 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+// The computeRoots function included in this is based on materials
+// covered by the following copyright and license:
+// 
+// Geometric Tools, LLC
+// Copyright (c) 1998-2010
+// Distributed under the Boost Software License, Version 1.0.
+// 
+// Permission is hereby granted, free of charge, to any person or organization
+// obtaining a copy of the software and accompanying documentation covered by
+// this license (the "Software") to use, reproduce, display, distribute,
+// execute, and transmit the Software, and to prepare derivative works of the
+// Software, and to permit third-parties to whom the Software is furnished to
+// do so, all subject to the following:
+// 
+// The copyright notices in the Software and this entire statement, including
+// the above license grant, this restriction and the following disclaimer,
+// must be included in all copies of the Software, in whole or in part, and
+// all derivative works of the Software, unless such copies or derivative
+// works are solely in the form of machine-executable object code generated by
+// a source language processor.
+// 
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
+// SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
+// FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
+// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+// DEALINGS IN THE SOFTWARE.
+
 #include <iostream>
 #include <Eigen/Core>
 #include <Eigen/Eigenvalues>
@@ -8,56 +61,49 @@ using namespace Eigen;
 using namespace std;
 
 template<typename Matrix, typename Roots>
-inline void computeRoots (const Matrix& rkA, Roots& adRoot)
+inline void computeRoots(const Matrix& m, Roots& roots)
 {
   typedef typename Matrix::Scalar Scalar;
-  const Scalar msInv3 = 1.0/3.0;
-  const Scalar msRoot3 = ei_sqrt(Scalar(3.0));
-
-  Scalar dA00 = rkA(0,0);
-  Scalar dA01 = rkA(0,1);
-  Scalar dA02 = rkA(0,2);
-  Scalar dA11 = rkA(1,1);
-  Scalar dA12 = rkA(1,2);
-  Scalar dA22 = rkA(2,2);
+  const Scalar s_inv3 = 1.0/3.0;
+  const Scalar s_sqrt3 = ei_sqrt(Scalar(3.0));
 
   // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0.  The
   // eigenvalues are the roots to this equation, all guaranteed to be
   // real-valued, because the matrix is symmetric.
-  Scalar dC0 = dA00*dA11*dA22 + Scalar(2)*dA01*dA02*dA12 - dA00*dA12*dA12 - dA11*dA02*dA02 - dA22*dA01*dA01;
-  Scalar dC1 = dA00*dA11 - dA01*dA01 + dA00*dA22 - dA02*dA02 + dA11*dA22 - dA12*dA12;
-  Scalar dC2 = dA00 + dA11 + dA22;
+  Scalar c0 = m(0,0)*m(1,1)*m(2,2) + Scalar(2)*m(0,1)*m(0,2)*m(1,2) - m(0,0)*m(1,2)*m(1,2) - m(1,1)*m(0,2)*m(0,2) - m(2,2)*m(0,1)*m(0,1);
+  Scalar c1 = m(0,0)*m(1,1) - m(0,1)*m(0,1) + m(0,0)*m(2,2) - m(0,2)*m(0,2) + m(1,1)*m(2,2) - m(1,2)*m(1,2);
+  Scalar c2 = m(0,0) + m(1,1) + m(2,2);
 
   // Construct the parameters used in classifying the roots of the equation
   // and in solving the equation for the roots in closed form.
-  Scalar dC2Div3 = dC2*msInv3;
-  Scalar dADiv3 = (dC1 - dC2*dC2Div3)*msInv3;
-  if (dADiv3 > Scalar(0))
-    dADiv3 = Scalar(0);
+  Scalar c2_over_3 = c2*s_inv3;
+  Scalar a_over_3 = (c1 - c2*c2_over_3)*s_inv3;
+  if (a_over_3 > Scalar(0))
+    a_over_3 = Scalar(0);
 
-  Scalar dMBDiv2 = Scalar(0.5)*(dC0 + dC2Div3*(Scalar(2)*dC2Div3*dC2Div3 - dC1));
+  Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
 
-  Scalar dQ = dMBDiv2*dMBDiv2 + dADiv3*dADiv3*dADiv3;
-  if (dQ > Scalar(0))
-    dQ = Scalar(0);
+  Scalar q = half_b*half_b + a_over_3*a_over_3*a_over_3;
+  if (q > Scalar(0))
+    q = Scalar(0);
 
   // Compute the eigenvalues by solving for the roots of the polynomial.
-  Scalar dMagnitude = ei_sqrt(-dADiv3);
-  Scalar dAngle = std::atan2(ei_sqrt(-dQ),dMBDiv2)*msInv3;
-  Scalar dCos = ei_cos(dAngle);
-  Scalar dSin = ei_sin(dAngle);
-  adRoot(0) = dC2Div3 + 2.f*dMagnitude*dCos;
-  adRoot(1) = dC2Div3 - dMagnitude*(dCos + msRoot3*dSin);
-  adRoot(2) = dC2Div3 - dMagnitude*(dCos - msRoot3*dSin);
+  Scalar rho = ei_sqrt(-a_over_3);
+  Scalar theta = std::atan2(ei_sqrt(-q),half_b)*s_inv3;
+  Scalar cos_theta = ei_cos(theta);
+  Scalar sin_theta = ei_sin(theta);
+  roots(0) = c2_over_3 + Scalar(2)*rho*cos_theta;
+  roots(1) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta);
+  roots(2) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta);
 
   // Sort in increasing order.
-  if (adRoot(0) >= adRoot(1))
-    std::swap(adRoot(0),adRoot(1));
-  if (adRoot(1) >= adRoot(2))
+  if (roots(0) >= roots(1))
+    std::swap(roots(0),roots(1));
+  if (roots(1) >= roots(2))
   {
-    std::swap(adRoot(1),adRoot(2));
-    if (adRoot(0) >= adRoot(1))
-      std::swap(adRoot(0),adRoot(1));
+    std::swap(roots(1),roots(2));
+    if (roots(0) >= roots(1))
+      std::swap(roots(0),roots(1));
   }
 }
 
@@ -65,21 +111,21 @@ template<typename Matrix, typename Vector>
 void eigen33(const Matrix& mat, Matrix& evecs, Vector& evals)
 {
   typedef typename Matrix::Scalar Scalar;
-    // Scale the matrix so its entries are in [-1,1].  The scaling is applied
-    // only when at least one matrix entry has magnitude larger than 1.
+  // Scale the matrix so its entries are in [-1,1].  The scaling is applied
+  // only when at least one matrix entry has magnitude larger than 1.
 
-    Scalar scale = mat.cwiseAbs()/*.template triangularView<Lower>()*/.maxCoeff();
-    scale = std::max(scale,Scalar(1));
-    Matrix scaledMat = mat / scale;
+  Scalar scale = mat.cwiseAbs()/*.template triangularView<Lower>()*/.maxCoeff();
+  scale = std::max(scale,Scalar(1));
+  Matrix scaledMat = mat / scale;
 
-    // Compute the eigenvalues
-//     scaledMat.setZero();
-    computeRoots(scaledMat,evals);
+  // Compute the eigenvalues
+//   scaledMat.setZero();
+  computeRoots(scaledMat,evals);
 
-    // compute the eigen vectors
-    // here we assume 3 differents eigenvalues
+  // compute the eigen vectors
+  // **here we assume 3 differents eigenvalues**
 
-    // "optimized version" which appears to be slower with gcc!
+  // "optimized version" which appears to be slower with gcc!
 //     Vector base;
 //     Scalar alpha, beta;
 //     base <<   scaledMat(1,0) * scaledMat(2,1),
@@ -93,22 +139,22 @@ void eigen33(const Matrix& mat, Matrix& evecs, Vector& evals)
 //     }
 //     evecs.col(2) = evecs.col(0).cross(evecs.col(1)).normalized();
 
-    // naive version
-    Matrix tmp;
-    tmp = scaledMat;
-    tmp.diagonal().array() -= evals(0);
-    evecs.col(0) = tmp.row(0).cross(tmp.row(1)).normalized();
-
-    tmp = scaledMat;
-    tmp.diagonal().array() -= evals(1);
-    evecs.col(1) = tmp.row(0).cross(tmp.row(1)).normalized();
-
-    tmp = scaledMat;
-    tmp.diagonal().array() -= evals(2);
-    evecs.col(2) = tmp.row(0).cross(tmp.row(1)).normalized();
-    
-    // Rescale back to the original size.
-    evals *= scale;
+  // naive version
+  Matrix tmp;
+  tmp = scaledMat;
+  tmp.diagonal().array() -= evals(0);
+  evecs.col(0) = tmp.row(0).cross(tmp.row(1)).normalized();
+
+  tmp = scaledMat;
+  tmp.diagonal().array() -= evals(1);
+  evecs.col(1) = tmp.row(0).cross(tmp.row(1)).normalized();
+
+  tmp = scaledMat;
+  tmp.diagonal().array() -= evals(2);
+  evecs.col(2) = tmp.row(0).cross(tmp.row(1)).normalized();
+  
+  // Rescale back to the original size.
+  evals *= scale;
 }
 
 int main()