1 /*
2 * goal-seek.c: A generic root finder.
3 *
4 * Author:
5 * Morten Welinder (terra@gnome.org)
6 *
7 */
8
9 #undef DEBUG_GOAL_SEEK
10 #ifdef STANDALONE
11 #define DEBUG_GOAL_SEEK
12 #endif
13
14 #include <gnumeric-config.h>
15 #include <numbers.h>
16 #include <gnumeric.h>
17 #include <tools/goal-seek.h>
18 #include <gnm-random.h>
19 #include <value.h>
20 #include <cell.h>
21 #include <sheet.h>
22
23 #include <stdlib.h>
24 #include <math.h>
25 #include <limits.h>
26
27
28 static gboolean
update_data(gnm_float x,gnm_float y,GnmGoalSeekData * data)29 update_data (gnm_float x, gnm_float y, GnmGoalSeekData *data)
30 {
31 if (!gnm_finite (y))
32 return FALSE;
33
34 if (y > 0) {
35 if (data->havexpos) {
36 if (data->havexneg) {
37 /*
38 * When we have pos and neg, prefer the new point only
39 * if it makes the pos-neg x-interval smaller.
40 */
41 if (gnm_abs (x - data->xneg) < gnm_abs (data->xpos - data->xneg)) {
42 data->xpos = x;
43 data->ypos = y;
44 }
45 } else if (y < data->ypos) {
46 /* We have pos only and our neg y is closer to zero. */
47 data->xpos = x;
48 data->ypos = y;
49 }
50 } else {
51 data->xpos = x;
52 data->ypos = y;
53 data->havexpos = TRUE;
54 }
55 return FALSE;
56 } else if (y < 0) {
57 if (data->havexneg) {
58 if (data->havexpos) {
59 /*
60 * When we have pos and neg, prefer the new point only
61 * if it makes the pos-neg x-interval smaller.
62 */
63 if (gnm_abs (x - data->xpos) < gnm_abs (data->xpos - data->xneg)) {
64 data->xneg = x;
65 data->yneg = y;
66 }
67 } else if (-y < -data->yneg) {
68 /* We have neg only and our neg y is closer to zero. */
69 data->xneg = x;
70 data->yneg = y;
71 }
72
73 } else {
74 data->xneg = x;
75 data->yneg = y;
76 data->havexneg = TRUE;
77 }
78 return FALSE;
79 } else {
80 /* Lucky guess... */
81 data->have_root = TRUE;
82 data->root = x;
83 #ifdef DEBUG_GOAL_SEEK
84 g_print ("update_data: got root %.20" GNM_FORMAT_g "\n", x);
85 #endif
86 return TRUE;
87 }
88 }
89
90
91 /*
92 * Calculate a reasonable approximation to the derivative of a function
93 * in a single point.
94 */
95 static GnmGoalSeekStatus
fake_df(GnmGoalSeekFunction f,gnm_float x,gnm_float * dfx,gnm_float xstep,GnmGoalSeekData * data,void * user_data)96 fake_df (GnmGoalSeekFunction f, gnm_float x, gnm_float *dfx, gnm_float xstep,
97 GnmGoalSeekData *data, void *user_data)
98 {
99 gnm_float xl, xr, yl, yr;
100 GnmGoalSeekStatus status;
101
102 #ifdef DEBUG_GOAL_SEEK
103 g_print ("fake_df (x=%.20" GNM_FORMAT_g ", xstep=%.20" GNM_FORMAT_g ")\n",
104 x, xstep);
105 #endif
106
107 xl = x - xstep;
108 if (xl < data->xmin)
109 xl = x;
110
111 xr = x + xstep;
112 if (xr > data->xmax)
113 xr = x;
114
115 if (xl == xr) {
116 #ifdef DEBUG_GOAL_SEEK
117 g_print ("==> xl == xr\n");
118 #endif
119 return GOAL_SEEK_ERROR;
120 }
121
122 status = f (xl, &yl, user_data);
123 if (status != GOAL_SEEK_OK) {
124 #ifdef DEBUG_GOAL_SEEK
125 g_print ("==> failure at xl\n");
126 #endif
127 return status;
128 }
129 #ifdef DEBUG_GOAL_SEEK
130 g_print ("==> xl=%.20" GNM_FORMAT_g "; yl=%.20" GNM_FORMAT_g "\n",
131 xl, yl);
132 #endif
133
134 status = f (xr, &yr, user_data);
135 if (status != GOAL_SEEK_OK) {
136 #ifdef DEBUG_GOAL_SEEK
137 g_print ("==> failure at xr\n");
138 #endif
139 return status;
140 }
141 #ifdef DEBUG_GOAL_SEEK
142 g_print ("==> xr=%.20" GNM_FORMAT_g "; yr=%.20" GNM_FORMAT_g "\n",
143 xr, yr);
144 #endif
145
146 *dfx = (yr - yl) / (xr - xl);
147 #ifdef DEBUG_GOAL_SEEK
148 g_print ("==> %.20" GNM_FORMAT_g "\n", *dfx);
149 #endif
150 return gnm_finite (*dfx) ? GOAL_SEEK_OK : GOAL_SEEK_ERROR;
151 }
152
153 void
goal_seek_initialize(GnmGoalSeekData * data)154 goal_seek_initialize (GnmGoalSeekData *data)
155 {
156 data->havexpos = data->havexneg = data->have_root = FALSE;
157 data->xpos = data->xneg = data->root = gnm_nan;
158 data->ypos = data->yneg = gnm_nan;
159 data->xmin = -1e10;
160 data->xmax = +1e10;
161 data->precision = 1e-10;
162 }
163
164
165 /**
166 * goal_seek_point:
167 * @f: (scope call): object function
168 * @data: #GnmGoalSeekData state
169 * @user_data: user data for @f
170 * @x0: root guess
171 *
172 * Seek a goal using a single point.
173 *
174 * Returns:
175 */
176 GnmGoalSeekStatus
goal_seek_point(GnmGoalSeekFunction f,GnmGoalSeekData * data,void * user_data,gnm_float x0)177 goal_seek_point (GnmGoalSeekFunction f, GnmGoalSeekData *data,
178 void *user_data, gnm_float x0)
179 {
180 GnmGoalSeekStatus status;
181 gnm_float y0;
182
183 if (data->have_root)
184 return GOAL_SEEK_OK;
185
186 #ifdef DEBUG_GOAL_SEEK
187 g_print ("goal_seek_point\n");
188 #endif
189
190 if (x0 < data->xmin || x0 > data->xmax)
191 return GOAL_SEEK_ERROR;
192
193 status = f (x0, &y0, user_data);
194 if (status != GOAL_SEEK_OK)
195 return status;
196
197 if (update_data (x0, y0, data))
198 return GOAL_SEEK_OK;
199
200 return GOAL_SEEK_ERROR;
201 }
202
203
204 static GnmGoalSeekStatus
goal_seek_newton_polish(GnmGoalSeekFunction f,GnmGoalSeekFunction df,GnmGoalSeekData * data,void * user_data,gnm_float x0,gnm_float y0)205 goal_seek_newton_polish (GnmGoalSeekFunction f, GnmGoalSeekFunction df,
206 GnmGoalSeekData *data, void *user_data,
207 gnm_float x0, gnm_float y0)
208 {
209 int iterations;
210 gnm_float last_df0 = 1;
211 gboolean try_newton = TRUE;
212 gboolean try_square = x0 != 0 && gnm_abs (x0) < 1e10;
213
214 #ifdef DEBUG_GOAL_SEEK
215 g_print ("goal_seek_newton_polish\n");
216 #endif
217
218 for (iterations = 0; iterations < 20; iterations++) {
219 if (try_square) {
220 gnm_float x1 = x0 * gnm_abs (x0);
221 gnm_float y1, r;
222 GnmGoalSeekStatus status = f (x1, &y1, user_data);
223 if (status != GOAL_SEEK_OK)
224 goto nomore_square;
225
226 if (update_data (x1, y1, data))
227 return GOAL_SEEK_OK;
228
229 r = gnm_abs (y1 / y0);
230 if (r >= 1)
231 goto nomore_square;
232
233 x0 = x1;
234 #ifdef DEBUG_GOAL_SEEK
235 g_print ("polish square: x0=%.20" GNM_FORMAT_g "\n",
236 x0);
237 #endif
238 if (r > 0.5)
239 goto nomore_square;
240
241 continue;
242
243 nomore_square:
244 try_square = FALSE;
245 }
246
247 if (try_newton) {
248 gnm_float df0, r, x1, y1;
249 GnmGoalSeekStatus status = df
250 ? df (x0, &df0, user_data)
251 : fake_df (f, x0, &df0, gnm_abs (x0) / 1e6, data, user_data);
252 if (status != GOAL_SEEK_OK || df0 == 0)
253 df0 = last_df0; /* Bogus */
254 else
255 last_df0 = df0;
256
257 x1 = x0 - y0 / df0;
258 if (x1 < data->xmin || x1 > data->xmax)
259 goto nomore_newton;
260
261 status = f (x1, &y1, user_data);
262 if (status != GOAL_SEEK_OK)
263 goto nomore_newton;
264
265 if (update_data (x1, y1, data))
266 return GOAL_SEEK_OK;
267
268 r = gnm_abs (y1 / y0);
269 if (r >= 1)
270 goto nomore_newton;
271
272 x0 = x1;
273 #ifdef DEBUG_GOAL_SEEK
274 g_print ("polish Newton: x0=%.20" GNM_FORMAT_g "\n",
275 x0);
276 #endif
277 if (r > 0.5)
278 goto nomore_newton;
279
280 continue;
281
282 nomore_newton:
283 try_newton = FALSE;
284 }
285
286 /* Nothing left to try. */
287 break;
288 }
289
290 if (goal_seek_bisection (f, data, user_data) == GOAL_SEEK_OK)
291 return GOAL_SEEK_OK;
292
293 data->root = x0;
294 data->have_root = TRUE;
295 return GOAL_SEEK_OK;
296 }
297
298
299 /**
300 * goal_seek_newton:
301 * @f: (scope call): object function
302 * @df: (scope call) (nullable): object function derivative
303 * @data: #GnmGoalSeekData state
304 * @user_data: user data for @f and @df
305 * @x0: root guess
306 *
307 * Seek a goal (root) using Newton's iterative method.
308 *
309 * The supplied function must (should) be continuously differentiable in
310 * the supplied interval. If @df is %NULL, this function will
311 * estimate the derivative.
312 *
313 * This method will find a root rapidly provided the initial guess, x0,
314 * is sufficiently close to the root. (The number of significant digits
315 * (asymptotically) goes like i^2 unless the root is a multiple root in
316 * which case it is only like c*i.)
317 */
318 GnmGoalSeekStatus
goal_seek_newton(GnmGoalSeekFunction f,GnmGoalSeekFunction df,GnmGoalSeekData * data,void * user_data,gnm_float x0)319 goal_seek_newton (GnmGoalSeekFunction f, GnmGoalSeekFunction df,
320 GnmGoalSeekData *data, void *user_data, gnm_float x0)
321 {
322 int iterations;
323 gnm_float precision = data->precision / 2;
324 gnm_float last_df0 = 1;
325 gnm_float step_factor = 1e-6;
326
327 if (data->have_root)
328 return GOAL_SEEK_OK;
329
330 #ifdef DEBUG_GOAL_SEEK
331 g_print ("goal_seek_newton\n");
332 #endif
333
334 for (iterations = 0; iterations < 100; iterations++) {
335 gnm_float x1, y0, df0, stepsize;
336 GnmGoalSeekStatus status;
337 gboolean flat;
338
339 #ifdef DEBUG_GOAL_SEEK
340 g_print ("x0 = %.20" GNM_FORMAT_g " (i=%d)\n", x0, iterations);
341 #endif
342
343 /* Check whether we have left the valid interval. */
344 if (x0 < data->xmin || x0 > data->xmax)
345 return GOAL_SEEK_ERROR;
346
347 status = f (x0, &y0, user_data);
348 if (status != GOAL_SEEK_OK)
349 return status;
350
351 #ifdef DEBUG_GOAL_SEEK
352 g_print (" y0 = %.20" GNM_FORMAT_g "\n", y0);
353 #endif
354
355 if (update_data (x0, y0, data))
356 return GOAL_SEEK_OK;
357
358 if (df)
359 status = df (x0, &df0, user_data);
360 else {
361 gnm_float xstep;
362
363 if (gnm_abs (x0) < 1e-10)
364 if (data->havexneg && data->havexpos)
365 xstep = gnm_abs (data->xpos - data->xneg) / 1e6;
366 else
367 xstep = (data->xmax - data->xmin) / 1e6;
368 else
369 xstep = step_factor * gnm_abs (x0);
370
371 status = fake_df (f, x0, &df0, xstep, data, user_data);
372 }
373 if (status != GOAL_SEEK_OK)
374 return status;
375
376 /* If we hit a flat spot, we are in trouble. */
377 flat = (df0 == 0);
378 if (flat) {
379 last_df0 /= 2;
380 if (gnm_abs (last_df0) <= GNM_MIN)
381 return GOAL_SEEK_ERROR;
382 df0 = last_df0; /* Might be utterly bogus. */
383 } else
384 last_df0 = df0;
385
386 if (data->havexpos && data->havexneg)
387 x1 = x0 - y0 / df0;
388 else
389 /*
390 * Overshoot slightly to prevent us from staying on
391 * just one side of the root.
392 */
393 x1 = x0 - 1.000001 * y0 / df0;
394
395 stepsize = gnm_abs (x1 - x0) / (gnm_abs (x0) + gnm_abs (x1));
396
397 #ifdef DEBUG_GOAL_SEEK
398 g_print (" df0 = %.20" GNM_FORMAT_g "\n", df0);
399 g_print (" ss = %.20" GNM_FORMAT_g "\n", stepsize);
400 #endif
401
402 if (stepsize < precision) {
403 goal_seek_newton_polish (f, df, data, user_data, x0, y0);
404 return GOAL_SEEK_OK;
405 }
406
407 if (flat && iterations > 0) {
408 /*
409 * Verify that we made progress using our
410 * potentially bogus df0.
411 */
412 gnm_float y1;
413
414 if (x1 < data->xmin || x1 > data->xmax)
415 return GOAL_SEEK_ERROR;
416
417 status = f (x1, &y1, user_data);
418 if (status != GOAL_SEEK_OK)
419 return status;
420
421 #ifdef DEBUG_GOAL_SEEK
422 g_print (" y1 = %.20" GNM_FORMAT_g "\n", y1);
423 #endif
424 if (gnm_abs (y1) >= 0.9 * gnm_abs (y0))
425 return GOAL_SEEK_ERROR;
426 }
427
428 if (stepsize < step_factor)
429 step_factor = stepsize;
430
431 x0 = x1;
432 }
433
434 return GOAL_SEEK_ERROR;
435 }
436
437 /**
438 * goal_seek_bisection:
439 * @f: (scope call): object function
440 * @data: #GnmGoalSeekData state.
441 * @user_data: user data for @f.
442 *
443 * Seek a goal (root) using bisection methods.
444 *
445 * The supplied function must (should) be continuous over the interval.
446 *
447 * Caller must have located a positive and a negative point.
448 *
449 * This method will find a root steadily using bisection to narrow the
450 * interval in which a root lies.
451 *
452 * It alternates between mid-point-bisection (semi-slow, but guaranteed
453 * progress), secant-bisection (usually quite fast, but sometimes gets
454 * nowhere), and Ridder's Method (usually fast, harder to fool than
455 * the secant method).
456 */
457 GnmGoalSeekStatus
goal_seek_bisection(GnmGoalSeekFunction f,GnmGoalSeekData * data,void * user_data)458 goal_seek_bisection (GnmGoalSeekFunction f, GnmGoalSeekData *data,
459 void *user_data)
460 {
461 int iterations;
462 gnm_float stepsize;
463 int newton_submethod = 0;
464
465 if (data->have_root)
466 return GOAL_SEEK_OK;
467
468 #ifdef DEBUG_GOAL_SEEK
469 g_print ("goal_seek_bisection\n");
470 #endif
471
472 if (!data->havexpos || !data->havexneg)
473 return GOAL_SEEK_ERROR;
474
475 stepsize = gnm_abs (data->xpos - data->xneg)
476 / (gnm_abs (data->xpos) + gnm_abs (data->xneg));
477
478 /* log_2 (10) = 3.3219 < 4. */
479 for (iterations = 0; iterations < 100 + GNM_DIG * 4; iterations++) {
480 gnm_float xmid, ymid;
481 GnmGoalSeekStatus status;
482 enum { M_SECANT, M_RIDDER, M_NEWTON, M_MIDPOINT } method;
483
484 method = (iterations % 4 == 0)
485 ? M_RIDDER
486 : ((iterations % 4 == 2)
487 ? M_NEWTON
488 : M_MIDPOINT);
489
490 again:
491 switch (method) {
492 default:
493 abort ();
494
495 case M_SECANT:
496 xmid = data->xpos - data->ypos *
497 ((data->xneg - data->xpos) /
498 (data->yneg - data->ypos));
499 break;
500
501 case M_RIDDER: {
502 gnm_float det;
503
504 xmid = (data->xpos + data->xneg) / 2;
505 status = f (xmid, &ymid, user_data);
506 if (status != GOAL_SEEK_OK)
507 continue;
508 if (ymid == 0) {
509 update_data (xmid, ymid, data);
510 return GOAL_SEEK_OK;
511 }
512
513 det = gnm_sqrt (ymid * ymid - data->ypos * data->yneg);
514 if (det == 0)
515 /* This might happen with underflow, I guess. */
516 continue;
517
518 xmid += (xmid - data->xpos) * ymid / det;
519 break;
520 }
521
522 case M_MIDPOINT:
523 xmid = (data->xpos + data->xneg) / 2;
524 break;
525
526 case M_NEWTON: {
527 gnm_float x0, y0, xstep, df0;
528
529 /* This method is only effective close-in. */
530 if (stepsize > 0.1) {
531 method = M_MIDPOINT;
532 goto again;
533 }
534
535 switch (newton_submethod++ % 4) {
536 case 0: x0 = data->xpos; y0 = data->ypos; break;
537 case 2: x0 = data->xneg; y0 = data->yneg; break;
538 default:
539 case 3:
540 case 1:
541 x0 = (data->xpos + data->xneg) / 2;
542
543 status = f (x0, &y0, user_data);
544 if (status != GOAL_SEEK_OK)
545 continue;
546 }
547
548 xstep = gnm_abs (data->xpos - data->xneg) / 1e6;
549 status = fake_df (f, x0, &df0, xstep, data, user_data);
550 if (status != GOAL_SEEK_OK)
551 continue;
552
553 if (df0 == 0)
554 continue;
555
556 /*
557 * Overshoot by 1% to prevent us from staying on
558 * just one side of the root.
559 */
560 xmid = x0 - 1.01 * y0 / df0;
561 }
562 }
563
564 if ((xmid < data->xpos && xmid < data->xneg) ||
565 (xmid > data->xpos && xmid > data->xneg)) {
566 /* We left the interval. */
567 xmid = (data->xpos + data->xneg) / 2;
568 method = M_MIDPOINT;
569 }
570
571 status = f (xmid, &ymid, user_data);
572 if (status != GOAL_SEEK_OK)
573 continue;
574
575 #ifdef DEBUG_GOAL_SEEK
576 {
577 const char *themethod;
578 switch (method) {
579 case M_MIDPOINT: themethod = "midpoint"; break;
580 case M_RIDDER: themethod = "Ridder"; break;
581 case M_SECANT: themethod = "secant"; break;
582 case M_NEWTON: themethod = "Newton"; break;
583 default: themethod = "?";
584 }
585
586 g_print ("xmid = %.20" GNM_FORMAT_g " (%s)\n", xmid, themethod);
587 g_print (" ymid = %.20" GNM_FORMAT_g "\n", ymid);
588 }
589 #endif
590
591 if (update_data (xmid, ymid, data)) {
592 return GOAL_SEEK_OK;
593 }
594
595 stepsize = gnm_abs (data->xpos - data->xneg)
596 / (gnm_abs (data->xpos) + gnm_abs (data->xneg));
597
598 #ifdef DEBUG_GOAL_SEEK
599 g_print (" ss = %.20" GNM_FORMAT_g "\n", stepsize);
600 #endif
601
602 if (stepsize < GNM_EPSILON) {
603 if (data->yneg < ymid)
604 ymid = data->yneg, xmid = data->xneg;
605
606 if (data->ypos < ymid)
607 ymid = data->ypos, xmid = data->xpos;
608
609 data->have_root = TRUE;
610 data->root = xmid;
611 return GOAL_SEEK_OK;
612 }
613 }
614 return GOAL_SEEK_ERROR;
615 }
616
617 #undef SECANT_P
618 #undef RIDDER_P
619
620 /**
621 * goal_seek_trawl_uniformly:
622 * @f: (scope call): object function
623 * @data: #GnmGoalSeekData state
624 * @user_data: user data for @f
625 * @xmin: lower search bound
626 * @xmax: upper search bound
627 * @points: number of points to try.
628 */
629 GnmGoalSeekStatus
goal_seek_trawl_uniformly(GnmGoalSeekFunction f,GnmGoalSeekData * data,void * user_data,gnm_float xmin,gnm_float xmax,int points)630 goal_seek_trawl_uniformly (GnmGoalSeekFunction f,
631 GnmGoalSeekData *data, void *user_data,
632 gnm_float xmin, gnm_float xmax,
633 int points)
634 {
635 int i;
636
637 if (data->have_root)
638 return GOAL_SEEK_OK;
639
640 #ifdef DEBUG_GOAL_SEEK
641 g_print ("goal_seek_trawl_uniformly\n");
642 #endif
643
644 if (xmin > xmax || xmin < data->xmin || xmax > data->xmax)
645 return GOAL_SEEK_ERROR;
646
647 for (i = 0; i < points; i++) {
648 gnm_float x, y;
649 GnmGoalSeekStatus status;
650
651 if (data->havexpos && data->havexneg)
652 break;
653
654 x = xmin + (xmax - xmin) * random_01 ();
655 status = f (x, &y, user_data);
656 if (status != GOAL_SEEK_OK)
657 /* We are not depending on the result, so go on. */
658 continue;
659
660 #ifdef DEBUG_GOAL_SEEK
661 g_print ("x = %.20" GNM_FORMAT_g "\n", x);
662 g_print (" y = %.20" GNM_FORMAT_g "\n", y);
663 #endif
664
665 if (update_data (x, y, data))
666 return GOAL_SEEK_OK;
667 }
668
669 /* We were not (extremely) lucky, so we did not actually hit the
670 root. We report this as an error. */
671 return GOAL_SEEK_ERROR;
672 }
673
674 /**
675 * goal_seek_trawl_normally:
676 * @f: (scope call): object function
677 * @data: #GnmGoalSeekData state
678 * @user_data: user data for @f
679 * @mu: search mean
680 * @sigma: search standard deviation
681 * @points: number of points to try.
682 */
683 GnmGoalSeekStatus
goal_seek_trawl_normally(GnmGoalSeekFunction f,GnmGoalSeekData * data,void * user_data,gnm_float mu,gnm_float sigma,int points)684 goal_seek_trawl_normally (GnmGoalSeekFunction f,
685 GnmGoalSeekData *data, void *user_data,
686 gnm_float mu, gnm_float sigma,
687 int points)
688 {
689 int i;
690
691 if (data->have_root)
692 return GOAL_SEEK_OK;
693
694 #ifdef DEBUG_GOAL_SEEK
695 g_print ("goal_seek_trawl_normally\n");
696 #endif
697
698 if (sigma <= 0 || mu < data->xmin || mu > data->xmax)
699 return GOAL_SEEK_ERROR;
700
701 for (i = 0; i < points; i++) {
702 gnm_float x, y;
703 GnmGoalSeekStatus status;
704
705 if (data->havexpos && data->havexneg)
706 break;
707
708 x = mu + sigma * random_normal ();
709 if (x < data->xmin || x > data->xmax)
710 continue;
711
712 status = f (x, &y, user_data);
713 if (status != GOAL_SEEK_OK)
714 /* We are not depending on the result, so go on. */
715 continue;
716
717 #ifdef DEBUG_GOAL_SEEK
718 g_print ("x = %.20" GNM_FORMAT_g "\n", x);
719 g_print (" y = %.20" GNM_FORMAT_g "\n", y);
720 #endif
721
722 if (update_data (x, y, data))
723 return GOAL_SEEK_OK;
724 }
725
726 /* We were not (extremely) lucky, so we did not actually hit the
727 root. We report this as an error. */
728 return GOAL_SEEK_ERROR;
729 }
730
731 /**
732 * gnm_goal_seek_eval_cell:
733 * @x: x-value for which to evaluate
734 * @y: (out): location to store result
735 * @data: user data
736 *
737 * Returns: An status indicating whether evaluation went ok.
738 */
739 GnmGoalSeekStatus
gnm_goal_seek_eval_cell(gnm_float x,gnm_float * y,gpointer data_)740 gnm_goal_seek_eval_cell (gnm_float x, gnm_float *y, gpointer data_)
741 {
742 GnmGoalSeekCellData const *data = data_;
743 GnmValue *v = value_new_float (x);
744
745 gnm_cell_set_value (data->xcell, v);
746 cell_queue_recalc (data->xcell);
747 gnm_cell_eval (data->ycell);
748
749 if (data->ycell->value &&
750 VALUE_IS_NUMBER (data->ycell->value)) {
751 *y = value_get_as_float (data->ycell->value) - data->ytarget;
752 if (gnm_finite (*y))
753 return GOAL_SEEK_OK;
754 }
755
756 return GOAL_SEEK_ERROR;
757 }
758
759 GnmGoalSeekStatus
gnm_goal_seek_cell(GnmGoalSeekData * data,GnmGoalSeekCellData * celldata)760 gnm_goal_seek_cell (GnmGoalSeekData *data,
761 GnmGoalSeekCellData *celldata)
762 {
763 GnmGoalSeekStatus status;
764 gboolean hadold;
765 gnm_float oldx;
766 GnmValue *v;
767
768 hadold = !VALUE_IS_EMPTY_OR_ERROR (celldata->xcell->value);
769 oldx = hadold ? value_get_as_float (celldata->xcell->value) : 0;
770
771 /* PLAN A: Newton's iterative method from initial or midpoint. */
772 {
773 gnm_float x0;
774
775 if (hadold && oldx >= data->xmin && oldx <= data->xmax)
776 x0 = oldx;
777 else
778 x0 = (data->xmin + data->xmax) / 2;
779
780 status = goal_seek_newton (gnm_goal_seek_eval_cell, NULL,
781 data, celldata,
782 x0);
783 if (status == GOAL_SEEK_OK)
784 goto DONE;
785 }
786
787 /* PLAN B: Trawl uniformly. */
788 if (!data->havexpos || !data->havexneg) {
789 status = goal_seek_trawl_uniformly (gnm_goal_seek_eval_cell,
790 data, celldata,
791 data->xmin, data->xmax,
792 100);
793 if (status == GOAL_SEEK_OK)
794 goto DONE;
795 }
796
797 /* PLAN C: Trawl normally from middle. */
798 if (!data->havexpos || !data->havexneg) {
799 gnm_float sigma, mu;
800 int i;
801
802 sigma = MIN (data->xmax - data->xmin, 1e6);
803 mu = (data->xmax + data->xmin) / 2;
804
805 for (i = 0; i < 5; i++) {
806 sigma /= 10;
807 status = goal_seek_trawl_normally (gnm_goal_seek_eval_cell,
808 data, celldata,
809 mu, sigma, 30);
810 if (status == GOAL_SEEK_OK)
811 goto DONE;
812 }
813 }
814
815 /* PLAN D: Trawl normally from left. */
816 if (!data->havexpos || !data->havexneg) {
817 gnm_float sigma, mu;
818 int i;
819
820 sigma = MIN (data->xmax - data->xmin, 1e6);
821 mu = data->xmin;
822
823 for (i = 0; i < 5; i++) {
824 sigma /= 10;
825 status = goal_seek_trawl_normally (gnm_goal_seek_eval_cell,
826 data, celldata,
827 mu, sigma, 20);
828 if (status == GOAL_SEEK_OK)
829 goto DONE;
830 }
831 }
832
833 /* PLAN E: Trawl normally from right. */
834 if (!data->havexpos || !data->havexneg) {
835 gnm_float sigma, mu;
836 int i;
837
838 sigma = MIN (data->xmax - data->xmin, 1e6);
839 mu = data->xmax;
840
841 for (i = 0; i < 5; i++) {
842 sigma /= 10;
843 status = goal_seek_trawl_normally (gnm_goal_seek_eval_cell,
844 data, celldata,
845 mu, sigma, 20);
846 if (status == GOAL_SEEK_OK)
847 goto DONE;
848 }
849 }
850
851 /* PLAN F: Newton iteration with uniform net of starting points. */
852 if (!data->havexpos || !data->havexneg) {
853 int i;
854 const int N = 10;
855
856 for (i = 1; i <= N; i++) {
857 gnm_float x0 = data->xmin +
858 (data->xmax - data->xmin) / (N + 1) * i;
859
860 status = goal_seek_newton (gnm_goal_seek_eval_cell, NULL,
861 data, celldata,
862 x0);
863 if (status == GOAL_SEEK_OK)
864 goto DONE;
865 }
866 }
867
868 /* PLAN Z: Bisection. */
869 {
870 status = goal_seek_bisection (gnm_goal_seek_eval_cell,
871 data, celldata);
872 if (status == GOAL_SEEK_OK)
873 goto DONE;
874 }
875
876 DONE:
877 if (status == GOAL_SEEK_OK)
878 v = value_new_float (data->root);
879 else if (hadold)
880 v = value_new_float (oldx);
881 else
882 v = value_new_empty ();
883 sheet_cell_set_value (celldata->xcell, v);
884
885 return status;
886 }
887