From 45e077f716a669a05587bfcfebb531bba4b7e044 Mon Sep 17 00:00:00 2001 From: Andreas Degert Date: Sat, 16 Nov 2019 10:04:27 +0100 Subject: [PATCH 1/3] Add missing overflow checks to rational_in_float() Overflow leads to wrong results and can create non-terminating loops in the postgresql backend. Changes: - return ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE when input is too large to be represented as rational - end refinement when the resulting numerator or denominator would not fit into an int32. - added some test cases --- expected/pg_rational_test.out | 38 +++++++++++++++++++++++++++++++++++ pg_rational.c | 29 +++++++++++++++++++------- sql/pg_rational_test.sql | 7 +++++++ 3 files changed, 67 insertions(+), 7 deletions(-) diff --git a/expected/pg_rational_test.out b/expected/pg_rational_test.out index 2cc784f..8ebd323 100644 --- a/expected/pg_rational_test.out +++ b/expected/pg_rational_test.out @@ -68,6 +68,44 @@ 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 +------------------ + 10000000/9999999 +(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 -- to float select '1/2'::rational::float; float8 diff --git a/pg_rational.c b/pg_rational.c index 2b21838..f28067a 100644 --- a/pg_rational.c +++ b/pg_rational.c @@ -111,15 +111,17 @@ rational_in_float(PG_FUNCTION_ARGS) { float8 target = PG_GETARG_FLOAT8(0), z, - prev_denom, + fnumer, + fdenom, error; int32 temp, + prev_denom, sign; Rational *result = palloc(sizeof(Rational)); - if (target == floor(target)) + if (target == (int32)target) { - result->numer = floor(target); + result->numer = (int32)target; result->denom = 1; PG_RETURN_POINTER(result); } @@ -127,16 +129,29 @@ rational_in_float(PG_FUNCTION_ARGS) sign = target < 0.0 ? -1 : 1; target = fabs(target); + if (target > INT32_MAX) { + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("value too large for rational"))); + } z = target; prev_denom = 0; + result->numer = (int32)round(target); 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); + fdenom = result->denom * floor(z) + prev_denom; + if (fdenom > INT32_MAX) { + break; + } + fnumer = round(target * fdenom); + if (fnumer > INT32_MAX) { + break; + } + prev_denom = result->denom; + result->numer = (int32)fnumer; + result->denom = (int32)fdenom; error = fabs(target - ((float8) result->numer / (float8) result->denom)); } while (z != floor(z) && error >= 1e-12); diff --git a/sql/pg_rational_test.sql b/sql/pg_rational_test.sql index 4ad74f4..fe3ef3e 100644 --- a/sql/pg_rational_test.sql +++ b/sql/pg_rational_test.sql @@ -23,6 +23,13 @@ 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; -- to float select '1/2'::rational::float; From ce6ae729bd0c2d759a67bfe591cd52d7e39a0e32 Mon Sep 17 00:00:00 2001 From: Andreas Degert Date: Tue, 19 Nov 2019 13:53:35 +0100 Subject: [PATCH 2/3] Fix range check to exclude NaN's - use double negation because conditions with NaN's always evaluate to false - add test case with NaN - some cosmetic changes --- pg_rational.c | 11 +++-------- sql/pg_rational_test.sql | 1 + 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/pg_rational.c b/pg_rational.c index f28067a..dbd3ddd 100644 --- a/pg_rational.c +++ b/pg_rational.c @@ -114,8 +114,7 @@ rational_in_float(PG_FUNCTION_ARGS) fnumer, fdenom, error; - int32 temp, - prev_denom, + int32 prev_denom, sign; Rational *result = palloc(sizeof(Rational)); @@ -129,7 +128,7 @@ rational_in_float(PG_FUNCTION_ARGS) sign = target < 0.0 ? -1 : 1; target = fabs(target); - if (target > INT32_MAX) { + if (!(target <= INT32_MAX)) { // also excludes NaN's ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("value too large for rational"))); @@ -142,13 +141,9 @@ rational_in_float(PG_FUNCTION_ARGS) { z = 1.0 / (z - floor(z)); fdenom = result->denom * floor(z) + prev_denom; - if (fdenom > INT32_MAX) { - break; - } fnumer = round(target * fdenom); - if (fnumer > INT32_MAX) { + if (fnumer > INT32_MAX || fdenom > INT32_MAX ) break; - } prev_denom = result->denom; result->numer = (int32)fnumer; result->denom = (int32)fdenom; diff --git a/sql/pg_rational_test.sql b/sql/pg_rational_test.sql index fe3ef3e..c6359db 100644 --- a/sql/pg_rational_test.sql +++ b/sql/pg_rational_test.sql @@ -30,6 +30,7 @@ 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; From c5818d752cc4b8e7420c400c2c1b94d92200b135 Mon Sep 17 00:00:00 2001 From: Andreas Degert Date: Tue, 19 Nov 2019 14:40:35 +0100 Subject: [PATCH 3/3] Exact float to rationial conversion This version of rational_in_float() converts the float into an int64/int64 rational and calculates the int32/int32 rational with the least error (where possible reproducing it exactly) using continued fractions. Additionally a funktion rational_limit_denominator() is definined which finds the best approximation for a rational with the denominator not exceeding a given value. One use case is to recover fractions that are stored with limited precision: rational_limit_denominator(0.3333::float, 1000) -> 1/3 --- expected/pg_rational_test.out | 14 ++-- pg_rational--0.0.1.sql | 5 ++ pg_rational.c | 151 +++++++++++++++++++++++++--------- 3 files changed, 124 insertions(+), 46 deletions(-) diff --git a/expected/pg_rational_test.out b/expected/pg_rational_test.out index 8ebd323..0aaff0e 100644 --- a/expected/pg_rational_test.out +++ b/expected/pg_rational_test.out @@ -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; @@ -75,9 +75,9 @@ select 1.000001::float::rational; (1 row) select 1.0000001::float::rational; - rational ------------------- - 10000000/9999999 + rational +------------------- + 10000001/10000000 (1 row) select 1.00000001::float::rational; @@ -106,6 +106,8 @@ select 2147483647::float::rational; 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 diff --git a/pg_rational--0.0.1.sql b/pg_rational--0.0.1.sql index 5be3671..5b8ebb8 100644 --- a/pg_rational--0.0.1.sql +++ b/pg_rational--0.0.1.sql @@ -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' diff --git a/pg_rational.c b/pg_rational.c index dbd3ddd..ac540c6 100644 --- a/pg_rational.c +++ b/pg_rational.c @@ -4,6 +4,7 @@ #include "libpq/pqformat.h" /* needed for send/recv functions */ #include #include +#include PG_MODULE_MAGIC; @@ -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 *); @@ -102,56 +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, - fnumer, - fdenom, - error; - int32 prev_denom, - 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 == (int32)target) - { - result->numer = (int32)target; - result->denom = 1; - PG_RETURN_POINTER(result); - } - - sign = target < 0.0 ? -1 : 1; - target = fabs(target); - - if (!(target <= INT32_MAX)) { // also excludes NaN's + if (!(fabs(target) <= max_numerator)) // also excludes NaN's ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("value too large for rational"))); - } - z = target; - prev_denom = 0; - result->numer = (int32)round(target); - result->denom = 1; - do - { - z = 1.0 / (z - floor(z)); - fdenom = result->denom * floor(z) + prev_denom; - fnumer = round(target * fdenom); - if (fnumer > INT32_MAX || fdenom > INT32_MAX ) - break; - prev_denom = result->denom; - result->numer = (int32)fnumer; - result->denom = (int32)fdenom; - error = fabs(target - ((float8) result->numer / (float8) result->denom)); - } while (z != floor(z) && error >= 1e-12); - - result->numer *= sign; + // 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); } @@ -236,6 +218,7 @@ 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); @@ -243,6 +226,19 @@ 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) { @@ -474,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'|. +*/ +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) {