/* File: sqrt5.c

Calculate the square root by successive approximation.

Let's put our square root algorithm into a function.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

double mySqrt(double x, double epsilon) ;


int main() {
   double input, error ;

   printf("This program computes the square root of x.\n") ;
   printf("x =? ") ;
   scanf("%lf", &input) ;
   printf("error =? ") ;
   scanf("%lf", &error) ;


   printf("\n\nThe square root of %.8f is approximately %.8f (within %.8f)\n",
     input,  mySqrt(input, error), error) ;

   printf("The \"real\" answer is: %.8f\n", sqrt(input) ) ;

   return 0 ;
}


/* Function: 

   mySqrt(x, epsilon) computes the square root of x within
   an error of epsilon.

*/

double mySqrt(double x, double epsilon) {

   double highEnd, lowEnd, newGuess ;

   if (x < 0) {
      printf("Negative numbers do not have a square root!\n") ;
      abort() ;   // quit program
   }

   if (epsilon <= 0) {
      printf("Epsilon must be positive!!\n") ;
      abort() ;   // quit program
   }

   // Remember that for x,  0 < x < 1,  the sqrt(x) > x.
   if (x >= 1) {
      lowEnd = 1.0 ;
      highEnd = x ;
   } else {
      highEnd = 1.0 ;
      lowEnd = x ;
   }

   // invariant: lowEnd * lowEnd <= x && x <= highEnd * highEnd

   while ( highEnd - lowEnd > epsilon) {

      // printf("\nhighEnd = %.8f, lowEnd = %.8f\n", highEnd, lowEnd) ;
      assert ( lowEnd * lowEnd <= x && x <= highEnd * highEnd ) ;

      newGuess = (highEnd + lowEnd) / 2.0 ; 

      if (newGuess * newGuess == x) {
         // printf("%.8f is the answer! lucky guess!\n", newGuess) ;
         highEnd = newGuess ;
         lowEnd = newGuess ;

      } else {
         if (newGuess * newGuess > x) {
            // printf("%.8f is too big\n", newGuess) ;
            highEnd = newGuess ;
         } else {
            // printf("%.8f is too low\n", newGuess) ;
            lowEnd = newGuess ;
         }
      }
   
   }

   return (highEnd + lowEnd) / 2.0  ;
}