/*  Newton's method for finding sqrt 2 using the mpfr library
    Written 2016 by Eric Olson for Math/CS 466/666 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <mpfr.h>

// N is bits of precision
#define N 100

mpfr_t rf,rdf,rg;
mpfr_ptr f(mpfr_t x){
    mpfr_mul(rf,x,x,MPFR_RNDN);
    mpfr_sub_ui(rf,rf,2,MPFR_RNDN);
    return rf;
}
mpfr_ptr df(mpfr_t x){
    mpfr_mul_ui(rdf,x,2,MPFR_RNDN);
    return rdf;
}
mpfr_ptr g(mpfr_t x){
    mpfr_ptr y=f(x), dy=df(x);
    mpfr_div(rg,y,dy,MPFR_RNDN);
    mpfr_sub(rg,x,rg,MPFR_RNDN);
    return rg;
}

int main(){
    mpfr_init2(rf,N); mpfr_init2(rdf,N);
    mpfr_init2(rg,N);

    mpfr_t xn,en,sqrt2;
    mpfr_init2(xn,N); mpfr_init2(en,N);
    mpfr_init2(sqrt2,N);
    mpfr_set_ui(xn,1,MPFR_RNDN);
    mpfr_sqrt_ui(sqrt2,2,MPFR_RNDN);

    int digits=N*log(2.0)/log(10.0);
    printf("Using %d digit arithmetic...\n",digits);
    for(int i=0;i<20;i++){
        mpfr_ptr xnp1=g(xn);
        if(mpfr_cmp(xnp1,xn)==0) {
            printf("...converged because xn=xnp1.\n");
            break;
        }
        mpfr_set(xn,xnp1,MPFR_RNDN);
        mpfr_sub(en,xn,sqrt2,MPFR_RNDN);
        mpfr_printf("en=%.4Re\n",en);
    }
    mpfr_printf("xn=%.*Rf\n",digits-1,xn);
    return 0;
}