ニュートン法

/*
newton.c: ニュートン法
*/
#include <stdio.h> // printf, fprintf, fgets, sscanf
#include <math.h> // fabs

double newton(double a,double eps)
{
    int n = 0;
    double x0, err, x;
    x = (a + 1.0) / 2.0;
    printf("# n, x, err\n");
    printf("%4d, %.15e\n", n, x);
    do {
        n++;
        x0 = x;
        x = 1.0 / 2 * (x0 + a / x0);
        err = fabs(x - x0);
        printf("%4d, %.15e, %.15e\n", n, x, err);
    } while (err >= eps);
    return x;
}

int main(void)
{
    int n = 0;
    double a = 0, x, x0, err, eps = 1.0e-10;
    char s[128];

    fprintf(stderr, " a = "); fgets(s, 128, stdin); sscanf(s, "%lf", &a);
    while (a <= 0.0) {
        fprintf(stderr, "'a' には正の数を入れてください。\n");
        fprintf(stderr, " a = "); fgets(s, 128, stdin); sscanf(s, "%lf", &a);
    }

    x = newton(a, eps);
    printf("\n# sqrt(%e) = %.15e\n", a, x);
    return 0;
}