rpp_quintic.cpp
Go to the documentation of this file.
00001 /* ========================================================================
00002  * PROJECT: ARToolKitPlus
00003  * ========================================================================
00004  *
00005  * Solution of a quintic equation by a hybrid method
00006  * Conversion from Fortran to C by Raoul Rausch (Oct 2004)
00007  *
00008  * Copyright of the derived and new portions of this work
00009  *     (C) 2006 Graz University of Technology
00010  *
00011  * This framework is free software; you can redistribute it and/or modify
00012  * it under the terms of the GNU General Public License as published by
00013  * the Free Software Foundation; either version 2 of the License, or
00014  * (at your option) any later version.
00015  *
00016  * This framework is distributed in the hope that it will be useful,
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  * GNU General Public License for more details.
00020  *
00021  * You should have received a copy of the GNU General Public License
00022  * along with this framework; if not, write to the Free Software
00023  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00024  *
00025  * For further information please contact 
00026  *   Dieter Schmalstieg
00027  *   <schmalstieg@icg.tu-graz.ac.at>
00028  *   Graz University of Technology, 
00029  *   Institut for Computer Graphics and Vision,
00030  *   Inffeldgasse 16a, 8010 Graz, Austria.
00031  * ========================================================================
00032  ** @author   Thomas Pintaric
00033  *
00034  * $Id: rpp_quintic.cpp 162 2006-04-19 21:28:10Z grabner $
00035  * @file
00036  * ======================================================================== */
00037 
00038 
00039 #ifndef _NO_LIBRPP_
00040 
00041 
00042 #include "math.h"
00043 #include "stdio.h"
00044 #define max(a,b) (a>b?a:b)
00045 #define min(a,b) (a<b?a:b)
00046 
00047 
00048 namespace rpp {
00049 
00050 
00051 int quintic(double [], double [], double [], int*, double);
00052 int quartic(double[], double[], double[], int* );
00053 int cubic(double[], double[], int*);
00054 int signR(double);
00055 double CBRT(double);
00056 
00057 /*-------------------- Global Function Description Block ----------------------
00058  *
00059  *     ***QUINTIC************************************************25.03.98
00060  *     Solution of a quintic equation by a hybrid method:
00061  *     first real solution is obtained numerically by the Newton method,
00062  *     the remaining four roots are obtained analytically by QUARTIC
00063  *     NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING
00064  *     ******************************************************************
00065  *     dd(0:4)     (i)  vector containing the polynomial coefficients
00066  *     sol(1:4)    (o)  results, real part
00067  *     soli(1:4)   (o)  results, imaginary part
00068  *     Nsol        (o)  number of real solutions 
00069  *
00070  *     17-Oct-2004 / Raoul Rausch
00071  *     Conversion from Fortran to C
00072  *
00073  *
00074  *-----------------------------------------------------------------------------
00075  */
00076 /*
00077 int quintic(double dd[6], double sol[5], double soli[5], int *Nsol, double xstart)
00078 {
00079         double  dd4[5], sol4[4], soli4[4], xnew, xs;//, soli4[4];dd[6], sol[5], soli[5],
00080         double sum, sum1, eC;
00081         const double eps = 1.e-8;
00082         int i, Nsol4;
00083 
00084         *Nsol = 0;
00085 
00086         //printf("\n Quintic!\n");
00087 
00088         if (dd[5] == 0.0)
00089         { 
00090                 //printf("\n ERROR: NOT A QUINTIC EQUATION");
00091                 return 0;
00092         }
00093 
00094         // Newton iteration of one real root
00095         xs= xstart;
00096         xnew = xstart;  //added rr
00097         do
00098         {
00099                 xs = xnew;      //added rr
00100                 sum = dd[0];
00101                 for (i=1;i<6;i++)       sum += dd[i]*pow(xs,i); // Don't know what ** means
00102                 sum1 = dd[1];
00103                 for (i=1;i<5;i++)       sum1 += (double)(i+1)*dd[i+1]*pow(xs,i);
00104                 xnew = xs - sum/sum1;
00105                 //if (fabs(xnew-xs) > eps)
00106                 //xs =xnew;
00107                 //printf("\n %f\t%f!", xs, xnew);
00108         }while (fabs(xnew-xs) > eps);
00109 
00110         eC = xnew;
00111         //
00112         // "eC" is one real root of quintic equation
00113         // reduce quintic to quartic equation using "eC"
00114         dd4[4] = dd[5];
00115         for (i=4;i>0;i--)       dd4[i-1] = dd[i] + eC*dd4[i];
00116 
00117         quartic(dd4, sol4, soli4, &Nsol4);
00118 
00119         
00120         sol[0] = eC;
00121         soli[0] = 0.0;
00122 
00123         for (i=0;i<4;i++)
00124         {
00125                 sol[i+1] =sol4[i];
00126                 soli[i+1] = soli4[i];
00127         }
00128         *Nsol = Nsol4 + 1;
00129 
00130         return 0;
00131 }
00132 */
00133 
00134 /*-------------------- Global Function Description Block ----------------------
00135  *
00136  *     ***QUARTIC************************************************25.03.98
00137  *     Solution of a quartic equation
00138  *     ref.: J. E. Hacke, Amer. Math. Monthly, Vol. 48, 327-328, (1941)
00139  *     NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING
00140  *     ******************************************************************
00141  *     dd(0:4)     (i)  vector containing the polynomial coefficients
00142  *     sol(1:4)    (o)  results, real part
00143  *     soli(1:4)   (o)  results, imaginary part
00144  *     Nsol        (o)  number of real solutions 
00145  *     ==================================================================
00146  *      17-Oct-2004 / Raoul Rausch
00147  *              Conversion from Fortran to C
00148  *
00149  *
00150  *-----------------------------------------------------------------------------
00151  */
00152 
00153  int quartic(double dd[5], double sol[4], double soli[4], int* Nsol)
00154  {
00155         double AA[4], z[3];
00156         double a, b, c, d, f, p, q, r, zsol, xK2, xL, xK, sqp, sqm;
00157         int ncube, i;
00158         *Nsol = 0;
00159 
00160         if (dd[4] == 0.0)
00161         {
00162                 //printf("\n ERROR: NOT A QUARTIC EQUATION");
00163                 return 0;
00164         }
00165 
00166         a = dd[4];
00167         b = dd[3];
00168         c = dd[2];
00169         d = dd[1];
00170         f = dd[0];
00171 
00172         p = (-3.0*pow(b,2) + 8.0 *a*c)/(8.0*pow(a,2));
00173         q = (pow(b,3) - 4.0*a*b*c + 8.0 *d*pow(a,2)) / (8.0*pow(a,3));
00174         r = (-3.0*pow(b,4) + 16.0 *a*pow(b,2)*c - 64.0 *pow(a,2)*b*d + 256.0 *pow(a,3)*f)/(256.0*pow(a,4));
00175         
00176         // Solve cubic resolvent
00177         AA[3] = 8.0;
00178         AA[2] = -4.0*p;
00179         AA[1] = -8.0*r;
00180         AA[0] = 4.0*p*r - pow(q,2);
00181 
00182         //printf("\n bcubic %.4e\t%.4e\t%.4e\t%.4e ", AA[0], AA[1], AA[2], AA[3]);
00183         cubic(AA, z, &ncube);
00184         //printf("\n acubic %.4e\t%.4e\t%.4e ", z[0], z[1], z[2]);
00185         
00186         zsol = - 1.e99;
00187         for (i=0;i<ncube;i++)   zsol = max(zsol, z[i]); //Not sure C has max fct
00188         z[0] =zsol;
00189         xK2 = 2.0*z[0] -p;
00190         xK = sqrt(xK2);
00191         xL = q/(2.0*xK);
00192         sqp = xK2 - 4.0 * (z[0] + xL);
00193         sqm = xK2 - 4.0 * (z[0] - xL);
00194 
00195         for (i=0;i<4;i++)       soli[i] = 0.0;
00196         if ( (sqp >= 0.0) && (sqm >= 0.0))
00197         {
00198                 //printf("\n case 1 ");
00199                 sol[0] = 0.5 * (xK + sqrt(sqp));
00200                 sol[1] = 0.5 * (xK - sqrt(sqp));
00201                 sol[2] = 0.5 * (-xK + sqrt(sqm));
00202                 sol[3] = 0.5 * (-xK - sqrt(sqm));
00203                 *Nsol = 4;
00204         }
00205         else if ( (sqp >= 0.0) && (sqm < 0.0))
00206         {
00207                 //printf("\n case 2 ");
00208                 sol[0] = 0.5 * (xK + sqrt(sqp));
00209                 sol[1] = 0.5 * (xK - sqrt(sqp));
00210                 sol[2] = -0.5 * xK;
00211                 sol[3] = -0.5 * xK;
00212                 soli[2] =  sqrt(-.25 * sqm);
00213                 soli[3] = -sqrt(-.25 * sqm);
00214                 *Nsol = 2;
00215         }
00216         else if ( (sqp < 0.0) && (sqm >= 0.0))
00217         {
00218                 //printf("\n case 3 ");
00219                 sol[0] = 0.5 * (-xK + sqrt(sqm));
00220                 sol[1] = 0.5 * (-xK - sqrt(sqm));
00221                 sol[2] = 0.5 * xK;
00222                 sol[3] = 0.5 * xK;
00223                 soli[2] =  sqrt(-0.25 * sqp);
00224                 soli[3] = -sqrt(-0.25 * sqp);
00225                 *Nsol = 2;
00226         }
00227         else if ( (sqp < 0.0) && (sqm < 0.0))
00228         {
00229                 //printf("\n case 4 ");
00230                 sol[0] = -0.5 * xK;
00231                 sol[1] = -0.5 * xK;
00232                 soli[0] =  sqrt(-0.25 * sqm);
00233                 soli[1] = -sqrt(-0.25 * sqm);
00234                 sol[2] = 0.5 * xK;
00235                 sol[3] = 0.5 * xK;
00236                 soli[2] =  sqrt(-0.25 * sqp);
00237                 soli[3] = -sqrt(-0.25 * sqp);
00238                 *Nsol = 0;
00239         }
00240         
00241         for (i=0;i<4;i++)       sol[i] -= b/(4.0*a);
00242         return 0;
00243  }
00244 
00245 
00246  /*-------------------- Global Function Description Block ----------------------
00247   *
00248   *     ***CUBIC************************************************08.11.1986
00249   *     Solution of a cubic equation
00250   *     Equations of lesser degree are solved by the appropriate formulas.
00251   *     The solutions are arranged in ascending order.
00252   *     NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING
00253   *     ******************************************************************
00254   *     A(0:3)      (i)  vector containing the polynomial coefficients
00255   *     X(1:L)      (o)  results
00256   *     L           (o)  number of valid solutions (beginning with X(1))
00257   *     ==================================================================
00258   *     17-Oct-2004 / Raoul Rausch
00259   *             Conversion from Fortran to C
00260   *
00261   *
00262   *-----------------------------------------------------------------------------
00263   */
00264 
00265 int cubic(double A[4], double X[3], int* L)
00266 {
00267         const double PI = 3.1415926535897932;
00268         const double THIRD = 1./3.;
00269         double U[3],W, P, Q, DIS, PHI;
00270         int i;
00271 
00272         //define cubic root as statement function
00273         // In C, the function is defined outside of the cubic fct
00274 
00275         // ====determine the degree of the polynomial ====
00276 
00277         if (A[3] != 0.0)
00278         {
00279                 //cubic problem
00280                 W = A[2]/A[3]*THIRD;
00281                 P = pow((A[1]/A[3]*THIRD - pow(W,2)),3);
00282                 Q = -.5*(2.0*pow(W,3)-(A[1]*W-A[0])/A[3] );
00283                 DIS = pow(Q,2)+P;
00284                 if ( DIS < 0.0 )
00285                 {
00286                         //three real solutions!
00287                         //Confine the argument of ACOS to the interval [-1;1]!
00288                         PHI = acos(min(1.0,max(-1.0,Q/sqrt(-P))));
00289                         P=2.0*pow((-P),(5.e-1*THIRD));
00290                         for (i=0;i<3;i++)       U[i] = P*cos((PHI+2*((double)i)*PI)*THIRD)-W;
00291                         X[0] = min(U[0], min(U[1], U[2]));
00292                         X[1] = max(min(U[0], U[1]),max( min(U[0], U[2]), min(U[1], U[2])));
00293                         X[2] = max(U[0], max(U[1], U[2]));
00294                         *L = 3;
00295                 }
00296                 else
00297                 {
00298                         // only one real solution!
00299                         DIS = sqrt(DIS);
00300                         X[0] = CBRT(Q+DIS)+CBRT(Q-DIS)-W;
00301                         *L=1;
00302                 }
00303         }
00304         else if (A[2] != 0.0)
00305         {
00306                 // quadratic problem
00307                 P = 0.5*A[1]/A[2];
00308                 DIS = pow(P,2)-A[0]/A[2];
00309                 if (DIS > 0.0)
00310                 {
00311                         // 2 real solutions
00312                         X[0] = -P - sqrt(DIS);
00313                         X[1] = -P + sqrt(DIS);
00314                         *L=2;
00315                 }
00316                 else
00317                 {
00318                         // no real solution
00319                         *L=0;
00320                 }
00321         }
00322         else if (A[1] != 0.0)
00323         {
00324                 //linear equation
00325                 X[0] =A[0]/A[1];
00326                 *L=1;
00327         }
00328         else
00329         {
00330                 //no equation
00331                 *L=0;
00332         }
00333  //
00334  //     ==== perform one step of a newton iteration in order to minimize
00335  //          round-off errors ====
00336  //
00337         for (i=0;i<*L;i++)
00338         {
00339                 X[i] = X[i] - (A[0]+X[i]*(A[1]+X[i]*(A[2]+X[i]*A[3])))/(A[1]+X[i]*(2.0*A[2]+X[i]*3.0*A[3]));
00340         //      printf("\n X inside cubic %.15e\n", X[i]);
00341         }
00342 
00343         return 0;
00344 }
00345 
00346 
00347 int signR(double Z)
00348 {
00349         int ret = 0;
00350         if (Z > 0.0)    ret = 1;
00351         if (Z < 0.0)    ret = -1;
00352         if (Z == 0.0)   ret =0;
00353 
00354         return ret;
00355 }
00356 
00357 double CBRT(double Z)
00358 {
00359         double ret;
00360         const double THIRD = 1./3.;
00361         //define cubic root as statement function
00362         //SIGN has different meanings in both C and Fortran
00363         // Was unable to use the sign command of C, so wrote my own
00364         // that why a new variable needs to be introduced that keeps track of the sign of
00365         // SIGN is supposed to return a 1, -1 or 0 depending on what the sign of the argument is
00366         ret = fabs(pow(fabs(Z),THIRD)) * signR(Z);
00367         return ret;
00368 }
00369 
00370 
00371 }  // namespace rpp
00372 
00373 
00374 #endif //_NO_LIBRPP_


tuw_artoolkitplus
Author(s): Markus Bader
autogenerated on Sun May 29 2016 02:50:12