Re: Issue related to mpf_set_default_prec()

[paul zimmermann]

   Hi,

I could not reproduce your problem. Please could you send the exact modified
program you use?

Paul Zimmermann
___
gmp-bugs mailing list
gmp-bugs@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-bugs

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 10

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 10 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 10 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 
*
* 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 
#include 
#include 
#include 

// 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,