MPFR дает NaN ответ, когда все входные данные верны - PullRequest
Новая
       173

MPFR дает NaN ответ, когда все входные данные верны

0 голосов
/ 31 марта 2020

У меня есть небольшая программа, написанная на C, которая использует MFPR и GMP для вычисления числа Пи во многих десятичных разрядах, используя алгоритм Чудновского .

Вот код: 

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

int main() {
  mpfr_t m, m1, m2;
  mpfr_t l;
  mpz_t x, x1;
  mpfr_t series, s1;
  mpfr_t c;

  mpfr_set_default_prec(1000000);
  mpfr_init_set_ui(m, 1, MPFR_RNDN);
  mpfr_init(m1);
  mpfr_init(m2);

  mpfr_init_set_ui(l, 13591409, MPFR_RNDN);

  mpz_init_set_ui(x, 1);
  mpz_init_set_str(x1, "-262537412640768000", 10);

  mpfr_init_set_ui(series, 0, MPFR_RNDN);
  mpfr_init(s1);

  mpfr_init_set_ui(c, 10005, MPFR_RNDN);
  mpfr_sqrt(c, c, MPFR_RNDN);
  mpfr_mul_ui(c, c, 426880, MPFR_RNDN);

  for(unsigned int k = 0; k < 1000; ++k) {
    // series += (m * l) / x
    mpfr_mul(s1, m, l, MPFR_RNDN);
    mpfr_printf("m * l: %Rf\n", s1);
    mpfr_div_z(s1, s1, x, MPFR_RNDN);
    mpfr_printf("(m * l) / x: %Rf\n", s1);

    mpfr_add(series, series, s1, MPFR_RNDN);
    mpfr_printf("series: %Rf\n", series);

    // (12k + 2)(12k + 6)(12k + 10)
    mpfr_set_ui(m1, 12, MPFR_RNDN);
    mpfr_mul_ui(m1, m1, k, MPFR_RNDN);
    mpfr_add_ui(m1, m1, 2, MPFR_RNDN);

    mpfr_set_ui(m2, 12, MPFR_RNDN);
    mpfr_mul_ui(m2, m2, k, MPFR_RNDN);
    mpfr_add_ui(m2, m2, 6, MPFR_RNDN);
    mpfr_mul(m1, m1, m2, MPFR_RNDN);

    mpfr_set_ui(m2, 12, MPFR_RNDN);
    mpfr_mul_ui(m2, m2, k, MPFR_RNDN);
    mpfr_add_ui(m2, m2, 10, MPFR_RNDN);
    mpfr_mul(m1, m1, m2, MPFR_RNDN);

    // (k + 1)^3
    mpfr_set_ui(m2, k + 1, MPFR_RNDN);
    mpfr_pow_ui(m2, m2, 3, MPFR_RNDN);

    // (12k + 2)(12k + 6)(12k + 10)
    // ----------------------------
    //          (k + 1)^3
    mpfr_div(m, m1, m2, MPFR_RNDN);

    // l += 545140134
    mpfr_add_ui(l, l, 545140134, MPFR_RNDN);

    // x *= -262537412640768000
    mpz_mul(x, x, x1);
  }

  mpfr_printf("m: %Rf\nl: %Rf\n", m, l);
  gmp_printf("x: %Zd\n", x);
  mpfr_printf("series: %Rf\n", series);

  mpfr_div(c, c, series, MPFR_RNDN);
  mpfr_printf("pi: %Rf\n", c);

  mpfr_clear(m);
  mpfr_clear(m1);
  mpfr_clear(m2);

  mpfr_clear(l);

  mpz_clear(x);
  mpz_clear(x1);

  mpfr_clear(series);
  mpfr_clear(s1);

  mpfr_clear(c);

  mpfr_free_cache();
  return 0;
}

Когда я запускаю программу, она печатает значения для m, l, x and series, но всегда печатает nan для s1.

Я пытался разместить операторы печати по всему коду, чтобы рисунок где он идет не так, и проблема, кажется, заключается в вызове mpfr_div_z(s1, s1, x, MPFR_RNDN), так как он продолжает чередоваться со значением -0 и 0. Я гуглил на предмет возможных решений, но ничего не нашел

...