home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / spline29 / zero.c < prev   
Encoding:
C/C++ Source or Header  |  1992-04-17  |  16.9 KB  |  653 lines
  1. /*
  2.  
  3.     zero - find a zero of a function evaluated by the calling program
  4.  
  5.     SYNOPSIS...
  6.         inversion = init_zero(&x, inversion);    -- initialize
  7.         advise_zero(x2, inversion);        -- suggest 2nd x value
  8.         ...
  9.         err = zero(y, inversion);        -- ask zero() to update x
  10.  
  11.     USAGE...
  12.         init_zero() performs initialization.  Its first argument points
  13.         to the independent variable.  init_zero() returns a pointer to
  14.         a record allocated from the heap which contains variables
  15.         required by zero().  If the second argument is non-NULL it is
  16.         assumed to be a pointer to a previously allocated record.  In
  17.         that case the old record is reused.  If the allocation fails,
  18.         init_zero() returns NULL.
  19.         
  20.         zero() performs the actual function inversion.  Its first
  21.         argument is the value of the function for the current value of
  22.         the independent variable.  The second argument is a pointer to
  23.         the variable record which was allocated and initialized by
  24.         init_zero().  zero() updates the independent variable in an
  25.         attempt to force the function value to zero.  It returns an
  26.         error code (see below for values).  
  27.  
  28.         The calling program may reset the independent variable before
  29.         evaluating the function.  This would be useful for initiation
  30.         (to "show" zero() points straddling the desired value) or if
  31.         zero() happened to set the independent variable just outside
  32.         its legal range.  However, use of the same value twice in a row
  33.         will result in an error (code 5).
  34.  
  35.     EXAMPLE...
  36.  
  37.     double target;            -- desired function value
  38.     double x=START_VALUE;    -- independent variable: initialized here, 
  39.                                                 updated by zero()
  40.     double y;                -- function value, evaluated here
  41.     int error;                -- error status, set nonzero by zero()...
  42.                                 0: ok so far
  43.                                 1: local extremum
  44.                                 2: function value not changing 
  45.                                 3: not converging fast enough 
  46.                                 4: repeated x value
  47.                                 5: null inversion record pointer
  48.  
  49.     char *inversion=0, *init_zero();
  50.  
  51.             ...
  52.  
  53.     inversion = init_zero(&x, inversion);    -- initialize
  54.     do
  55.         {y = (x + 1.)*x*(x - 1.);            -- evaluate the function
  56.         del = y - target;
  57.         if(fabs(del) < .00000001) break;
  58.         error = zero(del, inversion);        -- ask zero() to update x
  59.         } while(!error);
  60.  
  61.     or simply
  62.  
  63.     inversion = init_zero(&x, 0);            -- initialize
  64.     do
  65.         {y = (x + 1.)*x*(x - 1.);            -- evaluate the function
  66.         if(fabs(y - target) < .00000001) break;
  67.         } while(!zero(y - target, inversion));
  68.  
  69.     NOTES...
  70.         A programmer might state his problem as follows: "the result of
  71.         this computation is <value>, but I want to adjust <variable> to
  72.         move the result closer to <target>" In such a case, zero()
  73.         can be used to adjust <variable>.
  74.  
  75.         Note that only the interpolation logic, not the function
  76.         evaluation, is relegated to a subprogram.  Thus, zero() can
  77.         invert an expression without the programmer having to code it
  78.         as a separate function.  This can simplify the source code,
  79.         since the programmer need not create separate functions just
  80.         for the evaluation of simple expressions.
  81.  
  82.         The programmer also maintains full control - he codes the
  83.         stopping criterion and can even reject the values proposed by
  84.         zero().
  85.  
  86.         advise_zero() can be used to suggest a second value of x. 
  87.         This may speed the initial part of the search, which normally
  88.         uses small steps for ruggedness.
  89.  
  90.     METHOD...
  91.         The search procedes in two phases.  In the first phase, while
  92.         all y values are of the same sign, zero() uses two previous x,y
  93.         pairs to predict the desired x value (a secant search).  This
  94.         is an extrapolation and the x step size is restricted to help
  95.         prevent divergence.  At most, the step size can double on each
  96.         iteration.  The new x,y pair replaces the one further from the
  97.         target.
  98.  
  99.         As soon as both positive and negative y values are found,
  100.         zero() switches to interpolation which is much more reliable. 
  101.         In general, zero() performs quadratic interpolation using three
  102.         previous x,y pairs to speed convergence (usually doubling the
  103.         number of correct bits each iteration).  In this case, the new
  104.         x,y pair replaces the existing pair on the same side of the
  105.         target.  Therefore, zero() will never switch from interpolating
  106.         back to extrapolating.
  107.  
  108.     IMPLEMENTATION...
  109.         zero() must save its state between returns so it can pick up
  110.         where it left off.  Recursive calls may also be needed - when
  111.         the calling program evaluates its function, it may use zero()
  112.         to do so.  Therefore, zero() needs an independent collection
  113.         of variables for each inversion being performed.  Init_zero()
  114.         allocates storage for this, of type INVERSION_REC, in the heap. 
  115.         An INVERSION_REC may be considered an "object", init_zero()
  116.         its "constructor", and calls to zero() its "messages".  Its
  117.         "destructor" is free().  However, reuse is encouraged (see the
  118.         example).
  119.  
  120. */
  121.  
  122. #include <math.h>
  123. #ifdef MAIN
  124. #include <assert.h>
  125. #define __assertfail printf
  126. #include <stdio.h>
  127. /*
  128. #else
  129. #define assert(x) 0
  130. */
  131. #endif
  132.  
  133. #ifdef MAIN
  134. int verbose = 1;
  135. #endif
  136.  
  137. #define ITERATIONS 50
  138. #define MAX_REGISTERED 4
  139.  
  140. typedef struct
  141.     {
  142.     int state_rec, straddling_rec, quadratic_rec;
  143.     int i_rec, *error_rec;
  144.     double extrapolation_limit_rec;
  145.     double x1_rec, x2_rec, x3_rec;
  146.     double f1_rec, f2_rec, f3_rec;
  147.     double d1_rec, d2_rec, d3_rec;
  148.     double del_rec, dx_rec, slope_rec;
  149.     double *xp_rec, difference_rec;
  150.     int suggestion_cnt;
  151.     double suggestion[MAX_REGISTERED];
  152.     } INVERSION_REC;
  153.  
  154. enum error_status{OK, LOCAL_EXTREMUM, NOT_CHANGING, SLOW_CONVERGENCE,
  155.         REPEATED_X, NULL_POINTER, OUT_OF_MEMORY};
  156.  
  157. #define RR(x, s) /* printf("%20.5g  %s\n", x, s);  /**/
  158.  
  159. #define state                irp->state_rec
  160. #define straddling            irp->straddling_rec
  161. #define quadratic            irp->quadratic_rec
  162. #define i                    irp->i_rec
  163. #define error                irp->error_rec
  164. #define x1                    irp->x1_rec
  165. #define x2                    irp->x2_rec
  166. #define x3                    irp->x3_rec
  167. #define f1                    irp->f1_rec
  168. #define f2                    irp->f2_rec
  169. #define f3                    irp->f3_rec
  170. #define d1                    irp->d1_rec
  171. #define d2                    irp->d2_rec
  172. #define d3                    irp->d3_rec
  173. #define del                    irp->del_rec
  174. #define difference            irp->difference_rec
  175. #define dx                    irp->dx_rec
  176. #define extrapolation_limit    irp->extrapolation_limit_rec
  177. #define slope                irp->slope_rec
  178. #define xp                    irp->xp_rec
  179.  
  180. static int junk;
  181.  
  182. INVERSION_REC *
  183. init_zero(xp_, former_rec) 
  184. double *xp_;         /* ptr to user's variable which holds initial value of 
  185.                         function's independent variable, updated by zero() */
  186.  
  187. INVERSION_REC *former_rec;    /* points to previous record (which is reused), 
  188.                                 or NULL (so that a new record is malloc'ed) */
  189. {    INVERSION_REC *irp, *malloc();
  190.  
  191.     if(former_rec) irp = former_rec;
  192.     else irp = malloc(sizeof(INVERSION_REC));
  193.     if(irp)
  194.         {
  195.         xp = xp_;
  196.         straddling = 0;
  197.         extrapolation_limit = 2.;
  198.         state = 1;
  199.         irp->suggestion_cnt = 0;
  200.         }
  201.     return irp;
  202. }
  203.  
  204. advise_zero(value, irp) double value; INVERSION_REC *irp;
  205. {    if(irp->suggestion_cnt < MAX_REGISTERED)
  206.         irp->suggestion[irp->suggestion_cnt++] = value;    
  207. }
  208.  
  209. zero(value, irp) double value; INVERSION_REC *irp;
  210. {    double x31, x32, x312;
  211.     register double temp;
  212.  
  213.     if(!irp)
  214.         {return NULL_POINTER;            /* FAIL - null pointer */
  215.         }
  216.  
  217.     switch(state)
  218.         {
  219.         case 1:     goto loc1;
  220.         case 2:     goto loc2;
  221.         default:    goto loc3;
  222.         }
  223.  
  224. loc1:
  225.  
  226.     x1 = *xp;
  227.     f1 = value;
  228.     d1 = fabs(f1);
  229.     if(x1 == 0.) del = .01; else del = .01*x1;
  230.  
  231.  
  232.                                 /*
  233.                                     first, look for two x's with 
  234.                                     different function values 
  235.                                 */
  236.     for (i = 0; i < ITERATIONS; i++)    
  237.         {
  238.         if(irp->suggestion_cnt) 
  239.             {*xp = irp->suggestion[--irp->suggestion_cnt];
  240. #ifdef MAIN
  241.             if(verbose) printf("suggested --  \n");
  242. #endif
  243.             }
  244.         else 
  245.             {
  246.                                                     RR(del, "del (zero)");
  247.             *xp = x1 + del;
  248. #ifdef MAIN
  249.             if(verbose) printf("step --  \n");
  250. #endif
  251.             }
  252.                                                     RR(*xp, "*xp (zero)");
  253.  
  254.         state = 2;
  255.         return OK;                /* get 2nd function value */
  256. loc2:
  257.  
  258.                                                     RR(value, "value (zero)");
  259.         if(value != f1) break;
  260.         del *= -1.5;        /* look on both sides of the starting point */
  261.         }
  262.                                 /* save 2nd function value */
  263.     x2 = *xp;
  264.     f2 = value;
  265.     d2 = fabs(f2);
  266.     if(f1*f2 <= 0.) straddling = 1;
  267.     else if(d1 < d2)
  268.         {            /* swap points 1 and 2 so f2 will be closer to target */
  269.         temp = x1; x1 = x2; x2 = temp;
  270.         temp = f1; f1 = f2; f2 = temp;
  271.         temp = d1; d1 = d2; d2 = temp;
  272.         }
  273.  
  274.     del = fabs(del);
  275.     difference = d2;
  276.     quadratic = 0;
  277.     for (i = 0; i < ITERATIONS; i++)    
  278.         {
  279.         if(irp->suggestion_cnt) 
  280.             {*xp = irp->suggestion[--irp->suggestion_cnt];
  281.             del = fabs(*xp - x2);
  282. #ifdef MAIN
  283.             if(verbose) printf("suggested --  \n");
  284. #endif
  285.             }
  286.         else if(quadratic)
  287.             {
  288. #ifdef MAIN
  289.             if(verbose) printf("quadratic --  \n");
  290. #endif
  291.             del = fabs(*xp - x2);
  292.             quadratic = 0;
  293.             }
  294.         else
  295.             {
  296. #ifdef MAIN
  297.             if(verbose) printf("   linear --  \n");
  298. #endif
  299.             if(x1 == x2)
  300.                 {return REPEATED_X;  /* FAIL - same x value twice in a row */
  301.                 }
  302.  
  303.             slope = (f2 - f1)/(x2 - x1);
  304.  
  305.             if(slope == 0.)
  306.                 {return NOT_CHANGING;    /* FAIL - function value not changing */
  307.                 }
  308.  
  309.             dx = -f2/slope;
  310.  
  311.             if(!straddling && fabs(dx) > extrapolation_limit*del)
  312.                 dx *= extrapolation_limit*del/fabs(dx);    /* limit step size */
  313.         
  314.             del = fabs(dx);
  315.             *xp = x2 + dx;
  316.             }
  317.  
  318.         state = 3;
  319.         return OK;                /* get 3rd or following function value */
  320. loc3:
  321.  
  322.         x3 = *xp;
  323.         f3 = value;
  324.         d3 = f3;
  325.         if(straddling)
  326.             {
  327. #ifndef LINEAR
  328.             if(quadratic = ((x3-x1)*(x2-x3) > 0. && (f3-f1)*(f2-f3) > 0.))
  329.                 {
  330.             /*
  331.                 x3 is strictly between x1 and x2 (preventing division
  332.                 by 0), and f3 is strictly between f1 and f2 (permitting
  333.                 interpolation).  Use Aitkin's Iteration method to find
  334.                 a polynomial in y that passes through the three most
  335.                 recent points (see Abramowitz & Stegun 25.2.23). 
  336.                 Evaluate the polynomial at the target value of y to
  337.                 find a candidate value for x.
  338.             */
  339. #define aitkins(x0, f0, x1, f1, x) (f0*(x1-x) - f1*(x0-x))/(x1-x0)
  340.                 x31 = aitkins(f3, x3, f1, x1, 0.);
  341.                 x32 = aitkins(f3, x3, f2, x2, 0.);
  342.                 x312 = aitkins(f1, x31, f2, x32, 0.);
  343.                             /* insist that the new x be between x1 and x2 */
  344.                 if((x312 - x1)*(x2 - x312) > 0.) *xp = x312;
  345.                 else quadratic = 0;
  346.                 }
  347. #endif
  348. #define CASE(num) if(verbose) printf("--  case %d  ", num)
  349.             if(d2*d3 <= 0.)
  350.                 {
  351.                         /* we have:              2
  352.                                                /
  353.                                     ---------------  <--  y = 0
  354.                                            /
  355.                                          3
  356.                                        /
  357.                                      1
  358.  
  359.                             ...so we discard point 1
  360.                         */
  361. #ifdef MAIN
  362. CASE(1);                /* straddling, new point same side as point 1 */
  363. #endif
  364.                 d1 = d3;
  365.                 x1 = x3;
  366.                 f1 = f3;
  367.                 }
  368.             else
  369.                 {
  370.                         /* we have:              2
  371.                                                /
  372.                                              3
  373.                                            /
  374.                                     ---------------  <--  y = 0
  375.                                        /
  376.                                      1    
  377.  
  378.                             ...so we discard point 2
  379.                         */
  380. #ifdef MAIN
  381. CASE(2);                /* straddling, new point same side as point 2 */
  382. #endif
  383.                 d2 = d3;
  384.                 x2 = x3;
  385.                 f2 = f3;
  386.                 }
  387.             difference = fabs(d3);
  388.             }
  389.         else if(d3*f2 <= 0.)
  390.             {straddling = 1;                /* just now straddling target */
  391. #ifdef MAIN
  392.             if(verbose) printf("-- straddling --\n");
  393. #endif
  394.             x31 = aitkins(f3, x3, f1, x1, 0.);
  395.             x32 = aitkins(f3, x3, f2, x2, 0.);
  396.             x312 = aitkins(f1, x31, f2, x32, 0.);
  397.  
  398.             /* insist that the new x be between x3 and x2 */
  399.             if(quadratic = ((x312 - x3)*(x2 - x312) > 0.)) 
  400.                 *xp = x312;
  401.             if(d2 < d1)
  402.                 {
  403.                         /* we have:              1
  404.                                                /
  405.                                              2
  406.                                            /
  407.                                     ---------------  <--  y = 0
  408.                                        /
  409.                                      3    
  410.  
  411.                             ...so we discard point 1
  412.                         */
  413. #ifdef MAIN
  414. CASE(3);                /* just now straddling */
  415. #endif
  416.                 d1 = d3;
  417.                 x1 = x3;
  418.                 f1 = f3;
  419.                 d2 = f2;
  420.                 difference = fabs(d3);
  421.                 }
  422.             else
  423.                 {
  424.                         /* we have:              2
  425.                                                /
  426.                                              1
  427.                                            /
  428.                                     ---------------  <--  y = 0
  429.                                        /
  430.                                      3    
  431.  
  432.                             ...so we discard point 2
  433.                         */
  434. #ifdef MAIN
  435. CASE(4);                /* just now straddling */
  436. #endif
  437.                 d2 = d3;
  438.                 x2 = x3;
  439.                 f2 = f3;
  440.                 d1 = f1;
  441.                 difference = fabs(d3);
  442.                 }
  443.             }
  444.         else
  445.             {
  446.             d3 = fabs(d3);
  447.             if(d3 < d2)
  448.                 {
  449.                         /* we have:              1
  450.                                                /
  451.                                              2
  452.                                            /
  453.                                          3    
  454.                                        /
  455.                                     ---------------  <--  y = 0
  456.  
  457.                             ...so we discard point 1
  458.                         */
  459. #ifdef MAIN
  460. CASE(5);                /* still extrapolating, new point is closer */
  461. #endif
  462.                             /* this extrapolation will speed 
  463.                                 convergence for a high order zero,
  464.                                 but has little effect for normal cases */
  465.                 if((x2-x3)*(x1-x2) > 0 && (f2-f3)*(f1-f2) > 0)
  466.                     {
  467.                     x31 = aitkins(f3, x3, f1, x1, 0.);
  468.                     x32 = aitkins(f3, x3, f2, x2, 0.);
  469.                     x312 = aitkins(f1, x31, f2, x32, 0.);
  470.                     if(quadratic = ((x3-x312)*(x2-x3) > 0)
  471.                                 && (fabs(x312-x3) < extrapolation_limit*del))
  472.                         *xp = x312;
  473.                     }
  474.  
  475.                 d1 = d2; d2 = d3;
  476.                 x1 = x2; x2 = x3;
  477.                 f1 = f2; f2 = f3;
  478.                 }
  479.             else if(d3 < d1)
  480.                 {
  481.                         /* we have:                     1
  482.                                        3            /
  483.                                          \          2
  484.                                            \____/
  485.                                     ---------------  <--  y = 0
  486.                             ...so we discard point 1 but look between 2 and 3
  487.                         */
  488. #ifdef MAIN
  489. CASE(6);                /* still extrapolating, new point between 1 and 2 */
  490. #endif
  491.                 if(i > 7)
  492.                     {return LOCAL_EXTREMUM;    /* FAIL - local extremum */
  493.                     }
  494.                 d1 = d3;
  495.                 x1 = x3;
  496.                 f1 = f3;
  497.                 *xp = (x2 + x3)/2.;
  498.                 quadratic = 1;
  499.                 }
  500.             else 
  501.                 {
  502.                         /* we have:       3
  503.                                          \                  1
  504.                                            \            /
  505.                                              \          2
  506.                                                \____/
  507.                                         ---------------  <--  y = 0
  508.                             ...so we look closer to point 2
  509.                         */
  510. #ifdef MAIN
  511. CASE(7);                /* still extrapolating, new point worse than previous */
  512. #endif
  513.                 if(i > 7)
  514.                     {return LOCAL_EXTREMUM;    /* FAIL - local extremum */
  515.                     }
  516.                 extrapolation_limit = .5;
  517.                 }
  518.             difference = d2;
  519.  
  520. #ifdef MAIN
  521.             assert (d2<=d1);
  522. #endif
  523.             }
  524.         }
  525.     return SLOW_CONVERGENCE;        /* FAIL - not converging fast enough */
  526. }
  527.  
  528.  
  529. #ifdef MAIN
  530.  
  531. show_history(irp) INVERSION_REC *irp;
  532. {
  533.     printf("x1 =%11.7f    f1 =%11.7f\n", x1, f1);
  534.     printf("x2 =%11.7f    f2 =%11.7f\n", x2, f2);
  535.     printf("x3 =%11.7f    f3 =%11.7f\n", x3, f3);
  536. }
  537.  
  538. #undef state
  539. #undef straddling
  540. #undef i
  541. #undef error
  542. #undef x1
  543. #undef x2
  544. #undef x3
  545. #undef f1
  546. #undef f2
  547. #undef f3
  548. #undef d1
  549. #undef d2
  550. #undef d3
  551. #undef del
  552. #undef difference
  553. #undef dx
  554. #undef extrapolation_limit
  555. #undef slope
  556. #undef xp
  557.  
  558. double fun(x) double x;
  559. {    double y;
  560.             /*    this is x - x^3, which has extrema f(-.57735) = -.3849 
  561.                                             and f( .57735) =  .3849 */
  562.     y = -(x+1.)*x*(x-1.);
  563.     return y;
  564. }
  565.  
  566. double targets[]={-.5, -.4, -.3849, -.3, -.2, -.1, .0, .1};
  567.  
  568. main(int argc, char **argv)
  569. {
  570.     double x, y, del, target, bits, dif;
  571.     int error, eval, it;
  572.     INVERSION_REC *inversion=NULL, *init_zero();
  573.  
  574.     if(argc > 1 && strcmp(argv[1], "-s")==0) verbose = 0;
  575.  
  576.  
  577.     printf("; searching...\n");
  578.     inversion = init_zero(&x, inversion);
  579.     advise_zero(-0.9, inversion);
  580. /*
  581.     y = fun(x=-.9);
  582.     show_progress(x, y, target);
  583.     error = zero(y - target, inversion);
  584. */
  585.     x = -1.0;
  586.     target = -.1;
  587.     do
  588.         {y = fun(x);
  589.         del = y - target;
  590.         if(fabs(del) < .00000001) break;
  591.         show_progress(x, y, target);
  592.         error = zero(y - target, inversion);
  593.         } while(!error);
  594.     y = fun(x);
  595.     printf(";             target = %11.7f\n", target);
  596.     printf(";              error = %11.7f\n", y-target);
  597.     printf(";                del = %11.7f\n", del);
  598.     printf(";             status = %11d  (%s)\n\n", error, 
  599.                                             error?"FAILURE":"success");
  600.  
  601.  
  602. /*
  603.     the sample function reaches a local minimum near -0.3849 for 
  604.     x = -0.57735, so the first two inversion attempts will fail.
  605.     for (target = -.5; target < .3; target += .1)
  606. */
  607.     for (it = 0; it < sizeof(targets)/sizeof(double); it++)
  608.         {target = targets[it];
  609.         x = -0.7;
  610.         if(verbose) printf("; searching...\n");
  611.         inversion = init_zero(&x, inversion);
  612.         eval = 0;
  613.         do
  614.             {y = fun(x); eval++;
  615.             del = y - target;
  616.             if(fabs(del) < .00000001) break;
  617.             if(verbose) show_progress(x, y, target);
  618.             error = zero(del, inversion);
  619.             } while(fabs(del) > .00000001 && !error);
  620.         y = fun(x);
  621.         if(verbose)
  622.             {if(error) printf("\n");
  623.             printf(";             target = %11.7f\n", target);
  624.             printf(";              error = %11.7f\n", y-target);
  625.             printf(";                del = %11.7f\n", del);
  626.             printf(";             status = %11d  (%s)\n\n", error, 
  627.                                             error?"FAILURE":"success");
  628.             }
  629.         else 
  630.             if(error)
  631.                 printf("   fun(   nothing  ) = %11.7f  in %2d evaluations\n", 
  632.                                                                 target, eval);
  633.             else
  634.                 printf("   fun(%11.7f ) = %11.7f  in %2d evaluations\n", 
  635.                                                                 x, y, eval);
  636.         }
  637.  
  638. }
  639.  
  640. show_progress(x, y, target) double x, y, target;
  641. {    double dif, bits;
  642.     
  643.     dif = y-target;
  644.     if(!target) bits = 0.;
  645.     else if(dif) bits = -log(fabs(dif/target))/log(2.);
  646.     else bits = 60.;
  647.                 printf("   fun(%11.7f ) = %11.7f  %5.2f bits ok\n", x, y, bits);
  648. }
  649.  
  650. #endif /* MAIN */
  651.  
  652.  
  653.