problem when converting from double

[Robby Villegas]

To debug this, you can easily make GMP show you its value for the rational
number, or the values for numerator and denominator separately.

To write out the quotient, use:

  mpq_out_str(stdout, 10, val)

To write out the numerator and denominator separately as integers, use:

  mpz_out_str(stdout, 10, num)
  mpz_out_str(stdout, 10, den)

As an example, I did this for 0.3 (see program attached below), and got:

  GMP says the rational is:
  5404319552844595/18014398509481984
  This output has 34 bytes.

  Numerator = 5404319552844595
  Denominator = 18014398509481984

And there's your problem:  the integers for numerator and denominator are too
big to fit into a C int or long if the latter are only 32 bits.  Thus, I also
saw "junk" when I used your mpz_get_si calls.  The GMP documentation on
"Converting Integers" warns you of fitting too-large values into C integers:

  http://www.swox.com/gmp/manual/Converting-Integers.html#Converting%20Integers


As for why the rational representation uses such big integers, realize that C
doubles typically have 53 bits of precision, and the range of values is
typically DBL_MIN ~= 10^-308 and DBL_MAX ~= 10^308.  For instance, the binary
representation of 0.3 is typically:

  0.10011001100110011001100110011001100110011001100110011 x 2^-1

To write this as a ratio of integers, take the mantissa bits and make an
integer out of them; this yields 5404319552844595.  The denominator
is 2^54, which is 18014398509481984.

So, the output from GMP is not surprising.

You could truncate to the most significant 32 bits of the integers, *if* you
test that the number is in a reasonable range (1/INT_MAX < x < INT_MAX at the
outside).  But you won't be able to do this with all doubles, e.g. 1.0e20 is too
big for a quotient of 32-bit integers.


As others have said, you probably want a rationalizer library, from which
you can get nice quotients like 1/3 instead of bigthing/hugething.

Robby Villegas


- - - - - - - - - - - - - - - - - -  snip  - - - - - - - - - - - - - - - - -

/*
Sample usages of this command:

  thisCommand 0.3
  thisCommand 0.25
  thisCommand 3.14159265358979323e-290

*/


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

#include <gmp.h>


int main(int argc, char * argv[]) {
	mpq_t val;
	mpz_t num,
		den;
	const char * digits;

	double g;
	char * statusPtr;
	size_t bytes;
	
	if (argc != 2) {
		fprintf(stderr, "Usage:\n");
		fprintf(stderr, "  <thisCommand> <number>\n");
		exit(1);
	}
	digits = argv[1];

	mpq_init(val);
	mpz_init(num);
	mpz_init(den);
	
	g = strtod(digits, &statusPtr);
	if (statusPtr == digits) {
		fprintf(stderr, "Apparently, the given digits \"%s\" aren't parsible
as a floating point number by strtod().\n", digits);
		exit(1);
	}
	mpq_set_d(val, g);

	printf("GMP says the rational is:\n");
	bytes = mpq_out_str(stdout, 10, val);
	printf("\nThis output has %u bytes.\n", bytes);

	mpq_get_num(num, val);
	mpq_get_den(den, val);
	printf("Numerator = ");
	mpz_out_str(stdout, 10, num);
	printf("\n");
	printf("Denominator = ");
	mpz_out_str(stdout, 10, den);
	printf("\n\n");

	return 0;
}



%localhost (GMP) (49) > ./temp 0.25
GMP says the rational is:
1/4
This output has 3 bytes.
Numerator = 1
Denominator = 4
 


%localhost (GMP) (50) > ./temp 0.3
GMP says the rational is:
5404319552844595/18014398509481984
This output has 34 bytes.
Numerator = 5404319552844595
Denominator = 18014398509481984
 


%localhost (GMP) (51) > ./temp 3.14159265358979323e-290
GMP says the rational is:
5515253462785407/175555970201398037864189960037990696642380564349834626243584063630598316216309534309285622385163609395625111210811907575838661883607828732903171318983861449587663958422720200465138886329341888788528401320395513446131006525725061407689368272012526598792334483090416306874948482361796597953940777665648656384
This output has 323 bytes.
Numerator = 5515253462785407
Denominator = 175555970201398037864189960037990696642380564349834626243584063630598316216309534309285622385163609395625111210811907575838661883607828732903171318983861449587663958422720200465138886329341888788528401320395513446131006525725061407689368272012526598792334483090416306874948482361796597953940777665648656384
%localhost (GMP) (52) >