/* Copyright (C) 2014 by Alexandru Cojocaru */
/* This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see . */
/* Compile: gcc -std=c11 -lgmp */
#include
#include
#include
void
salamin_brent (int k, mpf_t pi)
{
mpf_t a, a_k_1, b, c, s, t;
mpf_inits (a, a_k_1, b, c, s, t, NULL);
mpf_set_ui (a_k_1, 1);
mpf_sqrt_ui (b, 2);
mpf_ui_div (b, 1, b);
mpf_set_ui (s, 2);
mpf_ui_div (s, 1, s);
for (int i = 1; i <= k; ++i) {
mpf_add (a, a_k_1, b);
mpf_div_ui (a, a, 2);
mpf_mul (b, a_k_1, b);
mpf_sqrt (b, b);
mpf_mul (c, a, a);
mpf_mul (t, b, b);
mpf_sub (c, c, t);
mpf_set_ui (t, 2);
mpf_pow_ui (t, t, i);
mpf_mul (t, t, c);
mpf_sub (s, s, t);
mpf_set (a_k_1, a);
}
mpf_set (pi, a);
mpf_pow_ui (pi, a, 2);
mpf_mul_ui (pi, pi, 2);
mpf_div (pi, pi, s);
}
int
main (void)
{
mpf_set_default_prec (140);
mpf_t pi;
mpf_init (pi);
salamin_brent (5, pi);
mpf_sub_ui (pi, pi, 3);
if (mpf_out_str (stdout, 10, 30, pi) == 0)
error (1, 0, "mpf_out_str failed");
}