489{
490 double a[MAXROUNDS][MAXROUNDS];
491
492
493 for (int k = 0; k < x.dim; k++)
494 {
495 double left, right;
496 double save_x = x[k];
497 double h = h_start;
498
499
500 x[k] = save_x - h;
501 left = Evaluate(x);
502
503
504 if (k == 0 || left < dfmin)
505 dfmin = left, dpmin = x;
506
507
508 x[k] = save_x + h;
509 right = Evaluate(x);
510
511
512 if (right < dfmin)
513 dfmin = left, dpmin = x;
514
515
516 a[0][0] = (right - left) / (2.0 * h);
517
518
519 double err = 1e30;
520
521
522
523 for (int i = 1; i < MAXROUNDS; i++)
524 {
525
526 h *= SQRT_HALF;
527
528
529 x[k] = save_x - h;
530 left = Evaluate(x);
531 if (left < dfmin) dfmin = left, dpmin = x;
532 x[k] = save_x + h;
533 right = Evaluate(x);
534 if (right < dfmin) dfmin = right, dpmin = x;
535
536
537 a[0][i] = (right - left) / (2.0 * h);
538
539
540 double factor = TWO;
541
542 for (int j = 1; j <= i; j++)
543 {
544 a[j][i] = (a[j-1][i] * factor - a[j-1][i-1])/(factor - 1.0);
545
546 factor *= TWO;
547
548 double error = max(fabs(a[j][i] - a[j-1][i]), fabs(a[j][i] - a[j-1][i-1]));
549
550
551 if (error < err)
552 {
553 err = error;
554 d[k] = a[j][i];
555 }
556 }
557
558
559 if (fabs(a[i][i] - a[i-1][i-1]) >= 2.0 * err)
560 {
561 x[k] = save_x;
562 break;
563 }
564 }
565
566 x[k] = save_x;
567 }
568}