Skip to content

Exact float conv #13

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Exact float conv #13

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
45e077f
Add missing overflow checks to rational_in_float()
Nov 16, 2019
ce6ae72
Fix range check to exclude NaN's
Nov 19, 2019
c5818d7
Exact float to rationial conversion
Nov 19, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 43 additions & 3 deletions expected/pg_rational_test.out
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,9 +51,9 @@ select 0.263157894737::float::rational;
(1 row)

select 3.141592625359::float::rational;
rational
-----------------
4712235/1499951
rational
--------------------
109795040/34948847
(1 row)

select 0.606557377049::float::rational;
Expand All @@ -68,6 +68,46 @@ select -0.5::float::rational;
-1/2
(1 row)

select 1.000001::float::rational;
rational
-----------------
1000001/1000000
(1 row)

select 1.0000001::float::rational;
rational
-------------------
10000001/10000000
(1 row)

select 1.00000001::float::rational;
rational
---------------------
100000001/100000000
(1 row)

select 1.000000001::float::rational;
rational
---------------------
999999918/999999917
(1 row)

select 1.0000000001::float::rational;
rational
----------
1/1
(1 row)

select 2147483647::float::rational;
rational
--------------
2147483647/1
(1 row)

select 2147483647.1::float::rational;
ERROR: value too large for rational
select 'NAN'::float::rational;
ERROR: value too large for rational
-- to float
select '1/2'::rational::float;
float8
Expand Down
5 changes: 5 additions & 0 deletions pg_rational--0.0.1.sql
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,6 +131,11 @@ RETURNS rational
AS '$libdir/pg_rational'
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;

CREATE FUNCTION rational_limit_denominator(rational, integer)
RETURNS rational
AS '$libdir/pg_rational'
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;

CREATE FUNCTION rational_intermediate(rational, rational)
RETURNS rational
AS '$libdir/pg_rational'
Expand Down
149 changes: 115 additions & 34 deletions pg_rational.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include "libpq/pqformat.h" /* needed for send/recv functions */
#include <limits.h>
#include <math.h>
#include <float.h>

PG_MODULE_MAGIC;

Expand All @@ -13,6 +14,7 @@ typedef struct
int32 denom;
} Rational;

static void limit_denominator(Rational *, int64, int64, int32);
static int32 gcd(int32, int32);
static bool simplify(Rational *);
static int32 cmp(Rational *, Rational *);
Expand Down Expand Up @@ -102,46 +104,36 @@ rational_in(PG_FUNCTION_ARGS)
PG_RETURN_POINTER(result);
}

/*
This function taken from John Kennedy's paper, "Algorithm To Convert a
Decimal to a Fraction." Translated from Pascal.
*/

Datum
rational_in_float(PG_FUNCTION_ARGS)
{
float8 target = PG_GETARG_FLOAT8(0),
z,
prev_denom,
error;
int32 temp,
sign;
float_part;
int exponent,
off;
int64 d, n;
Rational *result = palloc(sizeof(Rational));
const int32 max_denominator = INT32_MAX;
const int32 max_numerator = INT32_MAX;

if (target == floor(target))
{
result->numer = floor(target);
result->denom = 1;
PG_RETURN_POINTER(result);
}

sign = target < 0.0 ? -1 : 1;
target = fabs(target);

z = target;
prev_denom = 0;
result->denom = 1;
do
{
z = 1.0 / (z - floor(z));
temp = result->denom;
result->denom = result->denom * floor(z) + prev_denom;
prev_denom = temp;
result->numer = round(target * result->denom);

error = fabs(target - ((float8) result->numer / (float8) result->denom));
} while (z != floor(z) && error >= 1e-12);

result->numer *= sign;
if (!(fabs(target) <= max_numerator)) // also excludes NaN's
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("value too large for rational")));

// convert target into a fraction n/d (with d being a power of
// 2). It is exact as long as target isn't too small, then it
// looses precion because it's rounded below 2^-63.

float_part = frexp(target, &exponent);
exponent = DBL_MANT_DIG - exponent;
off = 0;
if (exponent >= 63)
off = exponent - 62;
n = round(ldexp(float_part, DBL_MANT_DIG-off));
d = (int64)1 << (exponent-off);
limit_denominator(result, n, d, max_denominator);
PG_RETURN_POINTER(result);
}

Expand Down Expand Up @@ -226,13 +218,27 @@ rational_send(PG_FUNCTION_ARGS)
************* ARITHMETIC **************
*/

PG_FUNCTION_INFO_V1(rational_limit_denominator);
PG_FUNCTION_INFO_V1(rational_simplify);
PG_FUNCTION_INFO_V1(rational_add);
PG_FUNCTION_INFO_V1(rational_sub);
PG_FUNCTION_INFO_V1(rational_mul);
PG_FUNCTION_INFO_V1(rational_div);
PG_FUNCTION_INFO_V1(rational_neg);

Datum
rational_limit_denominator(PG_FUNCTION_ARGS)
{
Rational *in = (Rational *) PG_GETARG_POINTER(0);
int32 limit = PG_GETARG_INT32(1);
Rational *out = palloc(sizeof(Rational));

limit_denominator(out, in->numer, in->denom, limit);

PG_RETURN_POINTER(out);
}


Datum
rational_simplify(PG_FUNCTION_ARGS)
{
Expand Down Expand Up @@ -464,6 +470,81 @@ rational_larger(PG_FUNCTION_ARGS)
************** INTERNAL ***************
*/

/*
limit_denomintor() uses continued fractions to convert the
rational n/d into the rational n'/d' with d' < max_denominator
and n' <= INT32_MAX and smallest |d/n-d'/n'|.
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a paper or a description of this algorithm you can link to in the comment?

*/
void
limit_denominator(Rational * r, int64 n, int64 d, int32 max_denominator)
{
float8 target,
error1,
error2,
df;
int neg, k, kn;
int64 a, d1;
int64 p0, q0, p1, q1, p2, q2;
const int32 max_numerator = INT32_MAX;

target = (float8)n / (float8)d;
neg = false;
if (n < 0)
{
neg = true;
n = -n;
}
p0 = 0;
q0 = 1;
p1 = 1;
q1 = 0;
while (true)
{
a = n / d;
q2 = q0 + a * q1;
if (q2 > max_denominator)
break;
p2 = p0 + a * p1;
if (p2 > max_numerator)
break;
d1 = n - a * d;
n = d;
d = d1;
p0 = p1;
q0 = q1;
p1 = p2;
q1 = q2;
if (d == 0 || target == (float8)p1 / (float8)q1)
break;
}
// calculate secondary convergent (reuse variables p2, q2)
// take largest possible k.
k = (max_denominator - q0) / q1;
if (p1 != 0) {
kn = (max_numerator - p0) / p1;
if (kn < k)
k = kn;
}
p2 = p0 + k * p1;
q2 = q0 + k * q1;
// select best of both solutions
error1 = fabs((float8)p1 / (float8)q1 - target);
error2 = fabs((float8)p2 / (float8)q2 - target);
df = error2 - error1;
if (df < 0 || (df == 0.0 && q2 < q1))
{
r->numer = p2;
r->denom = q2;
}
else
{
r->numer = p1;
r->denom = q1;
}
if (neg)
r->numer = -r->numer;
}

int32
gcd(int32 a, int32 b)
{
Expand Down
8 changes: 8 additions & 0 deletions sql/pg_rational_test.sql
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,14 @@ select 0.263157894737::float::rational;
select 3.141592625359::float::rational;
select 0.606557377049::float::rational;
select -0.5::float::rational;
select 1.000001::float::rational;
select 1.0000001::float::rational;
select 1.00000001::float::rational;
select 1.000000001::float::rational;
select 1.0000000001::float::rational;
select 2147483647::float::rational;
select 2147483647.1::float::rational;
select 'NAN'::float::rational;

-- to float
select '1/2'::rational::float;
Expand Down