{

    # saved scalars
    my ( $check, $fa, $fb, $ff, $ifconv, $iter, $maxit, $sa, $sb, $tol, $xx );

    # saved arrays
    my (@bvec);

    BEGIN {

        #
        #       Use Brent's procedure ZERO to find a root
        #       The function fofx contains the function whose root is sought
        #       This version looks for the root of sin(x)+x-1=0
        $tol   = 0.000001;
        $maxit = 25;
    }

    #
    #       find bounds for the root
    $sa = 0;
    $fa = fofx( \$sa );
    $sb = 1;
    $fb = fofx( \$sb );

    #
    #      normally we would check that the root is bounded by sa and sb,
    #      and keep looking for a bound if necessary.  But in this case
    #      the above values are known to bound the root.
    if ( $fa * $fb >= 0 ) {
        print 'sa=', '  ', $sa, '  ', 'fa=', '  ', $fa, '  ', 'sb=', '  ', $sb,
          '  ', 'fb=', '  ', $fb, "\n";
        print STDERR " STOP  'Bad starting states' statement executed\n";
        exit;
    }

    #     initialize for a search
    nbrenti( \$sa, \$fa, \$sb, \$fb, \$tol, \$xx, [ \@bvec, 0 ] );

    #     loop over iterations
    for ( $iter = 1 ; $iter <= $maxit ; $iter += 1 ) {

        $ff = fofx( \$xx );

        #     get next value of xx
        nbrentx( \$ff, \$xx, \$ifconv, [ \@bvec, 0 ] );
        if ( $ifconv != 0 ) {
            goto L700;
        }

        #       note - it is possible to make a convergence test here on the
        #       magnitude of ff as well

    }

    #     no convergence after maxit iterations -- cant happen if maxit
    #               is large enough for the given tol - see brents book

    #     converged - finish up
  L700:
    $check = fofx( \$xx );
    print 'last x is ', '  ', $xx, '  ', ' f(x)= ', '  ', $check, "\n";
    print STDERR " STOP  '**BYE**' statement executed\n";
    exit;
}

{    # sub fofx

    # unsaved scalars

    #
    #     This is the function whose root we seek.
    #     This version looks for the root of x+sin(x)=1
    sub fofx {

        # call arg declarations
        my ($R_x) = @_;
        use vars qw($x);

        local (*x) = ($R_x);
        my ($fofx);
        $fofx = $x + sin($x) - 1;
        return $fofx;
    }
}

{    # sub nbrenti

    # saved scalars
    my ( $fc, $ifconv, $sc );

    sub nbrenti {

        # call arg declarations
        my ( $R_sa, $R_fa, $R_sb, $R_fb, $R_tol, $R_xnext, $R_bvec ) = @_;
        use vars qw($sa $fa $sb $fb $tol $xnext);

        local ( *sa, *fa, *sb, *fb, *tol, *xnext ) =
          ( $R_sa, $R_fa, $R_sb, $R_fb, $R_tol, $R_xnext );
        my ( $ARRAY_REF_bvec, $OFFSET_bvec ) = @{$R_bvec};

        my $bvec = sub  : lvalue {
            $ARRAY_REF_bvec->[ $OFFSET_bvec + $_[0] - 1 ];
        };

        use vars qw($C_bvec);

        #
        #       initialization for brents root finder
        #       (see-algorithms for minimization without
        #       derivatives-,by richard p. brent, prentice-hall, 1973.)
        #
        #       input parameters -
        #       bounding states:
        #         sa - first x coordinate
        #         fa = f(sa)
        #         sb - second x coordinate
        #         fb = f(sb)
        #         **  must have fa*fb <= 0. **
        #       tolx = convergence tolerance on x
        #
        #       output parameter -
        #       xnext = next value of x to try
        #
        #       there are two subroutines:
        #               brenti must be called once to start an interation
        #               brentx must be called once for each step
        #
        #****************************************************************
        #
        #       sample coding: assume we have a root trapped between sa & sb
        #
        #       dimension bvec(7)
        #       data tol/.001/,maxit/25/
        #
        #     call nbrenti (sa,fa,sb,fb,tol,xx,bvec)
        #
        #     loop over iterations
        #     do 250 iter=1,maxit
        #
        #       calculate ff(xx)
        #
        #       get next value of xx
        #     call nbrentx (ff,xx,ifconv,bvec)
        #     if (ifconv.ne.0) go to 700
        #
        #       note - it is possible to make a convergence test here on the
        #       magnitude of ff as well
        #
        # 250 continue
        #
        #     no convergence after maxit iterations -- cant happen if maxit
        #               is large enough for the given tol - see brents book
        #
        #     converged - finish up
        # 700 ...
        #****************************************************************
        #
        $sc        = $sb;
        $fc        = $fb;
        $bvec->(1) = $sa;
        $bvec->(2) = $fa;
        $bvec->(3) = $sb;
        $bvec->(4) = $fb;
        $bvec->(5) = $sc;
        $bvec->(6) = $fc;
        $bvec->(7) = $tol;
        nbrentx( $R_fb, $R_xnext, \$ifconv, $R_bvec );
        return;
    }
}

{    # sub nbrentx

    # saved scalars
    my ( $bd, $be, $bm, $bp, $bq, $br, $bs, $fa, $fb, $fc, $sa, $sb, $sc,
        $tol );

    sub nbrentx {

        # call arg declarations
        my ( $R_flast, $R_xnext, $R_ifconv, $R_bvec ) = @_;
        use vars qw($flast $xnext $ifconv);

        local ( *flast, *xnext, *ifconv ) = ( $R_flast, $R_xnext, $R_ifconv );
        my ( $ARRAY_REF_bvec, $OFFSET_bvec ) = @{$R_bvec};

        my $bvec = sub  : lvalue {
            $ARRAY_REF_bvec->[ $OFFSET_bvec + $_[0] - 1 ];
        };

        use vars qw($C_bvec);

        #
        #     take one iteration step with brents method
        #     flast = value of f of previous step
        #     xnext = next value of x to try
        #     ifconv = 0 if not converged
        #            = 1 if converged
        #       bvec = state vector
        #
        #
        #       get the previous state
        $sa  = $bvec->(1);
        $fa  = $bvec->(2);
        $sb  = $bvec->(3);
        $fb  = $bvec->(4);
        $sc  = $bvec->(5);
        $fc  = $bvec->(6);
        $tol = $bvec->(7);

        #
        $fb = $flast;
        if ( ( $fb > 0. ) && ( $fc > 0. ) ) {
            goto L210;
        }
        if ( ( $fb <= 0. ) && ( $fc <= 0. ) ) {
            goto L210;
        }
        goto L220;

        #
      L210:
        $sc = $sa;
        $fc = $fa;
        $be = $sb - $sa;
        $bd = $be;
      L220:
        if ( abs($fc) >= abs($fb) ) {
            goto L230;
        }
        $sa = $sb;
        $sb = $sc;
        $sc = $sa;
        $fa = $fb;
        $fb = $fc;
        $fc = $fa;
      L230:
        $bm = 0.5 * ( $sc - $sb );

        if ( ( abs($bm) <= $tol ) || ( $fb == 0. ) ) {
            goto L340;
        }
        if ( ( abs($be) >= $tol ) && ( abs($fa) > abs($fb) ) ) {
            goto L240;
        }

        #
        #         bisection is forced
        $be = $bm;
        $bd = $be;
        goto L300;

        #
      L240:
        $bs = $fb / $fa;
        if ( $sa != $sc ) {
            goto L250;
        }

        #
        #         linear interpolation
        $bp = 2. * $bm * $bs;
        $bq = 1. - $bs;
        goto L260;

        #
        #         inverse quadratic interpolation
      L250:
        $bq = $fa / $fc;
        $br = $fb / $fc;
        $bp = $bs *
          ( 2. * $bm * $bq * ( $bq - $br ) - ( $sb - $sa ) * ( $br - 1. ) );
        $bq = ( $bq - 1. ) * ( $br - 1. ) * ( $bs - 1. );

        #
      L260:
        if ( $bp <= 0. ) {
            goto L270;
        }
        $bq = -$bq;
        goto L280;
      L270:
        $bp = -$bp;
      L280:
        $bs = $be;
        $be = $bd;

        if (   ( 2. * $bp >= 3. * $bm * $bq - abs( $tol * $bq ) )
            || ( $bp >= abs( 0.5 * $bs * $bq ) ) )
        {
            goto L290;
        }
        $bd = $bp / $bq;
        goto L300;
      L290:
        $be = $bm;
        $bd = $be;
      L300:
        $sa = $sb;
        $fa = $fb;

        if ( abs($bd) <= $tol ) {
            goto L310;
        }
        $sb = $sb + $bd;
        goto L330;
      L310:
        if ( $bm <= 0. ) {
            goto L320;
        }
        $sb = $sb + $tol;
        goto L330;
      L320:
        $sb = $sb - $tol;
      L330:
        $xnext  = $sb;
        $ifconv = 0;
        goto L900;
      L340:
        $ifconv = 1;

        #
        #       save the state
      L900:
        $bvec->(1) = $sa;
        $bvec->(2) = $fa;
        $bvec->(3) = $sb;
        $bvec->(4) = $fb;
        $bvec->(5) = $sc;
        $bvec->(6) = $fc;
        return;
    }
}