I believe this is based on an old Greek method...

  #define DECK 11

  double z() {
    const double s=2.0, c=0.1;
    double arr[DECK], low=1.0, up=s, inc=c;
    int j, k;
    for(j=0; j<15; ++j) {
      for(k=0; k<DECK; ++k) {
        *(arr+k)=low+(double)k*inc;
      }
      for(k=1; k<DECK; ++k) {
        if( *(arr+k-1)*(*(arr+k-1)) <s && *(arr+k)*(*(arr+k)) >s ) {
          low=*(arr+k-1);
          up=*(arr+k);
          break;
        }
      }
      inc*=c;
    }
    if( s-low*low < up*up-s )
      return low;
    else
      return up;
  }

This particular example only works to fifteen decimal places (roughly the limit on the precision of a double).