Issue related to mpf_set_default_prec()

[tsx-11 tsx-11]

Hi,

Reproducible with gcc/SunOS on Intel/SPARC with -O0, -O2, -O3.

Tested with Chudnovsky Algorithm: http://beej.us/blog/data/pi-chudnovsky-gmp/chudnovsky_c.txt

$ gcc -O2 -Wall -o chudnovsky chudnovsky.c -lgmp
$ chudnovsky 100000

Application uses mpf_set_z()/mpf_add()/mpf_mul()/mpf_div() with high precision (much bigger than mpf_get_default_prec()==64).


All mpf_ variables are initialized like:

  mpf_set_default_prec(precision_bits);
  mpf_inits(result, con, A, B, F, sum, NULL);

In this case, chudnovsky 100000 produces correct result.


Now, we change initialization of all mpf_ variables in a way that mpf_set_default_prec() is never called:

  // mpf_set_default_prec(precision_bits);
  // mpf_inits(result, con, A, B, F, sum, NULL);
  mpf_init2(result,precision_bits);
  mpf_init2(con,precision_bits);
  mpf_init2(A,precision_bits);
  mpf_init2(B,precision_bits);
  mpf_init2(F,precision_bits);
  mpf_init2(sum,precision_bits);

In this case, chudnovsky 100000 produces a wrong result (when comparing last digits).


==> it seems that mpf_set_z()/mpf_add()/mpf_mul()/mpf_div() somewhere internally anyway require mpf_set_default_prec() to be called before, and this is not good.

I am not 100% sure but it seems that precision is getting lost here:

//mpf_set_default_prec(precision_bits);
                mpf_set_z(B, c);
//mpf_set_default_prec(64);

Do I miss something here?

Thank you.



chudnovsky.c located on http://beej.us/blog/data/pi-chudnovsky-gmp/chudnovsky_c.txt:

/*
* Compute pi to a certain number of decimal digits, and print it.
*
*   gcc -O2 -Wall -o chudnovsky chudnovsky.c -lgmp
*
* WARNING: This is a demonstration program only, is not optimized for
* speed, and should not be used for serious work!
*
* The Chudnovsky Algorithm:
*                               _____
*                     426880 * /10005
*  pi = ---------------------------------------------
*         _inf_
*         \     (6*k)! * (13591409 + 545140134 * k)
*          \    -----------------------------------
*          /     (3*k)! * (k!)^3 * (-640320)^(3*k)
*         /____
*          k=0
*
* http://en.wikipedia.org/wiki/Pi#Rapidly_convergent_series
*
* First million digits: http://www.piday.org/million.php
*
* Copyright (c) 2012 Brian "Beej Jorgensen" Hall <beej at beej.us>
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
* CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
* TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
* SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>

// how many to display if the user doesn't specify:
#define DEFAULT_DIGITS 60

// how many decimal digits the algorithm generates per iteration:
#define DIGITS_PER_ITERATION 14.1816474627254776555

/**
 * Compute pi to the specified number of decimal digits using the
 * Chudnovsky Algorithm.
 *
 * http://en.wikipedia.org/wiki/Pi#Rapidly_convergent_series
 *
 * NOTE: this function returns a malloc()'d string!
 *
 * @param digits number of decimal digits to compute
 *
 * @return a malloc'd string result (with no decimal marker)
 */
char *chudnovsky(unsigned long digits)
{
    mpf_t result, con, A, B, F, sum;
    mpz_t a, b, c, d, e;
    char *output;
    mp_exp_t exp;
    double bits_per_digit;

    unsigned long int k, threek;
    unsigned long iterations = (digits/DIGITS_PER_ITERATION)+1;
    unsigned long precision_bits;

    // roughly compute how many bits of precision we need for
    // this many digit:
    bits_per_digit = 3.32192809488736234789; // log2(10)
    precision_bits = (digits * bits_per_digit) + 1;

    mpf_set_default_prec(precision_bits);

    // allocate GMP variables
    mpf_inits(result, con, A, B, F, sum, NULL);
    mpz_inits(a, b, c, d, e, NULL);

    mpf_set_ui(sum, 0); // sum already zero at this point, so just FYI

    // first the constant sqrt part
    mpf_sqrt_ui(con, 10005);
    mpf_mul_ui(con, con, 426880);

    // now the fun bit
    for (k = 0; k < iterations; k++) {
        threek = 3*k;

        mpz_fac_ui(a, 6*k);  // (6k)!

        mpz_set_ui(b, 545140134); // 13591409 + 545140134k
        mpz_mul_ui(b, b, k);
        mpz_add_ui(b, b, 13591409);

        mpz_fac_ui(c, threek);  // (3k)!

        mpz_fac_ui(d, k);  // (k!)^3
        mpz_pow_ui(d, d, 3);

        mpz_ui_pow_ui(e, 640320, threek); // -640320^(3k)
        if ((threek&1) == 1) { mpz_neg(e, e); }

        // numerator (in A)
        mpz_mul(a, a, b);
        mpf_set_z(A, a);

        // denominator (in B)
        mpz_mul(c, c, d);
        mpz_mul(c, c, e);
        mpf_set_z(B, c);

        // result
        mpf_div(F, A, B);

        // add on to sum
        mpf_add(sum, sum, F);
    }

    // final calculations (solve for pi)
    mpf_ui_div(sum, 1, sum); // invert result
    mpf_mul(sum, sum, con); // multiply by constant sqrt part

    // get result base-10 in a string:
    output = mpf_get_str(NULL, &exp, 10, digits, sum); // calls malloc()

    // free GMP variables
    mpf_clears(result, con, A, B, F, sum, NULL);
    mpz_clears(a, b, c, d, e, NULL);

    return output;
}

/**
 * Print a usage message and exit
 */
void usage_exit(void)
{
    fprintf(stderr, "usage: chudnovsky [digits]\n");
    exit(1);
}

/**
 * MAIN
 *
 * See usage_exit() for usage.
 */
int main(int argc, char **argv)
{
    char *pi, *endptr;
    long digits;

    switch (argc) {
        case 1:
            digits = DEFAULT_DIGITS;
            break;

        case 2:
            digits = strtol(argv[1], &endptr, 10);
            if (*endptr != '\0') { usage_exit(); }
            break;

        default:
            usage_exit();
    }

    if (digits < 1) { usage_exit(); }

    pi = chudnovsky(digits);

    // since there's no decimal point in the result, we'll print the
    // first digit, then the rest of it, with the expectation that the
    // decimal will appear after "3", as per usual:
    printf("%.1s.%s\n", pi, pi+1);

    // chudnovsky() malloc()s the result string, so let's be proper:
    free(pi);

    return 0;
}