~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/math/mathNum_eqn.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  * 
  6  * Redistribution and use in source and binary forms, with or without modification, 
  7  * are permitted provided that the following conditions are met:
  8  * 
  9  *   . Redistributions of source code must retain the above copyright notice, 
 10  *     this list of conditions and the following disclaimer.
 11  *    
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation 
 14  *     and/or other materials provided with the distribution.
 15  * 
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this 
 18  *     software without specific prior written permission.
 19  * 
 20  * 
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathNum_eqn.c,v $
 37  * Date    :  $Date: 2003/09/23 11:32:20 $
 38  * Version :  $Revision: 1.4 $
 39  * CVS Id  :  $Id: mathNum_eqn.c,v 1.4 2003/09/23 11:32:20 matts Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes : Quadratic & cubic roots
 44  *
 45  *********
 46 */
 47 /** 
 48  *  @file
 49  *  @brief Routines for finding the quadratic and cubic roots of a function.    
 50  *
 51  *
 52  * 
 53 */
 54 
 55 #include "mathNum_eqn.h"
 56 
 57 #if HAVE_CONFIG_H
 58 #include <config.h>
 59 #endif
 60 
 61 #include <math.h>
 62 #include <tina/math/math_GenDef.h>
 63 
 64 Bool            quadratic_roots(double a, double b, double c, double *x1, double *x2)
 65 {
 66         double          disc = b * b - 4.0 * a * c;
 67 
 68         /** assumes non-degenerate quadratic **/
 69 
 70         if (b < 0)
 71         {
 72                 a *= -1.0;
 73                 b *= -1.0;
 74                 c *= -1.0;
 75         }
 76         if (disc >= 0)
 77         {                       /**real roots **/
 78                 double          q = -0.5 * (b + sqrt(disc));
 79 
 80                 *x1 = q / a;
 81                 *x2 = c / q;
 82                 return (true);
 83         } else
 84         {
 85                 *x1 = -b / (2.0 * a);   /** real part **/
 86                 *x2 = sqrt(-disc) / (2.0 * a);  /** imag part **/
 87                 return (false);
 88         }
 89 }
 90 
 91 Bool            cubic_roots(double a, double b, double c, double d, double *x1, double *x2, double *x3)
 92 {
 93         double          disc, q, r;
 94 
 95         /** assumes non-degenerate cubic **/
 96         b /= a;
 97         c /= a;
 98         d /= a;
 99 
100         q = (b * b - 3.0 * c) / 9.0;
101         r = (2.0 * b * b * b - 9.0 * b * c + 27.0 * d) / 54.0;
102         disc = q * q * q - r * r;
103 
104         if (disc >= 0)
105         {
106                 double          theta;
107                 double          sq = sqrt(q), b3 = b / 3.0, sq3;
108 
109                 sq3 = sq * sq * sq;
110                 if (sq3 == 0.0)
111                         theta = 0.0;    /**arbitrary**/
112                 else
113                         theta = acos(r / sq3);
114 
115                 *x1 = -2.0 * sq * cos(theta / 3.0) - b3;
116                 *x2 = -2.0 * sq * cos((theta + TWOPI) / 3.0) - b3;
117                 *x3 = -2.0 * sq * cos((theta + 2.0 * TWOPI) / 3.0) - b3;
118                 return (true);
119         } else
120         {
121                 double          s = pow(fabs(r) + sqrt(-disc), 1.0 / 3.0);
122                 double          b3 = b / 3.0;
123                 double          x;
124 
125                 if (r >= 0.0)
126                         x = -(s + q / s);
127                 else
128                         x = (s + q / s);
129 
130                 x -= b3;
131 
132                 b += x;
133                 c += x * b;
134 
135                 *x1 = x;        /** real root **/
136                 return (quadratic_roots(1.0, b, c, x2, x3));    /** complex roots **/
137         }
138 }
139 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.