home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c129 / 1.ddi / ROOTS.C < prev    next >
Encoding:
C/C++ Source or Header  |  1988-01-21  |  5.6 KB  |  268 lines
  1. # include <math.h>
  2. float AFunction(float X)
  3. {
  4.    float _fret;
  5.  
  6.    _fret = (2.0 * X + 3.0) * (X - 3.0);
  7.    return(_fret);
  8. }
  9.  
  10. float AFunctionPrime(float X)
  11. {
  12.    float _fret;
  13.  
  14.    _fret = 4.0 * X - 3.0;
  15.    return(_fret);
  16. }
  17.  
  18. char RootBracketed(float x1, float x2)
  19. {
  20.    char _fret;
  21.    char result;
  22.  
  23.    if ( ((x1 > 0.0) && (x2 > 0.0)) || ((x1 < 0.0) && (x2 < 0.0)) ) {
  24.       result = 0;
  25.    }
  26.    else {
  27.       result = 1;
  28.    }
  29.    _fret = result;
  30.    return(_fret);
  31. }
  32.  
  33. float Minimum(float x1, float x2)
  34. {
  35.    float _fret;
  36.    float result;
  37.  
  38.    if ( x1 < x2 ) {
  39.       result = x1;
  40.    }
  41.    else {
  42.       result = x2;
  43.    }
  44.    _fret = result;
  45.    return(_fret);
  46. }
  47.  
  48. float BrentRoots(float x1,float x2,float Tolerance,
  49.                   int maxIterations,float *valueAtRoot,int *error)
  50. {
  51.    float _fret;
  52.    # define FPP 1.0e-11
  53.    # define nearzero 1.0e-10
  54.    float result;
  55.    float AA;
  56.    float BB;
  57.    float CC;
  58.    float DD;
  59.    float EE;
  60.    float FA;
  61.    float FB;
  62.    float FC;
  63.    float Tol1;
  64.    float PP;
  65.    float QQ;
  66.    float RR;
  67.    float SS;
  68.    float XM;
  69.    int i;
  70.    int j;
  71.    int k;
  72.    char done;
  73.  
  74.    i = 0;
  75.    done = 0;
  76.    (*error) = 0;
  77.    AA = x1;
  78.    BB = x2;
  79.    FA = AFunction(AA);
  80.    FB = AFunction(BB);
  81.    if ( ! (RootBracketed(FA,FB)) ) {
  82.       (*error) = 1;
  83.    }
  84.    else {
  85.       FC = FB;
  86.       do {
  87.          if ( ! (RootBracketed(FC,FB)) ) {
  88.             CC = AA;
  89.             FC = FA;
  90.             DD = BB - AA;
  91.             EE = DD;
  92.          }
  93.          if ( fabs(FC) < fabs(FB) ) {
  94.             AA = BB;
  95.             BB = CC;
  96.             CC = AA;
  97.             FA = FB;
  98.             FB = FC;
  99.             FC = FA;
  100.          }
  101.          Tol1 = 2.0 * FPP * fabs(BB) + 0.5 * Tolerance;
  102.          XM = 0.5 * (CC - BB);
  103.          if ( (fabs(XM) <= Tol1) || (fabs(FA) < nearzero) ) {
  104.             result = BB;
  105.             done = 1;
  106.             (*valueAtRoot) = AFunction(result);
  107.          }
  108.          else {
  109.             if ( (fabs(EE) >= Tol1) && (fabs(FA) > fabs(FB)) ) {
  110.                SS = FB / FA;
  111.                if ( fabs(AA - CC) < nearzero ) {
  112.                   PP = 2.0 * XM * SS;
  113.                   QQ = 1.0 - SS;
  114.                }
  115.                else {
  116.                   QQ = FA / FC;
  117.                   RR = FB / FC;
  118.                   PP = SS * (2.0 * XM * QQ * (QQ - RR) - (BB - AA) * (RR - 1.0));
  119.                   QQ = (QQ - 1.0) * (RR - 1.0) * (SS - 1.0);
  120.                }
  121.                if ( PP > nearzero ) {
  122.                   QQ = -QQ;
  123.                }
  124.                PP = fabs(PP);
  125.                if ( (2.0 * PP) < Minimum(3.0 * XM * QQ - fabs(Tol1 * QQ),fabs(EE * QQ)) ) {
  126.                   EE = DD;
  127.                   DD = PP / QQ;
  128.                }
  129.                else {
  130.                   DD = XM;
  131.                   EE = DD;
  132.                }
  133.             }
  134.             else {
  135.                DD = XM;
  136.                EE = DD;
  137.             }
  138.             AA = BB;
  139.             FA = FB;
  140.             if ( fabs(DD) > Tol1 ) {
  141.                BB = BB + DD;
  142.             }
  143.             else {
  144.                if ( XM > 0.0 ) {
  145.                   BB = BB + fabs(Tol1);
  146.                }
  147.                else {
  148.                   BB = BB - fabs(Tol1);
  149.                }
  150.             }
  151.             FB = AFunction(BB);
  152.             i = 2;
  153.          }
  154.       } while ( ! (done || (i == maxIterations)) );
  155.       if ( i == maxIterations ) {
  156.          (*error) = 2;
  157.       }
  158.    }
  159.    _fret = result;
  160.    return(_fret);
  161. }
  162.  
  163. float BisectionRoots(float x1, float x2, float Tolerance,
  164.                       int maxIterations,float *valueAtRoot,int *error)
  165. {
  166.    float _fret;
  167.    # define FPP 1.0e-11
  168.    # define nearzero 1.0e-10
  169.    float result;
  170.    float FA;
  171.    float FB;
  172.    float FC;
  173.    float XA;
  174.    float XB;
  175.    float XC;
  176.    int i;
  177.    int j;
  178.    int k;
  179.    char done;
  180.  
  181.    (*error) = 0;
  182.    i = 0;
  183.    done = 0;
  184.    XA = x1;
  185.    XB = x2;
  186.    FA = AFunction(XA);
  187.    FB = AFunction(XB);
  188.    if ( ! (RootBracketed(FA,FB)) ) {
  189.       (*error) = 1;
  190.    }
  191.    else {
  192.       do {
  193.          XC = (XA + XB) / 2.0;
  194.          FC = AFunction(XC);
  195.          if ( RootBracketed(FA,FC) ) {
  196.             XB = XC;
  197.          }
  198.          else {
  199.             XA = XC;
  200.             FA = FC;
  201.          }
  202.          if ( fabs(XB - XA) <= Tolerance ) {
  203.             done = 1;
  204.             result = (XA + XB) / 2.0;
  205.             (*valueAtRoot) = AFunction(result);
  206.          }
  207.          i = i + 1;
  208.       } while ( ! (done || (i == maxIterations)) );
  209.       if ( i == maxIterations ) {
  210.          (*error) = 2;
  211.       }
  212.    }
  213.    _fret = result;
  214.    return(_fret);
  215. }
  216.  
  217. float NewtonRoots(float x1,float Tolerance,float maxIterations,
  218.                    float *valueAtRoot, int *error)
  219. {
  220.    float _fret;
  221.    # define FPP 1.0e-11
  222.    # define nearzero 1.0e-10
  223.    float result;
  224.    float lastX;
  225.    float FX;
  226.    float FPX;
  227.    float X;
  228.    float XP;
  229.    int i;
  230.    int j;
  231.    int k;
  232.    char done;
  233.  
  234.    i = 0;
  235.    (*error) = 0;
  236.    done = 0;
  237.    X = x1;
  238.    FX = AFunction(X);
  239.    do {
  240.       if ( fabs(X) < nearzero ) {
  241.          (*error) = 1;
  242.       }
  243.       else {
  244.          FPX = AFunctionPrime(X);
  245.          if ( fabs(FPX) < nearzero ) {
  246.             (*error) = 2;
  247.          }
  248.          else {
  249.             lastX = X;
  250.             X = X - FX / FPX;
  251.             FX = AFunction(X);
  252.             if ( fabs(X - lastX) <= Tolerance ) {
  253.                done = 1;
  254.                result = X;
  255.                (*valueAtRoot) = AFunction(X);
  256.             }
  257.             i = i + 1;
  258.          }
  259.       }
  260.    } while ( ! (done || (i == maxIterations) || ((*error) != 0)) );
  261.    if ( i == maxIterations ) {
  262.       (*error) = 3;
  263.    }
  264.    _fret = result;
  265.    return(_fret);
  266. }
  267.  
  268.