source: trunk/Softwares/MiBench/src/c/basicmath-cubic.c @ 144

Last change on this file since 144 was 117, checked in by rosiere, 15 years ago

1) Platforms : add new organization for test
2) Load_Store_Unit : add array to count nb_check in store_queue
3) Issue_queue and Core_Glue : rewrite the issue network
4) Special_Register_Unit : add reset value to register CID
5) Softwares : add multicontext test
6) Softwares : add SPECINT
7) Softwares : add MiBench?
7) Read_queue : inhib access for r0
8) Change Core_Glue (network) - dont yet support priority and load balancing scheme
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 1.5 KB
Line 
1/* +++Date last modified: 05-Jul-1997 */
2
3/*
4**  CUBIC.C - Solve a cubic polynomial
5**  public domain by Ross Cottrell
6*/
7
8#include <math.h>
9#include <stdlib.h>
10#include "basicmath-snipmath.h"
11
12void SolveCubic(double  a,
13                double  b,
14                double  c,
15                double  d,
16                int    *solutions,
17                double *x)
18{
19      long double    a1 = b/a, a2 = c/a, a3 = d/a;
20      long double    Q = (a1*a1 - 3.0*a2)/9.0;
21      long double R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0;
22      double    R2_Q3 = R*R - Q*Q*Q;
23
24      double    theta;
25
26      if (R2_Q3 <= 0)
27      {
28            *solutions = 3;
29            theta = acos(R/sqrt(Q*Q*Q));
30            x[0] = -2.0*sqrt(Q)*cos(theta/3.0) - a1/3.0;
31            x[1] = -2.0*sqrt(Q)*cos((theta+2.0*PI)/3.0) - a1/3.0;
32            x[2] = -2.0*sqrt(Q)*cos((theta+4.0*PI)/3.0) - a1/3.0;
33      }
34      else
35      {
36            *solutions = 1;
37            x[0] = pow(sqrt(R2_Q3)+fabs(R), 1/3.0);
38            x[0] += Q/x[0];
39            x[0] *= (R < 0.0) ? 1 : -1;
40            x[0] -= a1/3.0;
41      }
42}
43
44#ifdef TEST
45
46int main(void)
47{
48      double  a1 = 1.0, b1 = -10.5, c1 = 32.0, d1 = -30.0;
49      double  a2 = 1.0, b2 = -4.5, c2 = 17.0, d2 = -30.0;
50      double  x[3];
51      int     solutions;
52
53      SolveCubic(a1, b1, c1, d1, &solutions, x);
54
55      /* should get 3 solutions: 2, 6 & 2.5   */
56
57      SolveCubic(a2, b2, c2, d2, &solutions, x);
58
59      /* should get 1 solution: 2.5           */
60
61      return 0;
62}
63
64#endif /* TEST */
Note: See TracBrowser for help on using the repository browser.