[ create a new paste ] login | about

Link: http://codepad.org/vtSpYtVV    [ raw code | output | fork | 1 comment ]

mohit_at_codepad - C++, pasted on Jun 1:
// Compare babylonian and newton raphson
#include <stdio.h>

int iter = 0;

double babylonian_sqrt(double arg) {
  double eps = 1E-16;
  double rt = arg / 4.0;
  double sq = rt * rt;
  while( (arg - sq > eps)  || (sq - arg > eps) ) {
    rt += arg / rt;
    rt /= 2;
    sq = rt * rt;
    ++iter;
  }
  return rt;
}

double nr_sqrt(double arg) {
  double eps = 1E-16;
  double rt = arg / 4.0;
  double sq = rt * rt;
  while( (arg - sq > eps)  || (sq - arg > eps) ) {
    double fx = sq - arg;
    double dx = 2 * rt;
    rt -= fx / dx;
    sq = rt * rt;
    ++iter;
  }
  return rt;
}

int main() {
  const double num = 63.;
  iter = 0;
  const double baby = babylonian_sqrt(num);
  printf( "Sqrt(%g) by Ancient Baylonian (After %03d iteration) = %.20g\n"
               , num
               , iter
               , baby);
  iter = 0;
  const double nr = nr_sqrt(num);
  printf( "Sqrt(%g) by Newton - Raphson  (After %03d iteration) = %.20g\n"
               , num
               , iter
               , nr);
  return 0;
}


Output:
1
2
Sqrt(63) by Ancient Baylonian (After 006 iteration) = 7.9372539331937721485
Sqrt(63) by Newton - Raphson  (After 006 iteration) = 7.9372539331937721485


Create a new paste based on this one


Comments:
posted by mohit_at_codepad on Jun 2
Babylonian =>
x1 = 0.5 (X0 + NUM/X0)
= X0 + 0.5 (-X0 + NUM/X0)
= 0.5 (-X0 ^ 2 + NUM) / X0
= (NUM - X0 ^ 2) / (2 * X0)
<= Newton raphson
Both are same. You can easily drive babylonian method from Newton raphson method.
reply