Wrong division ?

Torbjörn Granlund tg at gmplib.org
Sat Apr 2 09:09:05 UTC 2016


Pedro Gimeno <gmpdiscuss at formauri.es> writes:

  Indeed. I modified the OP's program like this:
  
  #include <stdio.h>
  #include <gmp.h>
  int main(void) {
     mpq_t xx;
     mpf_t result;
     mpq_init(xx);
     mpz_fac_ui(mpq_numref(xx),50);
     mpz_fac_ui(mpq_denref(xx),120);
     mpf_init2(result,50000);
     mpf_set_q(result, xx);
  
     /*-----Print-----*/
     printf("\n\n\n\n\n");
     gmp_printf("%4.1000Ff\n",result);
     mpq_clear(xx);
     mpf_clear(result);
     return 0;
  }
  
Innovative to use mpf_set_q!  I'd have naively set reasonable precisions
for the floats.  Perhaps something like this:

int main(void) {
  mpz_t x,y;
  mpf_t xx,yy,divisione;

  mpz_inits(x, y, NULL);
  mpf_init2(divisione,3322);  // 3322 = ceil(log2(10^1000))

  mpz_fac_ui(x,50);
  mpf_init2(xx, mpz_sizeinbase(x,2));
  mpf_set_z(xx,x);
  mpz_clear(x);

  mpz_fac_ui(y,120);
  mpf_init2(yy, mpz_sizeinbase(y,2));
  mpf_set_z(yy,y);
  mpz_clear(y);

  mpf_div(divisione,xx,yy);

  gmp_printf("%4.1000Ff\n",divisione);

  mpf_clears(xx, yy, divisione, NULL);
  return 0;
}

(Using just mpz_sizeinbase makes sure the values are represented
exactly, but one should make sure the precision is not much greater than
the ultimately required result.  Here it is fine as 120! << 10^1000.

  and the result matches the stated output. Unfortunately, %Qf is not
  yet implemented; I've been wanting it for a long time. The output of
  many exact calculations is a fraction, and having to use an
  intermediate mpf to store the result in order to display it is a waste
  of memory and time.
  
Wouldn't a prospective Qf in essence perform the same division that is
done for mpf_set_q, and thus use the same time and memory?

-- 
Torbjörn
Please encrypt, key id 0xC8601622


More information about the gmp-discuss mailing list