postgresql/src/backend/utils/adt/float.c

3915 lines
89 KiB
C
Raw Normal View History

/*-------------------------------------------------------------------------
*
* float.c
* Functions for the built-in floating-point types.
*
* Portions Copyright (c) 1996-2019, PostgreSQL Global Development Group
* Portions Copyright (c) 1994, Regents of the University of California
*
*
* IDENTIFICATION
2010-09-20 22:08:53 +02:00
* src/backend/utils/adt/float.c
*
*-------------------------------------------------------------------------
*/
#include "postgres.h"
#include <ctype.h>
#include <float.h>
#include <math.h>
#include <limits.h>
#include "catalog/pg_type.h"
#include "common/int.h"
Change floating-point output format for improved performance. Previously, floating-point output was done by rounding to a specific decimal precision; by default, to 6 or 15 decimal digits (losing information) or as requested using extra_float_digits. Drivers that wanted exact float values, and applications like pg_dump that must preserve values exactly, set extra_float_digits=3 (or sometimes 2 for historical reasons, though this isn't enough for float4). Unfortunately, decimal rounded output is slow enough to become a noticable bottleneck when dealing with large result sets or COPY of large tables when many floating-point values are involved. Floating-point output can be done much faster when the output is not rounded to a specific decimal length, but rather is chosen as the shortest decimal representation that is closer to the original float value than to any other value representable in the same precision. The recently published Ryu algorithm by Ulf Adams is both relatively simple and remarkably fast. Accordingly, change float4out/float8out to output shortest decimal representations if extra_float_digits is greater than 0, and make that the new default. Applications that need rounded output can set extra_float_digits back to 0 or below, and take the resulting performance hit. We make one concession to portability for systems with buggy floating-point input: we do not output decimal values that fall exactly halfway between adjacent representable binary values (which would rely on the reader doing round-to-nearest-even correctly). This is known to be a problem at least for VS2013 on Windows. Our version of the Ryu code originates from https://github.com/ulfjack/ryu/ at commit c9c3fb1979, but with the following (significant) modifications: - Output format is changed to use fixed-point notation for small exponents, as printf would, and also to use lowercase 'e', a minimum of 2 exponent digits, and a mandatory sign on the exponent, to keep the formatting as close as possible to previous output. - The output of exact midpoint values is disabled as noted above. - The integer fast-path code is changed somewhat (since we have fixed-point output and the upstream did not). - Our project style has been largely applied to the code with the exception of C99 declaration-after-statement, which has been retained as an exception to our present policy. - Most of upstream's debugging and conditionals are removed, and we use our own configure tests to determine things like uint128 availability. Changing the float output format obviously affects a number of regression tests. This patch uses an explicit setting of extra_float_digits=0 for test output that is not expected to be exactly reproducible (e.g. due to numerical instability or differing algorithms for transcendental functions). Conversions from floats to numeric are unchanged by this patch. These may appear in index expressions and it is not yet clear whether any change should be made, so that can be left for another day. This patch assumes that the only supported floating point format is now IEEE format, and the documentation is updated to reflect that. Code by me, adapting the work of Ulf Adams and other contributors. References: https://dl.acm.org/citation.cfm?id=3192369 Reviewed-by: Tom Lane, Andres Freund, Donald Dong Discussion: https://postgr.es/m/87r2el1bx6.fsf@news-spur.riddles.org.uk
2019-02-13 16:20:33 +01:00
#include "common/shortest_dec.h"
#include "libpq/pqformat.h"
Use a separate random seed for SQL random()/setseed() functions. Previously, the SQL random() function depended on libc's random(3), and setseed() invoked srandom(3). This results in interference between these functions and backend-internal uses of random(3). We'd never paid too much mind to that, but in the wake of commit 88bdbd3f7 which added log_statement_sample_rate, the interference arguably has a security consequence: if log_statement_sample_rate is active then an unprivileged user could probably control which if any of his SQL commands get logged, by issuing setseed() at the right times. That seems bad. To fix this reliably, we need random() and setseed() to use their own private random state variable. Standard random(3) isn't amenable to such usage, so let's switch to pg_erand48(). It's hard to say whether that's more or less "random" than any particular platform's version of random(3), but it does have a wider seed value and a longer period than are required by POSIX, so we can hope that this isn't a big downgrade. Also, we should now have uniform behavior of random() across platforms, which is worth something. While at it, upgrade the per-process seed initialization method to use pg_strong_random() if available, greatly reducing the predictability of the initial seed value. (I'll separately do something similar for the internal uses of random().) In addition to forestalling the possible security problem, this has a benefit in the other direction, which is that we can now document setseed() as guaranteeing a reproducible sequence of random() values. Previously, because of the possibility of internal calls of random(3), we could not promise any such thing. Discussion: https://postgr.es/m/3859.1545849900@sss.pgh.pa.us
2018-12-29 23:33:27 +01:00
#include "miscadmin.h"
#include "utils/array.h"
#include "utils/float.h"
#include "utils/fmgrprotos.h"
#include "utils/sortsupport.h"
Use a separate random seed for SQL random()/setseed() functions. Previously, the SQL random() function depended on libc's random(3), and setseed() invoked srandom(3). This results in interference between these functions and backend-internal uses of random(3). We'd never paid too much mind to that, but in the wake of commit 88bdbd3f7 which added log_statement_sample_rate, the interference arguably has a security consequence: if log_statement_sample_rate is active then an unprivileged user could probably control which if any of his SQL commands get logged, by issuing setseed() at the right times. That seems bad. To fix this reliably, we need random() and setseed() to use their own private random state variable. Standard random(3) isn't amenable to such usage, so let's switch to pg_erand48(). It's hard to say whether that's more or less "random" than any particular platform's version of random(3), but it does have a wider seed value and a longer period than are required by POSIX, so we can hope that this isn't a big downgrade. Also, we should now have uniform behavior of random() across platforms, which is worth something. While at it, upgrade the per-process seed initialization method to use pg_strong_random() if available, greatly reducing the predictability of the initial seed value. (I'll separately do something similar for the internal uses of random().) In addition to forestalling the possible security problem, this has a benefit in the other direction, which is that we can now document setseed() as guaranteeing a reproducible sequence of random() values. Previously, because of the possibility of internal calls of random(3), we could not promise any such thing. Discussion: https://postgr.es/m/3859.1545849900@sss.pgh.pa.us
2018-12-29 23:33:27 +01:00
#include "utils/timestamp.h"
Change floating-point output format for improved performance. Previously, floating-point output was done by rounding to a specific decimal precision; by default, to 6 or 15 decimal digits (losing information) or as requested using extra_float_digits. Drivers that wanted exact float values, and applications like pg_dump that must preserve values exactly, set extra_float_digits=3 (or sometimes 2 for historical reasons, though this isn't enough for float4). Unfortunately, decimal rounded output is slow enough to become a noticable bottleneck when dealing with large result sets or COPY of large tables when many floating-point values are involved. Floating-point output can be done much faster when the output is not rounded to a specific decimal length, but rather is chosen as the shortest decimal representation that is closer to the original float value than to any other value representable in the same precision. The recently published Ryu algorithm by Ulf Adams is both relatively simple and remarkably fast. Accordingly, change float4out/float8out to output shortest decimal representations if extra_float_digits is greater than 0, and make that the new default. Applications that need rounded output can set extra_float_digits back to 0 or below, and take the resulting performance hit. We make one concession to portability for systems with buggy floating-point input: we do not output decimal values that fall exactly halfway between adjacent representable binary values (which would rely on the reader doing round-to-nearest-even correctly). This is known to be a problem at least for VS2013 on Windows. Our version of the Ryu code originates from https://github.com/ulfjack/ryu/ at commit c9c3fb1979, but with the following (significant) modifications: - Output format is changed to use fixed-point notation for small exponents, as printf would, and also to use lowercase 'e', a minimum of 2 exponent digits, and a mandatory sign on the exponent, to keep the formatting as close as possible to previous output. - The output of exact midpoint values is disabled as noted above. - The integer fast-path code is changed somewhat (since we have fixed-point output and the upstream did not). - Our project style has been largely applied to the code with the exception of C99 declaration-after-statement, which has been retained as an exception to our present policy. - Most of upstream's debugging and conditionals are removed, and we use our own configure tests to determine things like uint128 availability. Changing the float output format obviously affects a number of regression tests. This patch uses an explicit setting of extra_float_digits=0 for test output that is not expected to be exactly reproducible (e.g. due to numerical instability or differing algorithms for transcendental functions). Conversions from floats to numeric are unchanged by this patch. These may appear in index expressions and it is not yet clear whether any change should be made, so that can be left for another day. This patch assumes that the only supported floating point format is now IEEE format, and the documentation is updated to reflect that. Code by me, adapting the work of Ulf Adams and other contributors. References: https://dl.acm.org/citation.cfm?id=3192369 Reviewed-by: Tom Lane, Andres Freund, Donald Dong Discussion: https://postgr.es/m/87r2el1bx6.fsf@news-spur.riddles.org.uk
2019-02-13 16:20:33 +01:00
/*
* Configurable GUC parameter
*
* If >0, use shortest-decimal format for output; this is both the default and
* allows for compatibility with clients that explicitly set a value here to
* get round-trip-accurate results. If 0 or less, then use the old, slow,
* decimal rounding method.
*/
int extra_float_digits = 1;
/* Cached constants for degree-based trig functions */
static bool degree_consts_set = false;
static float8 sin_30 = 0;
static float8 one_minus_cos_60 = 0;
static float8 asin_0_5 = 0;
static float8 acos_0_5 = 0;
static float8 atan_1_0 = 0;
static float8 tan_45 = 0;
static float8 cot_45 = 0;
/*
* These are intentionally not static; don't "fix" them. They will never
* be referenced by other files, much less changed; but we don't want the
* compiler to know that, else it might try to precompute expressions
* involving them. See comments for init_degree_constants().
*/
float8 degree_c_thirty = 30.0;
float8 degree_c_forty_five = 45.0;
float8 degree_c_sixty = 60.0;
float8 degree_c_one_half = 0.5;
float8 degree_c_one = 1.0;
Use a separate random seed for SQL random()/setseed() functions. Previously, the SQL random() function depended on libc's random(3), and setseed() invoked srandom(3). This results in interference between these functions and backend-internal uses of random(3). We'd never paid too much mind to that, but in the wake of commit 88bdbd3f7 which added log_statement_sample_rate, the interference arguably has a security consequence: if log_statement_sample_rate is active then an unprivileged user could probably control which if any of his SQL commands get logged, by issuing setseed() at the right times. That seems bad. To fix this reliably, we need random() and setseed() to use their own private random state variable. Standard random(3) isn't amenable to such usage, so let's switch to pg_erand48(). It's hard to say whether that's more or less "random" than any particular platform's version of random(3), but it does have a wider seed value and a longer period than are required by POSIX, so we can hope that this isn't a big downgrade. Also, we should now have uniform behavior of random() across platforms, which is worth something. While at it, upgrade the per-process seed initialization method to use pg_strong_random() if available, greatly reducing the predictability of the initial seed value. (I'll separately do something similar for the internal uses of random().) In addition to forestalling the possible security problem, this has a benefit in the other direction, which is that we can now document setseed() as guaranteeing a reproducible sequence of random() values. Previously, because of the possibility of internal calls of random(3), we could not promise any such thing. Discussion: https://postgr.es/m/3859.1545849900@sss.pgh.pa.us
2018-12-29 23:33:27 +01:00
/* State for drandom() and setseed() */
static bool drandom_seed_set = false;
static unsigned short drandom_seed[3] = {0, 0, 0};
/* Local function prototypes */
static double sind_q1(double x);
static double cosd_q1(double x);
static void init_degree_constants(void);
#ifndef HAVE_CBRT
/*
* Some machines (in particular, some versions of AIX) have an extern
* declaration for cbrt() in <math.h> but fail to provide the actual
* function, which causes configure to not set HAVE_CBRT. Furthermore,
* their compilers spit up at the mismatch between extern declaration
* and static definition. We work around that here by the expedient
* of a #define to make the actual name of the static function different.
*/
#define cbrt my_cbrt
static double cbrt(double x);
Phase 2 of pgindent updates. Change pg_bsd_indent to follow upstream rules for placement of comments to the right of code, and remove pgindent hack that caused comments following #endif to not obey the general rule. Commit e3860ffa4dd0dad0dd9eea4be9cc1412373a8c89 wasn't actually using the published version of pg_bsd_indent, but a hacked-up version that tried to minimize the amount of movement of comments to the right of code. The situation of interest is where such a comment has to be moved to the right of its default placement at column 33 because there's code there. BSD indent has always moved right in units of tab stops in such cases --- but in the previous incarnation, indent was working in 8-space tab stops, while now it knows we use 4-space tabs. So the net result is that in about half the cases, such comments are placed one tab stop left of before. This is better all around: it leaves more room on the line for comment text, and it means that in such cases the comment uniformly starts at the next 4-space tab stop after the code, rather than sometimes one and sometimes two tabs after. Also, ensure that comments following #endif are indented the same as comments following other preprocessor commands such as #else. That inconsistency turns out to have been self-inflicted damage from a poorly-thought-through post-indent "fixup" in pgindent. This patch is much less interesting than the first round of indent changes, but also bulkier, so I thought it best to separate the effects. Discussion: https://postgr.es/m/E1dAmxK-0006EE-1r@gemulon.postgresql.org Discussion: https://postgr.es/m/30527.1495162840@sss.pgh.pa.us
2017-06-21 21:18:54 +02:00
#endif /* HAVE_CBRT */
/*
* Returns -1 if 'val' represents negative infinity, 1 if 'val'
* represents (positive) infinity, and 0 otherwise. On some platforms,
* this is equivalent to the isinf() macro, but not everywhere: C99
* does not specify that isinf() needs to distinguish between positive
* and negative infinity.
*/
int
is_infinite(double val)
{
2004-08-29 07:07:03 +02:00
int inf = isinf(val);
if (inf == 0)
return 0;
else if (val > 0)
return 1;
else
return -1;
}
/* ========== USER I/O ROUTINES ========== */
/*
* float4in - converts "num" to float4
*
* Note that this code now uses strtof(), where it used to use strtod().
*
* The motivation for using strtof() is to avoid a double-rounding problem:
* for certain decimal inputs, if you round the input correctly to a double,
* and then round the double to a float, the result is incorrect in that it
* does not match the result of rounding the decimal value to float directly.
*
* One of the best examples is 7.038531e-26:
*
* 0xAE43FDp-107 = 7.03853069185120912085...e-26
* midpoint 7.03853100000000022281...e-26
* 0xAE43FEp-107 = 7.03853130814879132477...e-26
*
* making 0xAE43FDp-107 the correct float result, but if you do the conversion
* via a double, you get
*
* 0xAE43FD.7FFFFFF8p-107 = 7.03853099999999907487...e-26
* midpoint 7.03853099999999964884...e-26
* 0xAE43FD.80000000p-107 = 7.03853100000000022281...e-26
* 0xAE43FD.80000008p-107 = 7.03853100000000137076...e-26
*
* so the value rounds to the double exactly on the midpoint between the two
* nearest floats, and then rounding again to a float gives the incorrect
* result of 0xAE43FEp-107.
*
*/
Datum
float4in(PG_FUNCTION_ARGS)
{
char *num = PG_GETARG_CSTRING(0);
char *orig_num;
float val;
char *endptr;
/*
2005-10-15 04:49:52 +02:00
* endptr points to the first character _after_ the sequence we recognized
* as a valid floating point number. orig_num points to the original input
* string.
*/
orig_num = num;
/* skip leading whitespace */
while (*num != '\0' && isspace((unsigned char) *num))
num++;
/*
2005-10-15 04:49:52 +02:00
* Check for an empty-string input to begin with, to avoid the vagaries of
* strtod() on different platforms.
*/
if (*num == '\0')
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type %s: \"%s\"",
"real", orig_num)));
errno = 0;
val = strtof(num, &endptr);
/* did we not see anything that looks like a double? */
if (endptr == num || errno != 0)
{
int save_errno = errno;
/*
* C99 requires that strtof() accept NaN, [+-]Infinity, and [+-]Inf,
* but not all platforms support all of these (and some accept them
* but set ERANGE anyway...) Therefore, we check for these inputs
* ourselves if strtof() fails.
*
* Note: C99 also requires hexadecimal input as well as some extended
* forms of NaN, but we consider these forms unportable and don't try
* to support them. You can use 'em if your strtof() takes 'em.
*/
if (pg_strncasecmp(num, "NaN", 3) == 0)
{
val = get_float4_nan();
endptr = num + 3;
}
else if (pg_strncasecmp(num, "Infinity", 8) == 0)
{
val = get_float4_infinity();
endptr = num + 8;
}
else if (pg_strncasecmp(num, "+Infinity", 9) == 0)
{
val = get_float4_infinity();
endptr = num + 9;
}
else if (pg_strncasecmp(num, "-Infinity", 9) == 0)
{
2004-08-29 07:07:03 +02:00
val = -get_float4_infinity();
endptr = num + 9;
}
else if (pg_strncasecmp(num, "inf", 3) == 0)
{
val = get_float4_infinity();
endptr = num + 3;
}
else if (pg_strncasecmp(num, "+inf", 4) == 0)
{
val = get_float4_infinity();
endptr = num + 4;
}
else if (pg_strncasecmp(num, "-inf", 4) == 0)
{
val = -get_float4_infinity();
endptr = num + 4;
}
else if (save_errno == ERANGE)
{
/*
* Some platforms return ERANGE for denormalized numbers (those
* that are not zero, but are too close to zero to have full
* precision). We'd prefer not to throw error for that, so try to
* detect whether it's a "real" out-of-range condition by checking
* to see if the result is zero or huge.
*
* Use isinf() rather than HUGE_VALF on VS2013 because it
* generates a spurious overflow warning for -HUGE_VALF. Also use
* isinf() if HUGE_VALF is missing.
*/
if (val == 0.0 ||
#if !defined(HUGE_VALF) || (defined(_MSC_VER) && (_MSC_VER < 1900))
isinf(val)
#else
(val >= HUGE_VALF || val <= -HUGE_VALF)
#endif
)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("\"%s\" is out of range for type real",
orig_num)));
}
else
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type %s: \"%s\"",
"real", orig_num)));
}
#ifdef HAVE_BUGGY_SOLARIS_STRTOD
else
{
/*
2005-10-15 04:49:52 +02:00
* Many versions of Solaris have a bug wherein strtod sets endptr to
* point one byte beyond the end of the string when given "inf" or
* "infinity".
*/
if (endptr != num && endptr[-1] == '\0')
endptr--;
}
Phase 2 of pgindent updates. Change pg_bsd_indent to follow upstream rules for placement of comments to the right of code, and remove pgindent hack that caused comments following #endif to not obey the general rule. Commit e3860ffa4dd0dad0dd9eea4be9cc1412373a8c89 wasn't actually using the published version of pg_bsd_indent, but a hacked-up version that tried to minimize the amount of movement of comments to the right of code. The situation of interest is where such a comment has to be moved to the right of its default placement at column 33 because there's code there. BSD indent has always moved right in units of tab stops in such cases --- but in the previous incarnation, indent was working in 8-space tab stops, while now it knows we use 4-space tabs. So the net result is that in about half the cases, such comments are placed one tab stop left of before. This is better all around: it leaves more room on the line for comment text, and it means that in such cases the comment uniformly starts at the next 4-space tab stop after the code, rather than sometimes one and sometimes two tabs after. Also, ensure that comments following #endif are indented the same as comments following other preprocessor commands such as #else. That inconsistency turns out to have been self-inflicted damage from a poorly-thought-through post-indent "fixup" in pgindent. This patch is much less interesting than the first round of indent changes, but also bulkier, so I thought it best to separate the effects. Discussion: https://postgr.es/m/E1dAmxK-0006EE-1r@gemulon.postgresql.org Discussion: https://postgr.es/m/30527.1495162840@sss.pgh.pa.us
2017-06-21 21:18:54 +02:00
#endif /* HAVE_BUGGY_SOLARIS_STRTOD */
/* skip trailing whitespace */
while (*endptr != '\0' && isspace((unsigned char) *endptr))
endptr++;
/* if there is any junk left at the end of the string, bail out */
if (*endptr != '\0')
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type %s: \"%s\"",
"real", orig_num)));
PG_RETURN_FLOAT4(val);
}
/*
* float4out - converts a float4 number to a string
* using a standard output format
*/
Datum
float4out(PG_FUNCTION_ARGS)
{
float4 num = PG_GETARG_FLOAT4(0);
Improve snprintf.c's handling of NaN, Infinity, and minus zero. Up to now, float4out/float8out handled NaN and Infinity cases explicitly, and invoked psprintf only for ordinary float values. This was done because platform implementations of snprintf produce varying representations of these special cases. But now that we use snprintf.c always, it's better to give it the responsibility to produce a uniform representation of these cases, so that we have uniformity across the board not only in float4out/float8out. Hence, move that work into fmtfloat(). Also, teach fmtfloat() to recognize IEEE minus zero and handle it correctly. The previous coding worked only accidentally, and would fail for e.g. "%+f" format (it'd print "+-0.00000"). Now that we're using snprintf.c everywhere, it's not acceptable for it to do weird things in corner cases. (This incidentally avoids a portability problem we've seen on some really ancient platforms, that native sprintf does the wrong thing with minus zero.) Also, introduce a new entry point in snprintf.c to allow float[48]out to bypass the work of interpreting a well-known format spec, as well as bypassing the overhead of the psprintf layer. I modeled this API loosely on strfromd(). In my testing, this brings float[48]out back to approximately the same speed they had when using native snprintf, fixing one of the main performance issues caused by using snprintf.c. (There is some talk of more aggressive work to improve the speed of floating-point output conversion, but these changes seem to provide a better starting point for such work anyway.) Getting rid of the previous ad-hoc hack for Infinity/NaN in fmtfloat() allows removing <ctype.h> from snprintf.c's #includes. I also removed a few other #includes that I think are historical, though the buildfarm may expose that as wrong. Discussion: https://postgr.es/m/13178.1538794717@sss.pgh.pa.us
2018-10-08 18:19:20 +02:00
char *ascii = (char *) palloc(32);
int ndig = FLT_DIG + extra_float_digits;
Change floating-point output format for improved performance. Previously, floating-point output was done by rounding to a specific decimal precision; by default, to 6 or 15 decimal digits (losing information) or as requested using extra_float_digits. Drivers that wanted exact float values, and applications like pg_dump that must preserve values exactly, set extra_float_digits=3 (or sometimes 2 for historical reasons, though this isn't enough for float4). Unfortunately, decimal rounded output is slow enough to become a noticable bottleneck when dealing with large result sets or COPY of large tables when many floating-point values are involved. Floating-point output can be done much faster when the output is not rounded to a specific decimal length, but rather is chosen as the shortest decimal representation that is closer to the original float value than to any other value representable in the same precision. The recently published Ryu algorithm by Ulf Adams is both relatively simple and remarkably fast. Accordingly, change float4out/float8out to output shortest decimal representations if extra_float_digits is greater than 0, and make that the new default. Applications that need rounded output can set extra_float_digits back to 0 or below, and take the resulting performance hit. We make one concession to portability for systems with buggy floating-point input: we do not output decimal values that fall exactly halfway between adjacent representable binary values (which would rely on the reader doing round-to-nearest-even correctly). This is known to be a problem at least for VS2013 on Windows. Our version of the Ryu code originates from https://github.com/ulfjack/ryu/ at commit c9c3fb1979, but with the following (significant) modifications: - Output format is changed to use fixed-point notation for small exponents, as printf would, and also to use lowercase 'e', a minimum of 2 exponent digits, and a mandatory sign on the exponent, to keep the formatting as close as possible to previous output. - The output of exact midpoint values is disabled as noted above. - The integer fast-path code is changed somewhat (since we have fixed-point output and the upstream did not). - Our project style has been largely applied to the code with the exception of C99 declaration-after-statement, which has been retained as an exception to our present policy. - Most of upstream's debugging and conditionals are removed, and we use our own configure tests to determine things like uint128 availability. Changing the float output format obviously affects a number of regression tests. This patch uses an explicit setting of extra_float_digits=0 for test output that is not expected to be exactly reproducible (e.g. due to numerical instability or differing algorithms for transcendental functions). Conversions from floats to numeric are unchanged by this patch. These may appear in index expressions and it is not yet clear whether any change should be made, so that can be left for another day. This patch assumes that the only supported floating point format is now IEEE format, and the documentation is updated to reflect that. Code by me, adapting the work of Ulf Adams and other contributors. References: https://dl.acm.org/citation.cfm?id=3192369 Reviewed-by: Tom Lane, Andres Freund, Donald Dong Discussion: https://postgr.es/m/87r2el1bx6.fsf@news-spur.riddles.org.uk
2019-02-13 16:20:33 +01:00
if (extra_float_digits > 0)
{
float_to_shortest_decimal_buf(num, ascii);
PG_RETURN_CSTRING(ascii);
}
Improve snprintf.c's handling of NaN, Infinity, and minus zero. Up to now, float4out/float8out handled NaN and Infinity cases explicitly, and invoked psprintf only for ordinary float values. This was done because platform implementations of snprintf produce varying representations of these special cases. But now that we use snprintf.c always, it's better to give it the responsibility to produce a uniform representation of these cases, so that we have uniformity across the board not only in float4out/float8out. Hence, move that work into fmtfloat(). Also, teach fmtfloat() to recognize IEEE minus zero and handle it correctly. The previous coding worked only accidentally, and would fail for e.g. "%+f" format (it'd print "+-0.00000"). Now that we're using snprintf.c everywhere, it's not acceptable for it to do weird things in corner cases. (This incidentally avoids a portability problem we've seen on some really ancient platforms, that native sprintf does the wrong thing with minus zero.) Also, introduce a new entry point in snprintf.c to allow float[48]out to bypass the work of interpreting a well-known format spec, as well as bypassing the overhead of the psprintf layer. I modeled this API loosely on strfromd(). In my testing, this brings float[48]out back to approximately the same speed they had when using native snprintf, fixing one of the main performance issues caused by using snprintf.c. (There is some talk of more aggressive work to improve the speed of floating-point output conversion, but these changes seem to provide a better starting point for such work anyway.) Getting rid of the previous ad-hoc hack for Infinity/NaN in fmtfloat() allows removing <ctype.h> from snprintf.c's #includes. I also removed a few other #includes that I think are historical, though the buildfarm may expose that as wrong. Discussion: https://postgr.es/m/13178.1538794717@sss.pgh.pa.us
2018-10-08 18:19:20 +02:00
(void) pg_strfromd(ascii, 32, ndig, num);
PG_RETURN_CSTRING(ascii);
}
/*
* float4recv - converts external binary format to float4
*/
Datum
float4recv(PG_FUNCTION_ARGS)
{
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
PG_RETURN_FLOAT4(pq_getmsgfloat4(buf));
}
/*
* float4send - converts float4 to binary format
*/
Datum
float4send(PG_FUNCTION_ARGS)
{
float4 num = PG_GETARG_FLOAT4(0);
StringInfoData buf;
pq_begintypsend(&buf);
pq_sendfloat4(&buf, num);
PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}
/*
* float8in - converts "num" to float8
*/
Datum
float8in(PG_FUNCTION_ARGS)
{
char *num = PG_GETARG_CSTRING(0);
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
PG_RETURN_FLOAT8(float8in_internal(num, NULL, "double precision", num));
}
/* Convenience macro: set *have_error flag (if provided) or throw error */
#define RETURN_ERROR(throw_error, have_error) \
do { \
if (have_error) { \
*have_error = true; \
return 0.0; \
} else { \
throw_error; \
} \
} while (0)
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
/*
* float8in_internal_opt_error - guts of float8in()
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
*
* This is exposed for use by functions that want a reasonably
* platform-independent way of inputting doubles. The behavior is
* essentially like strtod + ereport on error, but note the following
* differences:
* 1. Both leading and trailing whitespace are skipped.
* 2. If endptr_p is NULL, we throw error if there's trailing junk.
* Otherwise, it's up to the caller to complain about trailing junk.
* 3. In event of a syntax error, the report mentions the given type_name
* and prints orig_string as the input; this is meant to support use of
* this function with types such as "box" and "point", where what we are
* parsing here is just a substring of orig_string.
*
* "num" could validly be declared "const char *", but that results in an
* unreasonable amount of extra casting both here and in callers, so we don't.
*
* When "*have_error" flag is provided, it's set instead of throwing an
* error. This is helpful when caller need to handle errors by itself.
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
*/
double
float8in_internal_opt_error(char *num, char **endptr_p,
const char *type_name, const char *orig_string,
bool *have_error)
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
{
double val;
char *endptr;
if (have_error)
*have_error = false;
/* skip leading whitespace */
while (*num != '\0' && isspace((unsigned char) *num))
num++;
/*
2005-10-15 04:49:52 +02:00
* Check for an empty-string input to begin with, to avoid the vagaries of
* strtod() on different platforms.
*/
if (*num == '\0')
RETURN_ERROR(ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type %s: \"%s\"",
type_name, orig_string))),
have_error);
errno = 0;
val = strtod(num, &endptr);
/* did we not see anything that looks like a double? */
if (endptr == num || errno != 0)
{
int save_errno = errno;
/*
* C99 requires that strtod() accept NaN, [+-]Infinity, and [+-]Inf,
* but not all platforms support all of these (and some accept them
* but set ERANGE anyway...) Therefore, we check for these inputs
* ourselves if strtod() fails.
*
* Note: C99 also requires hexadecimal input as well as some extended
* forms of NaN, but we consider these forms unportable and don't try
* to support them. You can use 'em if your strtod() takes 'em.
*/
if (pg_strncasecmp(num, "NaN", 3) == 0)
{
val = get_float8_nan();
endptr = num + 3;
}
else if (pg_strncasecmp(num, "Infinity", 8) == 0)
{
val = get_float8_infinity();
endptr = num + 8;
}
else if (pg_strncasecmp(num, "+Infinity", 9) == 0)
{
val = get_float8_infinity();
endptr = num + 9;
}
else if (pg_strncasecmp(num, "-Infinity", 9) == 0)
{
2004-08-29 07:07:03 +02:00
val = -get_float8_infinity();
endptr = num + 9;
}
else if (pg_strncasecmp(num, "inf", 3) == 0)
{
val = get_float8_infinity();
endptr = num + 3;
}
else if (pg_strncasecmp(num, "+inf", 4) == 0)
{
val = get_float8_infinity();
endptr = num + 4;
}
else if (pg_strncasecmp(num, "-inf", 4) == 0)
{
val = -get_float8_infinity();
endptr = num + 4;
}
else if (save_errno == ERANGE)
{
/*
* Some platforms return ERANGE for denormalized numbers (those
* that are not zero, but are too close to zero to have full
* precision). We'd prefer not to throw error for that, so try to
* detect whether it's a "real" out-of-range condition by checking
* to see if the result is zero or huge.
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
*
* On error, we intentionally complain about double precision not
* the given type name, and we print only the part of the string
* that is the current number.
*/
if (val == 0.0 || val >= HUGE_VAL || val <= -HUGE_VAL)
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
{
char *errnumber = pstrdup(num);
errnumber[endptr - num] = '\0';
RETURN_ERROR(ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("\"%s\" is out of range for type double precision",
errnumber))),
have_error);
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
}
}
else
RETURN_ERROR(ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type "
"%s: \"%s\"",
type_name, orig_string))),
have_error);
}
#ifdef HAVE_BUGGY_SOLARIS_STRTOD
else
{
/*
2005-10-15 04:49:52 +02:00
* Many versions of Solaris have a bug wherein strtod sets endptr to
* point one byte beyond the end of the string when given "inf" or
* "infinity".
*/
if (endptr != num && endptr[-1] == '\0')
endptr--;
}
Phase 2 of pgindent updates. Change pg_bsd_indent to follow upstream rules for placement of comments to the right of code, and remove pgindent hack that caused comments following #endif to not obey the general rule. Commit e3860ffa4dd0dad0dd9eea4be9cc1412373a8c89 wasn't actually using the published version of pg_bsd_indent, but a hacked-up version that tried to minimize the amount of movement of comments to the right of code. The situation of interest is where such a comment has to be moved to the right of its default placement at column 33 because there's code there. BSD indent has always moved right in units of tab stops in such cases --- but in the previous incarnation, indent was working in 8-space tab stops, while now it knows we use 4-space tabs. So the net result is that in about half the cases, such comments are placed one tab stop left of before. This is better all around: it leaves more room on the line for comment text, and it means that in such cases the comment uniformly starts at the next 4-space tab stop after the code, rather than sometimes one and sometimes two tabs after. Also, ensure that comments following #endif are indented the same as comments following other preprocessor commands such as #else. That inconsistency turns out to have been self-inflicted damage from a poorly-thought-through post-indent "fixup" in pgindent. This patch is much less interesting than the first round of indent changes, but also bulkier, so I thought it best to separate the effects. Discussion: https://postgr.es/m/E1dAmxK-0006EE-1r@gemulon.postgresql.org Discussion: https://postgr.es/m/30527.1495162840@sss.pgh.pa.us
2017-06-21 21:18:54 +02:00
#endif /* HAVE_BUGGY_SOLARIS_STRTOD */
/* skip trailing whitespace */
while (*endptr != '\0' && isspace((unsigned char) *endptr))
endptr++;
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
/* report stopping point if wanted, else complain if not end of string */
if (endptr_p)
*endptr_p = endptr;
else if (*endptr != '\0')
RETURN_ERROR(ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type "
"%s: \"%s\"",
type_name, orig_string))),
have_error);
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
return val;
}
/*
* Interface to float8in_internal_opt_error() without "have_error" argument.
*/
double
float8in_internal(char *num, char **endptr_p,
const char *type_name, const char *orig_string)
{
return float8in_internal_opt_error(num, endptr_p, type_name,
orig_string, NULL);
}
/*
* float8out - converts float8 number to a string
* using a standard output format
*/
Datum
float8out(PG_FUNCTION_ARGS)
{
float8 num = PG_GETARG_FLOAT8(0);
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
PG_RETURN_CSTRING(float8out_internal(num));
}
/*
* float8out_internal - guts of float8out()
*
* This is exposed for use by functions that want a reasonably
* platform-independent way of outputting doubles.
* The result is always palloc'd.
*/
char *
float8out_internal(double num)
{
Improve snprintf.c's handling of NaN, Infinity, and minus zero. Up to now, float4out/float8out handled NaN and Infinity cases explicitly, and invoked psprintf only for ordinary float values. This was done because platform implementations of snprintf produce varying representations of these special cases. But now that we use snprintf.c always, it's better to give it the responsibility to produce a uniform representation of these cases, so that we have uniformity across the board not only in float4out/float8out. Hence, move that work into fmtfloat(). Also, teach fmtfloat() to recognize IEEE minus zero and handle it correctly. The previous coding worked only accidentally, and would fail for e.g. "%+f" format (it'd print "+-0.00000"). Now that we're using snprintf.c everywhere, it's not acceptable for it to do weird things in corner cases. (This incidentally avoids a portability problem we've seen on some really ancient platforms, that native sprintf does the wrong thing with minus zero.) Also, introduce a new entry point in snprintf.c to allow float[48]out to bypass the work of interpreting a well-known format spec, as well as bypassing the overhead of the psprintf layer. I modeled this API loosely on strfromd(). In my testing, this brings float[48]out back to approximately the same speed they had when using native snprintf, fixing one of the main performance issues caused by using snprintf.c. (There is some talk of more aggressive work to improve the speed of floating-point output conversion, but these changes seem to provide a better starting point for such work anyway.) Getting rid of the previous ad-hoc hack for Infinity/NaN in fmtfloat() allows removing <ctype.h> from snprintf.c's #includes. I also removed a few other #includes that I think are historical, though the buildfarm may expose that as wrong. Discussion: https://postgr.es/m/13178.1538794717@sss.pgh.pa.us
2018-10-08 18:19:20 +02:00
char *ascii = (char *) palloc(32);
int ndig = DBL_DIG + extra_float_digits;
Change floating-point output format for improved performance. Previously, floating-point output was done by rounding to a specific decimal precision; by default, to 6 or 15 decimal digits (losing information) or as requested using extra_float_digits. Drivers that wanted exact float values, and applications like pg_dump that must preserve values exactly, set extra_float_digits=3 (or sometimes 2 for historical reasons, though this isn't enough for float4). Unfortunately, decimal rounded output is slow enough to become a noticable bottleneck when dealing with large result sets or COPY of large tables when many floating-point values are involved. Floating-point output can be done much faster when the output is not rounded to a specific decimal length, but rather is chosen as the shortest decimal representation that is closer to the original float value than to any other value representable in the same precision. The recently published Ryu algorithm by Ulf Adams is both relatively simple and remarkably fast. Accordingly, change float4out/float8out to output shortest decimal representations if extra_float_digits is greater than 0, and make that the new default. Applications that need rounded output can set extra_float_digits back to 0 or below, and take the resulting performance hit. We make one concession to portability for systems with buggy floating-point input: we do not output decimal values that fall exactly halfway between adjacent representable binary values (which would rely on the reader doing round-to-nearest-even correctly). This is known to be a problem at least for VS2013 on Windows. Our version of the Ryu code originates from https://github.com/ulfjack/ryu/ at commit c9c3fb1979, but with the following (significant) modifications: - Output format is changed to use fixed-point notation for small exponents, as printf would, and also to use lowercase 'e', a minimum of 2 exponent digits, and a mandatory sign on the exponent, to keep the formatting as close as possible to previous output. - The output of exact midpoint values is disabled as noted above. - The integer fast-path code is changed somewhat (since we have fixed-point output and the upstream did not). - Our project style has been largely applied to the code with the exception of C99 declaration-after-statement, which has been retained as an exception to our present policy. - Most of upstream's debugging and conditionals are removed, and we use our own configure tests to determine things like uint128 availability. Changing the float output format obviously affects a number of regression tests. This patch uses an explicit setting of extra_float_digits=0 for test output that is not expected to be exactly reproducible (e.g. due to numerical instability or differing algorithms for transcendental functions). Conversions from floats to numeric are unchanged by this patch. These may appear in index expressions and it is not yet clear whether any change should be made, so that can be left for another day. This patch assumes that the only supported floating point format is now IEEE format, and the documentation is updated to reflect that. Code by me, adapting the work of Ulf Adams and other contributors. References: https://dl.acm.org/citation.cfm?id=3192369 Reviewed-by: Tom Lane, Andres Freund, Donald Dong Discussion: https://postgr.es/m/87r2el1bx6.fsf@news-spur.riddles.org.uk
2019-02-13 16:20:33 +01:00
if (extra_float_digits > 0)
{
double_to_shortest_decimal_buf(num, ascii);
return ascii;
}
Improve snprintf.c's handling of NaN, Infinity, and minus zero. Up to now, float4out/float8out handled NaN and Infinity cases explicitly, and invoked psprintf only for ordinary float values. This was done because platform implementations of snprintf produce varying representations of these special cases. But now that we use snprintf.c always, it's better to give it the responsibility to produce a uniform representation of these cases, so that we have uniformity across the board not only in float4out/float8out. Hence, move that work into fmtfloat(). Also, teach fmtfloat() to recognize IEEE minus zero and handle it correctly. The previous coding worked only accidentally, and would fail for e.g. "%+f" format (it'd print "+-0.00000"). Now that we're using snprintf.c everywhere, it's not acceptable for it to do weird things in corner cases. (This incidentally avoids a portability problem we've seen on some really ancient platforms, that native sprintf does the wrong thing with minus zero.) Also, introduce a new entry point in snprintf.c to allow float[48]out to bypass the work of interpreting a well-known format spec, as well as bypassing the overhead of the psprintf layer. I modeled this API loosely on strfromd(). In my testing, this brings float[48]out back to approximately the same speed they had when using native snprintf, fixing one of the main performance issues caused by using snprintf.c. (There is some talk of more aggressive work to improve the speed of floating-point output conversion, but these changes seem to provide a better starting point for such work anyway.) Getting rid of the previous ad-hoc hack for Infinity/NaN in fmtfloat() allows removing <ctype.h> from snprintf.c's #includes. I also removed a few other #includes that I think are historical, though the buildfarm may expose that as wrong. Discussion: https://postgr.es/m/13178.1538794717@sss.pgh.pa.us
2018-10-08 18:19:20 +02:00
(void) pg_strfromd(ascii, 32, ndig, num);
Improve portability of I/O behavior for the geometric types. Formerly, the geometric I/O routines such as box_in and point_out relied directly on strtod() and sprintf() for conversion of the float8 component values of their data types. However, the behavior of those functions is pretty platform-dependent, especially for edge-case values such as infinities and NaNs. This was exposed by commit acdf2a8b372aec1d, which added test cases involving boxes with infinity endpoints, and immediately failed on Windows and AIX buildfarm members. We solved these problems years ago in the main float8in and float8out functions, so let's fix it by making the geometric types use that code instead of depending directly on the platform-supplied functions. To do this, refactor the float8in code so that it can be used to parse just part of a string, and as a convenience make the guts of float8out usable without going through DirectFunctionCall. While at it, get rid of geo_ops.c's fairly shaky assumptions about the maximum output string length for a double, by having it build results in StringInfo buffers instead of fixed-length strings. In passing, convert all the "invalid input syntax for type foo" messages in this area of the code into "invalid input syntax for type %s" to reduce the number of distinct translatable strings, per recent discussion. We would have needed a fair number of the latter anyway for code-sharing reasons, so we might as well just go whole hog. Note: this patch is by no means intended to guarantee that the geometric types uniformly behave sanely for infinity or NaN component values. But any bugs we have in that line were there all along, they were just harder to reach in a platform-independent way.
2016-03-30 23:25:03 +02:00
return ascii;
}
/*
* float8recv - converts external binary format to float8
*/
Datum
float8recv(PG_FUNCTION_ARGS)
{
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
PG_RETURN_FLOAT8(pq_getmsgfloat8(buf));
}
/*
* float8send - converts float8 to binary format
*/
Datum
float8send(PG_FUNCTION_ARGS)
{
float8 num = PG_GETARG_FLOAT8(0);
StringInfoData buf;
pq_begintypsend(&buf);
pq_sendfloat8(&buf, num);
PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}
/* ========== PUBLIC ROUTINES ========== */
/*
* ======================
* FLOAT4 BASE OPERATIONS
* ======================
*/
/*
* float4abs - returns |arg1| (absolute value)
*/
Datum
float4abs(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
PG_RETURN_FLOAT4((float4) fabs(arg1));
}
/*
* float4um - returns -arg1 (unary minus)
*/
Datum
float4um(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 result;
result = -arg1;
PG_RETURN_FLOAT4(result);
}
Datum
float4up(PG_FUNCTION_ARGS)
{
float4 arg = PG_GETARG_FLOAT4(0);
PG_RETURN_FLOAT4(arg);
}
Datum
float4larger(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
float4 result;
if (float4_gt(arg1, arg2))
result = arg1;
else
result = arg2;
PG_RETURN_FLOAT4(result);
}
Datum
float4smaller(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
float4 result;
if (float4_lt(arg1, arg2))
result = arg1;
else
result = arg2;
PG_RETURN_FLOAT4(result);
}
/*
* ======================
* FLOAT8 BASE OPERATIONS
* ======================
*/
/*
* float8abs - returns |arg1| (absolute value)
*/
Datum
float8abs(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
PG_RETURN_FLOAT8(fabs(arg1));
}
/*
* float8um - returns -arg1 (unary minus)
*/
Datum
float8um(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
result = -arg1;
PG_RETURN_FLOAT8(result);
}
Datum
float8up(PG_FUNCTION_ARGS)
{
float8 arg = PG_GETARG_FLOAT8(0);
PG_RETURN_FLOAT8(arg);
}
Datum
float8larger(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
float8 result;
if (float8_gt(arg1, arg2))
result = arg1;
else
result = arg2;
PG_RETURN_FLOAT8(result);
}
Datum
float8smaller(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
float8 result;
if (float8_lt(arg1, arg2))
result = arg1;
else
result = arg2;
PG_RETURN_FLOAT8(result);
}
/*
* ====================
* ARITHMETIC OPERATORS
* ====================
*/
/*
* float4pl - returns arg1 + arg2
* float4mi - returns arg1 - arg2
* float4mul - returns arg1 * arg2
* float4div - returns arg1 / arg2
*/
Datum
float4pl(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
2007-11-15 22:14:46 +01:00
PG_RETURN_FLOAT4(float4_pl(arg1, arg2));
}
Datum
float4mi(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_FLOAT4(float4_mi(arg1, arg2));
}
Datum
float4mul(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_FLOAT4(float4_mul(arg1, arg2));
}
Datum
float4div(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_FLOAT4(float4_div(arg1, arg2));
}
/*
* float8pl - returns arg1 + arg2
* float8mi - returns arg1 - arg2
* float8mul - returns arg1 * arg2
* float8div - returns arg1 / arg2
*/
Datum
float8pl(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_FLOAT8(float8_pl(arg1, arg2));
}
Datum
float8mi(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_FLOAT8(float8_mi(arg1, arg2));
}
Datum
float8mul(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_FLOAT8(float8_mul(arg1, arg2));
}
Datum
float8div(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_FLOAT8(float8_div(arg1, arg2));
}
/*
* ====================
* COMPARISON OPERATORS
* ====================
*/
/*
* float4{eq,ne,lt,le,gt,ge} - float4/float4 comparison operations
*/
Fix GiST index build for NaN values in geometric types. GiST index build could go into an infinite loop when presented with boxes (or points, circles or polygons) containing NaN component values. This happened essentially because the code assumed that x == x is true for any "double" value x; but it's not true for NaNs. The looping behavior was not the only problem though: we also attempted to sort the items using simple double comparisons. Since NaNs violate the trichotomy law, qsort could (in principle at least) get arbitrarily confused and mess up the sorting of ordinary values as well as NaNs. And we based splitting choices on box size calculations that could produce NaNs, again resulting in undesirable behavior. To fix, replace all comparisons of doubles in this logic with float8_cmp_internal, which is NaN-aware and is careful to sort NaNs consistently, higher than any non-NaN. Also rearrange the box size calculation to not produce NaNs; instead it should produce an infinity for a box with NaN on one side and not-NaN on the other. I don't by any means claim that this solves all problems with NaNs in geometric values, but it should at least make GiST index insertion work reliably with such data. It's likely that the index search side of things still needs some work, and probably regular geometric operations too. But with this patch we're laying down a convention for how such cases ought to behave. Per bug #14238 from Guang-Dih Lei. Back-patch to 9.2; the code used before commit 7f3bd86843e5aad8 is quite different and doesn't lock up on my simple test case, nor on the submitter's dataset. Report: <20160708151747.1426.60150@wrigleys.postgresql.org> Discussion: <28685.1468246504@sss.pgh.pa.us>
2016-07-15 00:45:59 +02:00
int
float4_cmp_internal(float4 a, float4 b)
{
if (float4_gt(a, b))
return 1;
if (float4_lt(a, b))
return -1;
return 0;
}
Datum
float4eq(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float4_eq(arg1, arg2));
}
Datum
float4ne(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float4_ne(arg1, arg2));
}
Datum
float4lt(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float4_lt(arg1, arg2));
}
Datum
float4le(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float4_le(arg1, arg2));
}
Datum
float4gt(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float4_gt(arg1, arg2));
}
Datum
float4ge(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float4_ge(arg1, arg2));
}
Datum
btfloat4cmp(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_INT32(float4_cmp_internal(arg1, arg2));
}
static int
btfloat4fastcmp(Datum x, Datum y, SortSupport ssup)
{
float4 arg1 = DatumGetFloat4(x);
float4 arg2 = DatumGetFloat4(y);
return float4_cmp_internal(arg1, arg2);
}
Datum
btfloat4sortsupport(PG_FUNCTION_ARGS)
{
SortSupport ssup = (SortSupport) PG_GETARG_POINTER(0);
ssup->comparator = btfloat4fastcmp;
PG_RETURN_VOID();
}
/*
* float8{eq,ne,lt,le,gt,ge} - float8/float8 comparison operations
*/
Fix GiST index build for NaN values in geometric types. GiST index build could go into an infinite loop when presented with boxes (or points, circles or polygons) containing NaN component values. This happened essentially because the code assumed that x == x is true for any "double" value x; but it's not true for NaNs. The looping behavior was not the only problem though: we also attempted to sort the items using simple double comparisons. Since NaNs violate the trichotomy law, qsort could (in principle at least) get arbitrarily confused and mess up the sorting of ordinary values as well as NaNs. And we based splitting choices on box size calculations that could produce NaNs, again resulting in undesirable behavior. To fix, replace all comparisons of doubles in this logic with float8_cmp_internal, which is NaN-aware and is careful to sort NaNs consistently, higher than any non-NaN. Also rearrange the box size calculation to not produce NaNs; instead it should produce an infinity for a box with NaN on one side and not-NaN on the other. I don't by any means claim that this solves all problems with NaNs in geometric values, but it should at least make GiST index insertion work reliably with such data. It's likely that the index search side of things still needs some work, and probably regular geometric operations too. But with this patch we're laying down a convention for how such cases ought to behave. Per bug #14238 from Guang-Dih Lei. Back-patch to 9.2; the code used before commit 7f3bd86843e5aad8 is quite different and doesn't lock up on my simple test case, nor on the submitter's dataset. Report: <20160708151747.1426.60150@wrigleys.postgresql.org> Discussion: <28685.1468246504@sss.pgh.pa.us>
2016-07-15 00:45:59 +02:00
int
float8_cmp_internal(float8 a, float8 b)
{
if (float8_gt(a, b))
return 1;
if (float8_lt(a, b))
return -1;
return 0;
}
Datum
float8eq(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_eq(arg1, arg2));
}
Datum
float8ne(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_ne(arg1, arg2));
}
Datum
float8lt(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_lt(arg1, arg2));
}
Datum
float8le(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_le(arg1, arg2));
}
Datum
float8gt(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_gt(arg1, arg2));
}
Datum
float8ge(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_ge(arg1, arg2));
}
Datum
btfloat8cmp(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
}
static int
btfloat8fastcmp(Datum x, Datum y, SortSupport ssup)
{
float8 arg1 = DatumGetFloat8(x);
float8 arg2 = DatumGetFloat8(y);
return float8_cmp_internal(arg1, arg2);
}
Datum
btfloat8sortsupport(PG_FUNCTION_ARGS)
{
SortSupport ssup = (SortSupport) PG_GETARG_POINTER(0);
ssup->comparator = btfloat8fastcmp;
PG_RETURN_VOID();
}
Datum
btfloat48cmp(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
/* widen float4 to float8 and then compare */
PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
}
Datum
btfloat84cmp(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
/* widen float4 to float8 and then compare */
PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
}
/*
* in_range support function for float8.
*
* Note: we needn't supply a float8_float4 variant, as implicit coercion
* of the offset value takes care of that scenario just as well.
*/
Datum
in_range_float8_float8(PG_FUNCTION_ARGS)
{
float8 val = PG_GETARG_FLOAT8(0);
float8 base = PG_GETARG_FLOAT8(1);
float8 offset = PG_GETARG_FLOAT8(2);
bool sub = PG_GETARG_BOOL(3);
bool less = PG_GETARG_BOOL(4);
float8 sum;
/*
* Reject negative or NaN offset. Negative is per spec, and NaN is
* because appropriate semantics for that seem non-obvious.
*/
if (isnan(offset) || offset < 0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PRECEDING_OR_FOLLOWING_SIZE),
errmsg("invalid preceding or following size in window function")));
/*
* Deal with cases where val and/or base is NaN, following the rule that
* NaN sorts after non-NaN (cf float8_cmp_internal). The offset cannot
* affect the conclusion.
*/
if (isnan(val))
{
if (isnan(base))
PG_RETURN_BOOL(true); /* NAN = NAN */
else
PG_RETURN_BOOL(!less); /* NAN > non-NAN */
}
else if (isnan(base))
{
PG_RETURN_BOOL(less); /* non-NAN < NAN */
}
/*
* Deal with infinite offset (necessarily +inf, at this point). We must
* special-case this because if base happens to be -inf, their sum would
* be NaN, which is an overflow-ish condition we should avoid.
*/
if (isinf(offset))
{
PG_RETURN_BOOL(sub ? !less : less);
}
/*
* Otherwise it should be safe to compute base +/- offset. We trust the
* FPU to cope if base is +/-inf or the true sum would overflow, and
* produce a suitably signed infinity, which will compare properly against
* val whether or not that's infinity.
*/
if (sub)
sum = base - offset;
else
sum = base + offset;
if (less)
PG_RETURN_BOOL(val <= sum);
else
PG_RETURN_BOOL(val >= sum);
}
/*
* in_range support function for float4.
*
* We would need a float4_float8 variant in any case, so we supply that and
* let implicit coercion take care of the float4_float4 case.
*/
Datum
in_range_float4_float8(PG_FUNCTION_ARGS)
{
float4 val = PG_GETARG_FLOAT4(0);
float4 base = PG_GETARG_FLOAT4(1);
float8 offset = PG_GETARG_FLOAT8(2);
bool sub = PG_GETARG_BOOL(3);
bool less = PG_GETARG_BOOL(4);
float8 sum;
/*
* Reject negative or NaN offset. Negative is per spec, and NaN is
* because appropriate semantics for that seem non-obvious.
*/
if (isnan(offset) || offset < 0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PRECEDING_OR_FOLLOWING_SIZE),
errmsg("invalid preceding or following size in window function")));
/*
* Deal with cases where val and/or base is NaN, following the rule that
* NaN sorts after non-NaN (cf float8_cmp_internal). The offset cannot
* affect the conclusion.
*/
if (isnan(val))
{
if (isnan(base))
PG_RETURN_BOOL(true); /* NAN = NAN */
else
PG_RETURN_BOOL(!less); /* NAN > non-NAN */
}
else if (isnan(base))
{
PG_RETURN_BOOL(less); /* non-NAN < NAN */
}
/*
* Deal with infinite offset (necessarily +inf, at this point). We must
* special-case this because if base happens to be -inf, their sum would
* be NaN, which is an overflow-ish condition we should avoid.
*/
if (isinf(offset))
{
PG_RETURN_BOOL(sub ? !less : less);
}
/*
* Otherwise it should be safe to compute base +/- offset. We trust the
* FPU to cope if base is +/-inf or the true sum would overflow, and
* produce a suitably signed infinity, which will compare properly against
* val whether or not that's infinity.
*/
if (sub)
sum = base - offset;
else
sum = base + offset;
if (less)
PG_RETURN_BOOL(val <= sum);
else
PG_RETURN_BOOL(val >= sum);
}
/*
* ===================
* CONVERSION ROUTINES
* ===================
*/
/*
* ftod - converts a float4 number to a float8 number
*/
Datum
ftod(PG_FUNCTION_ARGS)
{
float4 num = PG_GETARG_FLOAT4(0);
PG_RETURN_FLOAT8((float8) num);
}
/*
* dtof - converts a float8 number to a float4 number
*/
Datum
dtof(PG_FUNCTION_ARGS)
{
float8 num = PG_GETARG_FLOAT8(0);
check_float4_val((float4) num, isinf(num), num == 0);
PG_RETURN_FLOAT4((float4) num);
}
/*
* dtoi4 - converts a float8 number to an int4 number
*/
Datum
dtoi4(PG_FUNCTION_ARGS)
{
float8 num = PG_GETARG_FLOAT8(0);
/*
* Get rid of any fractional part in the input. This is so we don't fail
* on just-out-of-range values that would round into range. Note
* assumption that rint() will pass through a NaN or Inf unchanged.
*/
num = rint(num);
/* Range check */
if (unlikely(isnan(num) || !FLOAT8_FITS_IN_INT32(num)))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("integer out of range")));
PG_RETURN_INT32((int32) num);
}
/*
* dtoi2 - converts a float8 number to an int2 number
*/
Datum
dtoi2(PG_FUNCTION_ARGS)
{
float8 num = PG_GETARG_FLOAT8(0);
/*
* Get rid of any fractional part in the input. This is so we don't fail
* on just-out-of-range values that would round into range. Note
* assumption that rint() will pass through a NaN or Inf unchanged.
*/
num = rint(num);
/* Range check */
if (unlikely(isnan(num) || !FLOAT8_FITS_IN_INT16(num)))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("smallint out of range")));
PG_RETURN_INT16((int16) num);
}
/*
* i4tod - converts an int4 number to a float8 number
*/
Datum
i4tod(PG_FUNCTION_ARGS)
{
int32 num = PG_GETARG_INT32(0);
PG_RETURN_FLOAT8((float8) num);
}
/*
* i2tod - converts an int2 number to a float8 number
*/
Datum
i2tod(PG_FUNCTION_ARGS)
{
int16 num = PG_GETARG_INT16(0);
PG_RETURN_FLOAT8((float8) num);
}
/*
* ftoi4 - converts a float4 number to an int4 number
*/
Datum
ftoi4(PG_FUNCTION_ARGS)
{
float4 num = PG_GETARG_FLOAT4(0);
/*
* Get rid of any fractional part in the input. This is so we don't fail
* on just-out-of-range values that would round into range. Note
* assumption that rint() will pass through a NaN or Inf unchanged.
*/
num = rint(num);
/* Range check */
if (unlikely(isnan(num) || !FLOAT4_FITS_IN_INT32(num)))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("integer out of range")));
PG_RETURN_INT32((int32) num);
}
/*
* ftoi2 - converts a float4 number to an int2 number
*/
Datum
ftoi2(PG_FUNCTION_ARGS)
{
float4 num = PG_GETARG_FLOAT4(0);
/*
* Get rid of any fractional part in the input. This is so we don't fail
* on just-out-of-range values that would round into range. Note
* assumption that rint() will pass through a NaN or Inf unchanged.
*/
num = rint(num);
/* Range check */
if (unlikely(isnan(num) || !FLOAT4_FITS_IN_INT16(num)))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("smallint out of range")));
PG_RETURN_INT16((int16) num);
}
/*
* i4tof - converts an int4 number to a float4 number
*/
Datum
i4tof(PG_FUNCTION_ARGS)
{
int32 num = PG_GETARG_INT32(0);
PG_RETURN_FLOAT4((float4) num);
}
/*
* i2tof - converts an int2 number to a float4 number
*/
Datum
i2tof(PG_FUNCTION_ARGS)
{
int16 num = PG_GETARG_INT16(0);
PG_RETURN_FLOAT4((float4) num);
}
/*
* =======================
* RANDOM FLOAT8 OPERATORS
* =======================
*/
/*
* dround - returns ROUND(arg1)
*/
Datum
dround(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
PG_RETURN_FLOAT8(rint(arg1));
}
/*
* dceil - returns the smallest integer greater than or
* equal to the specified float
*/
Datum
dceil(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
PG_RETURN_FLOAT8(ceil(arg1));
}
/*
* dfloor - returns the largest integer lesser than or
* equal to the specified float
*/
Datum
dfloor(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
PG_RETURN_FLOAT8(floor(arg1));
}
/*
* dsign - returns -1 if the argument is less than 0, 0
* if the argument is equal to 0, and 1 if the
* argument is greater than zero.
*/
Datum
dsign(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
if (arg1 > 0)
result = 1.0;
else if (arg1 < 0)
result = -1.0;
else
result = 0.0;
PG_RETURN_FLOAT8(result);
}
/*
* dtrunc - returns truncation-towards-zero of arg1,
* arg1 >= 0 ... the greatest integer less
* than or equal to arg1
* arg1 < 0 ... the least integer greater
* than or equal to arg1
*/
Datum
dtrunc(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
if (arg1 >= 0)
result = floor(arg1);
else
result = -floor(-arg1);
PG_RETURN_FLOAT8(result);
}
/*
* dsqrt - returns square root of arg1
*/
Datum
dsqrt(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
if (arg1 < 0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
errmsg("cannot take square root of a negative number")));
result = sqrt(arg1);
check_float8_val(result, isinf(arg1), arg1 == 0);
PG_RETURN_FLOAT8(result);
}
/*
* dcbrt - returns cube root of arg1
*/
Datum
dcbrt(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
result = cbrt(arg1);
check_float8_val(result, isinf(arg1), arg1 == 0);
PG_RETURN_FLOAT8(result);
}
/*
* dpow - returns pow(arg1,arg2)
*/
Datum
dpow(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
float8 result;
/*
* The POSIX spec says that NaN ^ 0 = 1, and 1 ^ NaN = 1, while all other
* cases with NaN inputs yield NaN (with no error). Many older platforms
* get one or more of these cases wrong, so deal with them via explicit
* logic rather than trusting pow(3).
*/
if (isnan(arg1))
{
if (isnan(arg2) || arg2 != 0.0)
PG_RETURN_FLOAT8(get_float8_nan());
PG_RETURN_FLOAT8(1.0);
}
if (isnan(arg2))
{
if (arg1 != 1.0)
PG_RETURN_FLOAT8(get_float8_nan());
PG_RETURN_FLOAT8(1.0);
}
/*
2005-10-15 04:49:52 +02:00
* The SQL spec requires that we emit a particular SQLSTATE error code for
* certain error conditions. Specifically, we don't return a
* divide-by-zero error code for 0 ^ -1.
*/
if (arg1 == 0 && arg2 < 0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
errmsg("zero raised to a negative power is undefined")));
if (arg1 < 0 && floor(arg2) != arg2)
ereport(ERROR,
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
errmsg("a negative number raised to a non-integer power yields a complex result")));
/*
* pow() sets errno only on some platforms, depending on whether it
2007-11-15 22:14:46 +01:00
* follows _IEEE_, _POSIX_, _XOPEN_, or _SVID_, so we try to avoid using
* errno. However, some platform/CPU combinations return errno == EDOM
* and result == NaN for negative arg1 and very large arg2 (they must be
2007-11-15 22:14:46 +01:00
* using something different from our floor() test to decide it's
* invalid). Other platforms (HPPA) return errno == ERANGE and a large
* (HUGE_VAL) but finite result to signal overflow.
*/
errno = 0;
result = pow(arg1, arg2);
if (errno == EDOM && isnan(result))
{
if ((fabs(arg1) > 1 && arg2 >= 0) || (fabs(arg1) < 1 && arg2 < 0))
/* The sign of Inf is not significant in this case. */
result = get_float8_infinity();
else if (fabs(arg1) != 1)
result = 0;
else
result = 1;
}
else if (errno == ERANGE && result != 0 && !isinf(result))
result = get_float8_infinity();
2007-11-15 22:14:46 +01:00
check_float8_val(result, isinf(arg1) || isinf(arg2), arg1 == 0);
PG_RETURN_FLOAT8(result);
}
/*
* dexp - returns the exponential function of arg1
*/
Datum
dexp(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
errno = 0;
result = exp(arg1);
if (errno == ERANGE && result != 0 && !isinf(result))
result = get_float8_infinity();
check_float8_val(result, isinf(arg1), false);
PG_RETURN_FLOAT8(result);
}
/*
* dlog1 - returns the natural logarithm of arg1
*/
Datum
dlog1(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/*
2005-10-15 04:49:52 +02:00
* Emit particular SQLSTATE error codes for ln(). This is required by the
* SQL standard.
*/
if (arg1 == 0.0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
errmsg("cannot take logarithm of zero")));
if (arg1 < 0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
errmsg("cannot take logarithm of a negative number")));
result = log(arg1);
check_float8_val(result, isinf(arg1), arg1 == 1);
PG_RETURN_FLOAT8(result);
}
/*
* dlog10 - returns the base 10 logarithm of arg1
*/
Datum
dlog10(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/*
2005-10-15 04:49:52 +02:00
* Emit particular SQLSTATE error codes for log(). The SQL spec doesn't
* define log(), but it does define ln(), so it makes sense to emit the
* same error code for an analogous error condition.
*/
if (arg1 == 0.0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
errmsg("cannot take logarithm of zero")));
if (arg1 < 0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
errmsg("cannot take logarithm of a negative number")));
result = log10(arg1);
check_float8_val(result, isinf(arg1), arg1 == 1);
PG_RETURN_FLOAT8(result);
}
/*
* dacos - returns the arccos of arg1 (radians)
*/
Datum
dacos(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/* Per the POSIX spec, return NaN if the input is NaN */
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
/*
* The principal branch of the inverse cosine function maps values in the
* range [-1, 1] to values in the range [0, Pi], so we should reject any
* inputs outside that range and the result will always be finite.
*/
if (arg1 < -1.0 || arg1 > 1.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
result = acos(arg1);
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* dasin - returns the arcsin of arg1 (radians)
*/
Datum
dasin(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/* Per the POSIX spec, return NaN if the input is NaN */
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
/*
* The principal branch of the inverse sine function maps values in the
* range [-1, 1] to values in the range [-Pi/2, Pi/2], so we should reject
* any inputs outside that range and the result will always be finite.
*/
if (arg1 < -1.0 || arg1 > 1.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
result = asin(arg1);
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* datan - returns the arctan of arg1 (radians)
*/
Datum
datan(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/* Per the POSIX spec, return NaN if the input is NaN */
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
/*
* The principal branch of the inverse tangent function maps all inputs to
* values in the range [-Pi/2, Pi/2], so the result should always be
* finite, even if the input is infinite.
*/
result = atan(arg1);
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* atan2 - returns the arctan of arg1/arg2 (radians)
*/
Datum
datan2(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
float8 result;
/* Per the POSIX spec, return NaN if either input is NaN */
if (isnan(arg1) || isnan(arg2))
PG_RETURN_FLOAT8(get_float8_nan());
/*
* atan2 maps all inputs to values in the range [-Pi, Pi], so the result
* should always be finite, even if the inputs are infinite.
*/
result = atan2(arg1, arg2);
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* dcos - returns the cosine of arg1 (radians)
*/
Datum
dcos(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/* Per the POSIX spec, return NaN if the input is NaN */
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
/*
* cos() is periodic and so theoretically can work for all finite inputs,
* but some implementations may choose to throw error if the input is so
* large that there are no significant digits in the result. So we should
* check for errors. POSIX allows an error to be reported either via
* errno or via fetestexcept(), but currently we only support checking
* errno. (fetestexcept() is rumored to report underflow unreasonably
* early on some platforms, so it's not clear that believing it would be a
* net improvement anyway.)
*
* For infinite inputs, POSIX specifies that the trigonometric functions
* should return a domain error; but we won't notice that unless the
* platform reports via errno, so also explicitly test for infinite
* inputs.
*/
errno = 0;
result = cos(arg1);
if (errno != 0 || isinf(arg1))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* dcot - returns the cotangent of arg1 (radians)
*/
Datum
dcot(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/* Per the POSIX spec, return NaN if the input is NaN */
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
/* Be sure to throw an error if the input is infinite --- see dcos() */
errno = 0;
result = tan(arg1);
if (errno != 0 || isinf(arg1))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
result = 1.0 / result;
check_float8_val(result, true /* cot(0) == Inf */ , true);
PG_RETURN_FLOAT8(result);
}
/*
* dsin - returns the sine of arg1 (radians)
*/
Datum
dsin(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/* Per the POSIX spec, return NaN if the input is NaN */
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
/* Be sure to throw an error if the input is infinite --- see dcos() */
errno = 0;
result = sin(arg1);
if (errno != 0 || isinf(arg1))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* dtan - returns the tangent of arg1 (radians)
*/
Datum
dtan(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/* Per the POSIX spec, return NaN if the input is NaN */
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
/* Be sure to throw an error if the input is infinite --- see dcos() */
errno = 0;
result = tan(arg1);
if (errno != 0 || isinf(arg1))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
check_float8_val(result, true /* tan(pi/2) == Inf */ , true);
PG_RETURN_FLOAT8(result);
}
/* ========== DEGREE-BASED TRIGONOMETRIC FUNCTIONS ========== */
/*
* Initialize the cached constants declared at the head of this file
* (sin_30 etc). The fact that we need those at all, let alone need this
* Rube-Goldberg-worthy method of initializing them, is because there are
* compilers out there that will precompute expressions such as sin(constant)
* using a sin() function different from what will be used at runtime. If we
* want exact results, we must ensure that none of the scaling constants used
* in the degree-based trig functions are computed that way. To do so, we
* compute them from the variables degree_c_thirty etc, which are also really
* constants, but the compiler cannot assume that.
*
* Other hazards we are trying to forestall with this kluge include the
* possibility that compilers will rearrange the expressions, or compute
* some intermediate results in registers wider than a standard double.
*
* In the places where we use these constants, the typical pattern is like
* volatile float8 sin_x = sin(x * RADIANS_PER_DEGREE);
* return (sin_x / sin_30);
* where we hope to get a value of exactly 1.0 from the division when x = 30.
* The volatile temporary variable is needed on machines with wide float
* registers, to ensure that the result of sin(x) is rounded to double width
* the same as the value of sin_30 has been. Experimentation with gcc shows
* that marking the temp variable volatile is necessary to make the store and
* reload actually happen; hopefully the same trick works for other compilers.
* (gcc's documentation suggests using the -ffloat-store compiler switch to
* ensure this, but that is compiler-specific and it also pessimizes code in
* many places where we don't care about this.)
*/
static void
init_degree_constants(void)
{
sin_30 = sin(degree_c_thirty * RADIANS_PER_DEGREE);
one_minus_cos_60 = 1.0 - cos(degree_c_sixty * RADIANS_PER_DEGREE);
asin_0_5 = asin(degree_c_one_half);
acos_0_5 = acos(degree_c_one_half);
atan_1_0 = atan(degree_c_one);
tan_45 = sind_q1(degree_c_forty_five) / cosd_q1(degree_c_forty_five);
cot_45 = cosd_q1(degree_c_forty_five) / sind_q1(degree_c_forty_five);
degree_consts_set = true;
}
#define INIT_DEGREE_CONSTANTS() \
do { \
if (!degree_consts_set) \
init_degree_constants(); \
} while(0)
/*
* asind_q1 - returns the inverse sine of x in degrees, for x in
* the range [0, 1]. The result is an angle in the
* first quadrant --- [0, 90] degrees.
*
* For the 3 special case inputs (0, 0.5 and 1), this
* function will return exact values (0, 30 and 90
* degrees respectively).
*/
static double
asind_q1(double x)
{
/*
* Stitch together inverse sine and cosine functions for the ranges [0,
* 0.5] and (0.5, 1]. Each expression below is guaranteed to return
* exactly 30 for x=0.5, so the result is a continuous monotonic function
* over the full range.
*/
if (x <= 0.5)
{
volatile float8 asin_x = asin(x);
return (asin_x / asin_0_5) * 30.0;
}
else
{
volatile float8 acos_x = acos(x);
return 90.0 - (acos_x / acos_0_5) * 60.0;
}
}
/*
* acosd_q1 - returns the inverse cosine of x in degrees, for x in
* the range [0, 1]. The result is an angle in the
* first quadrant --- [0, 90] degrees.
*
* For the 3 special case inputs (0, 0.5 and 1), this
* function will return exact values (0, 60 and 90
* degrees respectively).
*/
static double
acosd_q1(double x)
{
/*
* Stitch together inverse sine and cosine functions for the ranges [0,
* 0.5] and (0.5, 1]. Each expression below is guaranteed to return
* exactly 60 for x=0.5, so the result is a continuous monotonic function
* over the full range.
*/
if (x <= 0.5)
{
volatile float8 asin_x = asin(x);
return 90.0 - (asin_x / asin_0_5) * 30.0;
}
else
{
volatile float8 acos_x = acos(x);
return (acos_x / acos_0_5) * 60.0;
}
}
/*
* dacosd - returns the arccos of arg1 (degrees)
*/
Datum
dacosd(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/* Per the POSIX spec, return NaN if the input is NaN */
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
INIT_DEGREE_CONSTANTS();
/*
* The principal branch of the inverse cosine function maps values in the
* range [-1, 1] to values in the range [0, 180], so we should reject any
* inputs outside that range and the result will always be finite.
*/
if (arg1 < -1.0 || arg1 > 1.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
if (arg1 >= 0.0)
result = acosd_q1(arg1);
else
result = 90.0 + asind_q1(-arg1);
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* dasind - returns the arcsin of arg1 (degrees)
*/
Datum
dasind(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/* Per the POSIX spec, return NaN if the input is NaN */
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
INIT_DEGREE_CONSTANTS();
/*
* The principal branch of the inverse sine function maps values in the
* range [-1, 1] to values in the range [-90, 90], so we should reject any
* inputs outside that range and the result will always be finite.
*/
if (arg1 < -1.0 || arg1 > 1.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
if (arg1 >= 0.0)
result = asind_q1(arg1);
else
result = -asind_q1(-arg1);
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* datand - returns the arctan of arg1 (degrees)
*/
Datum
datand(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
volatile float8 atan_arg1;
/* Per the POSIX spec, return NaN if the input is NaN */
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
INIT_DEGREE_CONSTANTS();
/*
* The principal branch of the inverse tangent function maps all inputs to
* values in the range [-90, 90], so the result should always be finite,
* even if the input is infinite. Additionally, we take care to ensure
* than when arg1 is 1, the result is exactly 45.
*/
atan_arg1 = atan(arg1);
result = (atan_arg1 / atan_1_0) * 45.0;
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* atan2d - returns the arctan of arg1/arg2 (degrees)
*/
Datum
datan2d(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
float8 result;
volatile float8 atan2_arg1_arg2;
/* Per the POSIX spec, return NaN if either input is NaN */
if (isnan(arg1) || isnan(arg2))
PG_RETURN_FLOAT8(get_float8_nan());
INIT_DEGREE_CONSTANTS();
/*
* atan2d maps all inputs to values in the range [-180, 180], so the
* result should always be finite, even if the inputs are infinite.
*
* Note: this coding assumes that atan(1.0) is a suitable scaling constant
* to get an exact result from atan2(). This might well fail on us at
* some point, requiring us to decide exactly what inputs we think we're
* going to guarantee an exact result for.
*/
atan2_arg1_arg2 = atan2(arg1, arg2);
result = (atan2_arg1_arg2 / atan_1_0) * 45.0;
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* sind_0_to_30 - returns the sine of an angle that lies between 0 and
* 30 degrees. This will return exactly 0 when x is 0,
* and exactly 0.5 when x is 30 degrees.
*/
static double
sind_0_to_30(double x)
{
volatile float8 sin_x = sin(x * RADIANS_PER_DEGREE);
return (sin_x / sin_30) / 2.0;
}
/*
* cosd_0_to_60 - returns the cosine of an angle that lies between 0
* and 60 degrees. This will return exactly 1 when x
* is 0, and exactly 0.5 when x is 60 degrees.
*/
static double
cosd_0_to_60(double x)
{
volatile float8 one_minus_cos_x = 1.0 - cos(x * RADIANS_PER_DEGREE);
return 1.0 - (one_minus_cos_x / one_minus_cos_60) / 2.0;
}
/*
* sind_q1 - returns the sine of an angle in the first quadrant
* (0 to 90 degrees).
*/
static double
sind_q1(double x)
{
/*
* Stitch together the sine and cosine functions for the ranges [0, 30]
* and (30, 90]. These guarantee to return exact answers at their
* endpoints, so the overall result is a continuous monotonic function
* that gives exact results when x = 0, 30 and 90 degrees.
*/
if (x <= 30.0)
return sind_0_to_30(x);
else
return cosd_0_to_60(90.0 - x);
}
/*
* cosd_q1 - returns the cosine of an angle in the first quadrant
* (0 to 90 degrees).
*/
static double
cosd_q1(double x)
{
/*
* Stitch together the sine and cosine functions for the ranges [0, 60]
* and (60, 90]. These guarantee to return exact answers at their
* endpoints, so the overall result is a continuous monotonic function
* that gives exact results when x = 0, 60 and 90 degrees.
*/
if (x <= 60.0)
return cosd_0_to_60(x);
else
return sind_0_to_30(90.0 - x);
}
/*
* dcosd - returns the cosine of arg1 (degrees)
*/
Datum
dcosd(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
int sign = 1;
/*
* Per the POSIX spec, return NaN if the input is NaN and throw an error
* if the input is infinite.
*/
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
if (isinf(arg1))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
INIT_DEGREE_CONSTANTS();
/* Reduce the range of the input to [0,90] degrees */
arg1 = fmod(arg1, 360.0);
if (arg1 < 0.0)
{
/* cosd(-x) = cosd(x) */
arg1 = -arg1;
}
if (arg1 > 180.0)
{
/* cosd(360-x) = cosd(x) */
arg1 = 360.0 - arg1;
}
if (arg1 > 90.0)
{
/* cosd(180-x) = -cosd(x) */
arg1 = 180.0 - arg1;
sign = -sign;
}
result = sign * cosd_q1(arg1);
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* dcotd - returns the cotangent of arg1 (degrees)
*/
Datum
dcotd(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
volatile float8 cot_arg1;
int sign = 1;
/*
* Per the POSIX spec, return NaN if the input is NaN and throw an error
* if the input is infinite.
*/
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
if (isinf(arg1))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
INIT_DEGREE_CONSTANTS();
/* Reduce the range of the input to [0,90] degrees */
arg1 = fmod(arg1, 360.0);
if (arg1 < 0.0)
{
/* cotd(-x) = -cotd(x) */
arg1 = -arg1;
sign = -sign;
}
if (arg1 > 180.0)
{
/* cotd(360-x) = -cotd(x) */
arg1 = 360.0 - arg1;
sign = -sign;
}
if (arg1 > 90.0)
{
/* cotd(180-x) = -cotd(x) */
arg1 = 180.0 - arg1;
sign = -sign;
}
cot_arg1 = cosd_q1(arg1) / sind_q1(arg1);
result = sign * (cot_arg1 / cot_45);
/*
* On some machines we get cotd(270) = minus zero, but this isn't always
* true. For portability, and because the user constituency for this
* function probably doesn't want minus zero, force it to plain zero.
*/
if (result == 0.0)
result = 0.0;
check_float8_val(result, true /* cotd(0) == Inf */ , true);
PG_RETURN_FLOAT8(result);
}
/*
* dsind - returns the sine of arg1 (degrees)
*/
Datum
dsind(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
int sign = 1;
/*
* Per the POSIX spec, return NaN if the input is NaN and throw an error
* if the input is infinite.
*/
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
if (isinf(arg1))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
INIT_DEGREE_CONSTANTS();
/* Reduce the range of the input to [0,90] degrees */
arg1 = fmod(arg1, 360.0);
if (arg1 < 0.0)
{
/* sind(-x) = -sind(x) */
arg1 = -arg1;
sign = -sign;
}
if (arg1 > 180.0)
{
/* sind(360-x) = -sind(x) */
arg1 = 360.0 - arg1;
sign = -sign;
}
if (arg1 > 90.0)
{
/* sind(180-x) = sind(x) */
arg1 = 180.0 - arg1;
}
result = sign * sind_q1(arg1);
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* dtand - returns the tangent of arg1 (degrees)
*/
Datum
dtand(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
volatile float8 tan_arg1;
int sign = 1;
/*
* Per the POSIX spec, return NaN if the input is NaN and throw an error
* if the input is infinite.
*/
if (isnan(arg1))
PG_RETURN_FLOAT8(get_float8_nan());
if (isinf(arg1))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
INIT_DEGREE_CONSTANTS();
/* Reduce the range of the input to [0,90] degrees */
arg1 = fmod(arg1, 360.0);
if (arg1 < 0.0)
{
/* tand(-x) = -tand(x) */
arg1 = -arg1;
sign = -sign;
}
if (arg1 > 180.0)
{
/* tand(360-x) = -tand(x) */
arg1 = 360.0 - arg1;
sign = -sign;
}
if (arg1 > 90.0)
{
/* tand(180-x) = -tand(x) */
arg1 = 180.0 - arg1;
sign = -sign;
}
tan_arg1 = sind_q1(arg1) / cosd_q1(arg1);
result = sign * (tan_arg1 / tan_45);
/*
* On some machines we get tand(180) = minus zero, but this isn't always
* true. For portability, and because the user constituency for this
* function probably doesn't want minus zero, force it to plain zero.
*/
if (result == 0.0)
result = 0.0;
check_float8_val(result, true /* tand(90) == Inf */ , true);
PG_RETURN_FLOAT8(result);
}
/*
* degrees - returns degrees converted from radians
*/
Datum
degrees(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
PG_RETURN_FLOAT8(float8_div(arg1, RADIANS_PER_DEGREE));
}
/*
* dpi - returns the constant PI
*/
Datum
dpi(PG_FUNCTION_ARGS)
{
PG_RETURN_FLOAT8(M_PI);
}
/*
* radians - returns radians converted from degrees
*/
Datum
radians(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
PG_RETURN_FLOAT8(float8_mul(arg1, RADIANS_PER_DEGREE));
}
/* ========== HYPERBOLIC FUNCTIONS ========== */
/*
* dsinh - returns the hyperbolic sine of arg1
*/
Datum
dsinh(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
errno = 0;
result = sinh(arg1);
/*
* if an ERANGE error occurs, it means there is an overflow. For sinh,
* the result should be either -infinity or infinity, depending on the
* sign of arg1.
*/
if (errno == ERANGE)
{
if (arg1 < 0)
result = -get_float8_infinity();
else
result = get_float8_infinity();
}
check_float8_val(result, true, true);
PG_RETURN_FLOAT8(result);
}
/*
* dcosh - returns the hyperbolic cosine of arg1
*/
Datum
dcosh(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
errno = 0;
result = cosh(arg1);
/*
* if an ERANGE error occurs, it means there is an overflow. As cosh is
* always positive, it always means the result is positive infinity.
*/
if (errno == ERANGE)
result = get_float8_infinity();
check_float8_val(result, true, false);
PG_RETURN_FLOAT8(result);
}
/*
* dtanh - returns the hyperbolic tangent of arg1
*/
Datum
dtanh(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/*
* For tanh, we don't need an errno check because it never overflows.
*/
result = tanh(arg1);
check_float8_val(result, false, true);
PG_RETURN_FLOAT8(result);
}
/*
* dasinh - returns the inverse hyperbolic sine of arg1
*/
Datum
dasinh(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/*
* For asinh, we don't need an errno check because it never overflows.
*/
result = asinh(arg1);
check_float8_val(result, true, true);
PG_RETURN_FLOAT8(result);
}
/*
* dacosh - returns the inverse hyperbolic cosine of arg1
*/
Datum
dacosh(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/*
* acosh is only defined for inputs >= 1.0. By checking this ourselves,
* we need not worry about checking for an EDOM error, which is a good
* thing because some implementations will report that for NaN. Otherwise,
* no error is possible.
*/
if (arg1 < 1.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
result = acosh(arg1);
check_float8_val(result, true, true);
PG_RETURN_FLOAT8(result);
}
/*
* datanh - returns the inverse hyperbolic tangent of arg1
*/
Datum
datanh(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float8 result;
/*
* atanh is only defined for inputs between -1 and 1. By checking this
* ourselves, we need not worry about checking for an EDOM error, which is
* a good thing because some implementations will report that for NaN.
*/
if (arg1 < -1.0 || arg1 > 1.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("input is out of range")));
/*
* Also handle the infinity cases ourselves; this is helpful because old
* glibc versions may produce the wrong errno for this. All other inputs
* cannot produce an error.
*/
if (arg1 == -1.0)
result = -get_float8_infinity();
else if (arg1 == 1.0)
result = get_float8_infinity();
else
result = atanh(arg1);
check_float8_val(result, true, true);
PG_RETURN_FLOAT8(result);
}
/*
* drandom - returns a random number
*/
Datum
drandom(PG_FUNCTION_ARGS)
{
float8 result;
Use a separate random seed for SQL random()/setseed() functions. Previously, the SQL random() function depended on libc's random(3), and setseed() invoked srandom(3). This results in interference between these functions and backend-internal uses of random(3). We'd never paid too much mind to that, but in the wake of commit 88bdbd3f7 which added log_statement_sample_rate, the interference arguably has a security consequence: if log_statement_sample_rate is active then an unprivileged user could probably control which if any of his SQL commands get logged, by issuing setseed() at the right times. That seems bad. To fix this reliably, we need random() and setseed() to use their own private random state variable. Standard random(3) isn't amenable to such usage, so let's switch to pg_erand48(). It's hard to say whether that's more or less "random" than any particular platform's version of random(3), but it does have a wider seed value and a longer period than are required by POSIX, so we can hope that this isn't a big downgrade. Also, we should now have uniform behavior of random() across platforms, which is worth something. While at it, upgrade the per-process seed initialization method to use pg_strong_random() if available, greatly reducing the predictability of the initial seed value. (I'll separately do something similar for the internal uses of random().) In addition to forestalling the possible security problem, this has a benefit in the other direction, which is that we can now document setseed() as guaranteeing a reproducible sequence of random() values. Previously, because of the possibility of internal calls of random(3), we could not promise any such thing. Discussion: https://postgr.es/m/3859.1545849900@sss.pgh.pa.us
2018-12-29 23:33:27 +01:00
/* Initialize random seed, if not done yet in this process */
if (unlikely(!drandom_seed_set))
{
/*
* If possible, initialize the seed using high-quality random bits.
* Should that fail for some reason, we fall back on a lower-quality
* seed based on current time and PID.
*/
if (!pg_strong_random(drandom_seed, sizeof(drandom_seed)))
Use a separate random seed for SQL random()/setseed() functions. Previously, the SQL random() function depended on libc's random(3), and setseed() invoked srandom(3). This results in interference between these functions and backend-internal uses of random(3). We'd never paid too much mind to that, but in the wake of commit 88bdbd3f7 which added log_statement_sample_rate, the interference arguably has a security consequence: if log_statement_sample_rate is active then an unprivileged user could probably control which if any of his SQL commands get logged, by issuing setseed() at the right times. That seems bad. To fix this reliably, we need random() and setseed() to use their own private random state variable. Standard random(3) isn't amenable to such usage, so let's switch to pg_erand48(). It's hard to say whether that's more or less "random" than any particular platform's version of random(3), but it does have a wider seed value and a longer period than are required by POSIX, so we can hope that this isn't a big downgrade. Also, we should now have uniform behavior of random() across platforms, which is worth something. While at it, upgrade the per-process seed initialization method to use pg_strong_random() if available, greatly reducing the predictability of the initial seed value. (I'll separately do something similar for the internal uses of random().) In addition to forestalling the possible security problem, this has a benefit in the other direction, which is that we can now document setseed() as guaranteeing a reproducible sequence of random() values. Previously, because of the possibility of internal calls of random(3), we could not promise any such thing. Discussion: https://postgr.es/m/3859.1545849900@sss.pgh.pa.us
2018-12-29 23:33:27 +01:00
{
TimestampTz now = GetCurrentTimestamp();
uint64 iseed;
/* Mix the PID with the most predictable bits of the timestamp */
iseed = (uint64) now ^ ((uint64) MyProcPid << 32);
drandom_seed[0] = (unsigned short) iseed;
drandom_seed[1] = (unsigned short) (iseed >> 16);
drandom_seed[2] = (unsigned short) (iseed >> 32);
}
drandom_seed_set = true;
}
/* pg_erand48 produces desired result range [0.0 - 1.0) */
result = pg_erand48(drandom_seed);
PG_RETURN_FLOAT8(result);
}
/*
* setseed - set seed for the random number generator
*/
Datum
setseed(PG_FUNCTION_ARGS)
{
float8 seed = PG_GETARG_FLOAT8(0);
Use a separate random seed for SQL random()/setseed() functions. Previously, the SQL random() function depended on libc's random(3), and setseed() invoked srandom(3). This results in interference between these functions and backend-internal uses of random(3). We'd never paid too much mind to that, but in the wake of commit 88bdbd3f7 which added log_statement_sample_rate, the interference arguably has a security consequence: if log_statement_sample_rate is active then an unprivileged user could probably control which if any of his SQL commands get logged, by issuing setseed() at the right times. That seems bad. To fix this reliably, we need random() and setseed() to use their own private random state variable. Standard random(3) isn't amenable to such usage, so let's switch to pg_erand48(). It's hard to say whether that's more or less "random" than any particular platform's version of random(3), but it does have a wider seed value and a longer period than are required by POSIX, so we can hope that this isn't a big downgrade. Also, we should now have uniform behavior of random() across platforms, which is worth something. While at it, upgrade the per-process seed initialization method to use pg_strong_random() if available, greatly reducing the predictability of the initial seed value. (I'll separately do something similar for the internal uses of random().) In addition to forestalling the possible security problem, this has a benefit in the other direction, which is that we can now document setseed() as guaranteeing a reproducible sequence of random() values. Previously, because of the possibility of internal calls of random(3), we could not promise any such thing. Discussion: https://postgr.es/m/3859.1545849900@sss.pgh.pa.us
2018-12-29 23:33:27 +01:00
uint64 iseed;
Use a separate random seed for SQL random()/setseed() functions. Previously, the SQL random() function depended on libc's random(3), and setseed() invoked srandom(3). This results in interference between these functions and backend-internal uses of random(3). We'd never paid too much mind to that, but in the wake of commit 88bdbd3f7 which added log_statement_sample_rate, the interference arguably has a security consequence: if log_statement_sample_rate is active then an unprivileged user could probably control which if any of his SQL commands get logged, by issuing setseed() at the right times. That seems bad. To fix this reliably, we need random() and setseed() to use their own private random state variable. Standard random(3) isn't amenable to such usage, so let's switch to pg_erand48(). It's hard to say whether that's more or less "random" than any particular platform's version of random(3), but it does have a wider seed value and a longer period than are required by POSIX, so we can hope that this isn't a big downgrade. Also, we should now have uniform behavior of random() across platforms, which is worth something. While at it, upgrade the per-process seed initialization method to use pg_strong_random() if available, greatly reducing the predictability of the initial seed value. (I'll separately do something similar for the internal uses of random().) In addition to forestalling the possible security problem, this has a benefit in the other direction, which is that we can now document setseed() as guaranteeing a reproducible sequence of random() values. Previously, because of the possibility of internal calls of random(3), we could not promise any such thing. Discussion: https://postgr.es/m/3859.1545849900@sss.pgh.pa.us
2018-12-29 23:33:27 +01:00
if (seed < -1 || seed > 1 || isnan(seed))
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("setseed parameter %g is out of allowed range [-1,1]",
seed)));
/* Use sign bit + 47 fractional bits to fill drandom_seed[] */
iseed = (int64) (seed * (float8) UINT64CONST(0x7FFFFFFFFFFF));
drandom_seed[0] = (unsigned short) iseed;
drandom_seed[1] = (unsigned short) (iseed >> 16);
drandom_seed[2] = (unsigned short) (iseed >> 32);
drandom_seed_set = true;
PG_RETURN_VOID();
}
/*
* =========================
* FLOAT AGGREGATE OPERATORS
* =========================
*
* float8_accum - accumulate for AVG(), variance aggregates, etc.
* float4_accum - same, but input data is float4
* float8_avg - produce final result for float AVG()
* float8_var_samp - produce final result for float VAR_SAMP()
* float8_var_pop - produce final result for float VAR_POP()
* float8_stddev_samp - produce final result for float STDDEV_SAMP()
* float8_stddev_pop - produce final result for float STDDEV_POP()
*
* The naive schoolbook implementation of these aggregates works by
* accumulating sum(X) and sum(X^2). However, this approach suffers from
* large rounding errors in the final computation of quantities like the
* population variance (N*sum(X^2) - sum(X)^2) / N^2, since each of the
* intermediate terms is potentially very large, while the difference is often
* quite small.
*
* Instead we use the Youngs-Cramer algorithm [1] which works by accumulating
* Sx=sum(X) and Sxx=sum((X-Sx/N)^2), using a numerically stable algorithm to
* incrementally update those quantities. The final computations of each of
* the aggregate values is then trivial and gives more accurate results (for
* example, the population variance is just Sxx/N). This algorithm is also
* fairly easy to generalize to allow parallel execution without loss of
* precision (see, for example, [2]). For more details, and a comparison of
* this with other algorithms, see [3].
*
* The transition datatype for all these aggregates is a 3-element array
* of float8, holding the values N, Sx, Sxx in that order.
*
* Note that we represent N as a float to avoid having to build a special
* datatype. Given a reasonable floating-point implementation, there should
* be no accuracy loss unless N exceeds 2 ^ 52 or so (by which time the
* user will have doubtless lost interest anyway...)
*
* [1] Some Results Relevant to Choice of Sum and Sum-of-Product Algorithms,
* E. A. Youngs and E. M. Cramer, Technometrics Vol 13, No 3, August 1971.
*
* [2] Updating Formulae and a Pairwise Algorithm for Computing Sample
* Variances, T. F. Chan, G. H. Golub & R. J. LeVeque, COMPSTAT 1982.
*
* [3] Numerically Stable Parallel Computation of (Co-)Variance, Erich
* Schubert and Michael Gertz, Proceedings of the 30th International
* Conference on Scientific and Statistical Database Management, 2018.
*/
static float8 *
check_float8_array(ArrayType *transarray, const char *caller, int n)
{
/*
* We expect the input to be an N-element float array; verify that. We
2001-03-22 05:01:46 +01:00
* don't need to use deconstruct_array() since the array data is just
* going to look like a C array of N float8 values.
*/
if (ARR_NDIM(transarray) != 1 ||
ARR_DIMS(transarray)[0] != n ||
ARR_HASNULL(transarray) ||
ARR_ELEMTYPE(transarray) != FLOAT8OID)
elog(ERROR, "%s: expected %d-element float8 array", caller, n);
return (float8 *) ARR_DATA_PTR(transarray);
}
/*
* float8_combine
*
* An aggregate combine function used to combine two 3 fields
* aggregate transition data into a single transition data.
* This function is used only in two stage aggregation and
* shouldn't be called outside aggregate context.
*/
Datum
float8_combine(PG_FUNCTION_ARGS)
{
ArrayType *transarray1 = PG_GETARG_ARRAYTYPE_P(0);
ArrayType *transarray2 = PG_GETARG_ARRAYTYPE_P(1);
float8 *transvalues1;
float8 *transvalues2;
float8 N1,
Sx1,
Sxx1,
N2,
Sx2,
Sxx2,
tmp,
N,
Sx,
Sxx;
transvalues1 = check_float8_array(transarray1, "float8_combine", 3);
transvalues2 = check_float8_array(transarray2, "float8_combine", 3);
N1 = transvalues1[0];
Sx1 = transvalues1[1];
Sxx1 = transvalues1[2];
N2 = transvalues2[0];
Sx2 = transvalues2[1];
Sxx2 = transvalues2[2];
/*--------------------
* The transition values combine using a generalization of the
* Youngs-Cramer algorithm as follows:
*
* N = N1 + N2
* Sx = Sx1 + Sx2
* Sxx = Sxx1 + Sxx2 + N1 * N2 * (Sx1/N1 - Sx2/N2)^2 / N;
*
* It's worth handling the special cases N1 = 0 and N2 = 0 separately
* since those cases are trivial, and we then don't need to worry about
* division-by-zero errors in the general case.
*--------------------
*/
if (N1 == 0.0)
{
N = N2;
Sx = Sx2;
Sxx = Sxx2;
}
else if (N2 == 0.0)
{
N = N1;
Sx = Sx1;
Sxx = Sxx1;
}
else
{
N = N1 + N2;
Sx = float8_pl(Sx1, Sx2);
tmp = Sx1 / N1 - Sx2 / N2;
Sxx = Sxx1 + Sxx2 + N1 * N2 * tmp * tmp / N;
check_float8_val(Sxx, isinf(Sxx1) || isinf(Sxx2), true);
}
/*
* If we're invoked as an aggregate, we can cheat and modify our first
* parameter in-place to reduce palloc overhead. Otherwise we construct a
* new array with the updated transition data and return it.
*/
if (AggCheckCallContext(fcinfo, NULL))
{
transvalues1[0] = N;
transvalues1[1] = Sx;
transvalues1[2] = Sxx;
PG_RETURN_ARRAYTYPE_P(transarray1);
}
else
{
Datum transdatums[3];
ArrayType *result;
transdatums[0] = Float8GetDatumFast(N);
transdatums[1] = Float8GetDatumFast(Sx);
transdatums[2] = Float8GetDatumFast(Sxx);
result = construct_array(transdatums, 3,
FLOAT8OID,
sizeof(float8), FLOAT8PASSBYVAL, 'd');
PG_RETURN_ARRAYTYPE_P(result);
}
}
Datum
float8_accum(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 newval = PG_GETARG_FLOAT8(1);
float8 *transvalues;
float8 N,
Sx,
Sxx,
tmp;
transvalues = check_float8_array(transarray, "float8_accum", 3);
N = transvalues[0];
Sx = transvalues[1];
Sxx = transvalues[2];
/*
* Use the Youngs-Cramer algorithm to incorporate the new value into the
* transition values.
*/
N += 1.0;
Sx += newval;
if (transvalues[0] > 0.0)
{
tmp = newval * N - Sx;
Sxx += tmp * tmp / (N * transvalues[0]);
/*
* Overflow check. We only report an overflow error when finite
* inputs lead to infinite results. Note also that Sxx should be NaN
* if any of the inputs are infinite, so we intentionally prevent Sxx
* from becoming infinite.
*/
if (isinf(Sx) || isinf(Sxx))
{
if (!isinf(transvalues[1]) && !isinf(newval))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("value out of range: overflow")));
Sxx = get_float8_nan();
}
}
2007-11-15 22:14:46 +01:00
/*
* If we're invoked as an aggregate, we can cheat and modify our first
2005-10-15 04:49:52 +02:00
* parameter in-place to reduce palloc overhead. Otherwise we construct a
* new array with the updated transition data and return it.
*/
if (AggCheckCallContext(fcinfo, NULL))
{
transvalues[0] = N;
transvalues[1] = Sx;
transvalues[2] = Sxx;
PG_RETURN_ARRAYTYPE_P(transarray);
}
else
{
Datum transdatums[3];
ArrayType *result;
transdatums[0] = Float8GetDatumFast(N);
transdatums[1] = Float8GetDatumFast(Sx);
transdatums[2] = Float8GetDatumFast(Sxx);
result = construct_array(transdatums, 3,
FLOAT8OID,
sizeof(float8), FLOAT8PASSBYVAL, 'd');
PG_RETURN_ARRAYTYPE_P(result);
}
}
Datum
float4_accum(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
2007-11-15 22:14:46 +01:00
/* do computations as float8 */
float8 newval = PG_GETARG_FLOAT4(1);
float8 *transvalues;
float8 N,
Sx,
Sxx,
tmp;
transvalues = check_float8_array(transarray, "float4_accum", 3);
N = transvalues[0];
Sx = transvalues[1];
Sxx = transvalues[2];
/*
* Use the Youngs-Cramer algorithm to incorporate the new value into the
* transition values.
*/
N += 1.0;
Sx += newval;
if (transvalues[0] > 0.0)
{
tmp = newval * N - Sx;
Sxx += tmp * tmp / (N * transvalues[0]);
/*
* Overflow check. We only report an overflow error when finite
* inputs lead to infinite results. Note also that Sxx should be NaN
* if any of the inputs are infinite, so we intentionally prevent Sxx
* from becoming infinite.
*/
if (isinf(Sx) || isinf(Sxx))
{
if (!isinf(transvalues[1]) && !isinf(newval))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("value out of range: overflow")));
Sxx = get_float8_nan();
}
}
2007-11-15 22:14:46 +01:00
/*
* If we're invoked as an aggregate, we can cheat and modify our first
2005-10-15 04:49:52 +02:00
* parameter in-place to reduce palloc overhead. Otherwise we construct a
* new array with the updated transition data and return it.
*/
if (AggCheckCallContext(fcinfo, NULL))
{
transvalues[0] = N;
transvalues[1] = Sx;
transvalues[2] = Sxx;
PG_RETURN_ARRAYTYPE_P(transarray);
}
else
{
Datum transdatums[3];
ArrayType *result;
transdatums[0] = Float8GetDatumFast(N);
transdatums[1] = Float8GetDatumFast(Sx);
transdatums[2] = Float8GetDatumFast(Sxx);
result = construct_array(transdatums, 3,
FLOAT8OID,
sizeof(float8), FLOAT8PASSBYVAL, 'd');
PG_RETURN_ARRAYTYPE_P(result);
}
}
Datum
float8_avg(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
float8 N,
Sx;
transvalues = check_float8_array(transarray, "float8_avg", 3);
N = transvalues[0];
Sx = transvalues[1];
/* ignore Sxx */
/* SQL defines AVG of no values to be NULL */
if (N == 0.0)
PG_RETURN_NULL();
PG_RETURN_FLOAT8(Sx / N);
}
Datum
float8_var_pop(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
float8 N,
Sxx;
transvalues = check_float8_array(transarray, "float8_var_pop", 3);
N = transvalues[0];
/* ignore Sx */
Sxx = transvalues[2];
/* Population variance is undefined when N is 0, so return NULL */
if (N == 0.0)
PG_RETURN_NULL();
/* Note that Sxx is guaranteed to be non-negative */
PG_RETURN_FLOAT8(Sxx / N);
}
Datum
float8_var_samp(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
float8 N,
Sxx;
transvalues = check_float8_array(transarray, "float8_var_samp", 3);
N = transvalues[0];
/* ignore Sx */
Sxx = transvalues[2];
/* Sample variance is undefined when N is 0 or 1, so return NULL */
if (N <= 1.0)
PG_RETURN_NULL();
/* Note that Sxx is guaranteed to be non-negative */
PG_RETURN_FLOAT8(Sxx / (N - 1.0));
}
Datum
float8_stddev_pop(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
float8 N,
Sxx;
transvalues = check_float8_array(transarray, "float8_stddev_pop", 3);
N = transvalues[0];
/* ignore Sx */
Sxx = transvalues[2];
/* Population stddev is undefined when N is 0, so return NULL */
if (N == 0.0)
PG_RETURN_NULL();
/* Note that Sxx is guaranteed to be non-negative */
PG_RETURN_FLOAT8(sqrt(Sxx / N));
}
Datum
float8_stddev_samp(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
float8 N,
Sxx;
transvalues = check_float8_array(transarray, "float8_stddev_samp", 3);
N = transvalues[0];
/* ignore Sx */
Sxx = transvalues[2];
/* Sample stddev is undefined when N is 0 or 1, so return NULL */
if (N <= 1.0)
PG_RETURN_NULL();
/* Note that Sxx is guaranteed to be non-negative */
PG_RETURN_FLOAT8(sqrt(Sxx / (N - 1.0)));
}
/*
* =========================
* SQL2003 BINARY AGGREGATES
* =========================
*
* As with the preceding aggregates, we use the Youngs-Cramer algorithm to
* reduce rounding errors in the aggregate final functions.
*
* The transition datatype for all these aggregates is a 6-element array of
* float8, holding the values N, Sx=sum(X), Sxx=sum((X-Sx/N)^2), Sy=sum(Y),
* Syy=sum((Y-Sy/N)^2), Sxy=sum((X-Sx/N)*(Y-Sy/N)) in that order.
*
* Note that Y is the first argument to all these aggregates!
*
* It might seem attractive to optimize this by having multiple accumulator
* functions that only calculate the sums actually needed. But on most
* modern machines, a couple of extra floating-point multiplies will be
* insignificant compared to the other per-tuple overhead, so I've chosen
* to minimize code space instead.
*/
Datum
float8_regr_accum(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 newvalY = PG_GETARG_FLOAT8(1);
float8 newvalX = PG_GETARG_FLOAT8(2);
float8 *transvalues;
2006-10-04 02:30:14 +02:00
float8 N,
Sx,
Sxx,
Sy,
Syy,
Sxy,
tmpX,
tmpY,
scale;
transvalues = check_float8_array(transarray, "float8_regr_accum", 6);
N = transvalues[0];
Sx = transvalues[1];
Sxx = transvalues[2];
Sy = transvalues[3];
Syy = transvalues[4];
Sxy = transvalues[5];
/*
* Use the Youngs-Cramer algorithm to incorporate the new values into the
* transition values.
*/
N += 1.0;
Sx += newvalX;
Sy += newvalY;
if (transvalues[0] > 0.0)
{
tmpX = newvalX * N - Sx;
tmpY = newvalY * N - Sy;
scale = 1.0 / (N * transvalues[0]);
Sxx += tmpX * tmpX * scale;
Syy += tmpY * tmpY * scale;
Sxy += tmpX * tmpY * scale;
/*
* Overflow check. We only report an overflow error when finite
* inputs lead to infinite results. Note also that Sxx, Syy and Sxy
* should be NaN if any of the relevant inputs are infinite, so we
* intentionally prevent them from becoming infinite.
*/
if (isinf(Sx) || isinf(Sxx) || isinf(Sy) || isinf(Syy) || isinf(Sxy))
{
if (((isinf(Sx) || isinf(Sxx)) &&
!isinf(transvalues[1]) && !isinf(newvalX)) ||
((isinf(Sy) || isinf(Syy)) &&
!isinf(transvalues[3]) && !isinf(newvalY)) ||
(isinf(Sxy) &&
!isinf(transvalues[1]) && !isinf(newvalX) &&
!isinf(transvalues[3]) && !isinf(newvalY)))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("value out of range: overflow")));
if (isinf(Sxx))
Sxx = get_float8_nan();
if (isinf(Syy))
Syy = get_float8_nan();
if (isinf(Sxy))
Sxy = get_float8_nan();
}
}
2007-11-15 22:14:46 +01:00
/*
* If we're invoked as an aggregate, we can cheat and modify our first
* parameter in-place to reduce palloc overhead. Otherwise we construct a
* new array with the updated transition data and return it.
*/
if (AggCheckCallContext(fcinfo, NULL))
{
transvalues[0] = N;
transvalues[1] = Sx;
transvalues[2] = Sxx;
transvalues[3] = Sy;
transvalues[4] = Syy;
transvalues[5] = Sxy;
PG_RETURN_ARRAYTYPE_P(transarray);
}
else
{
Datum transdatums[6];
ArrayType *result;
transdatums[0] = Float8GetDatumFast(N);
transdatums[1] = Float8GetDatumFast(Sx);
transdatums[2] = Float8GetDatumFast(Sxx);
transdatums[3] = Float8GetDatumFast(Sy);
transdatums[4] = Float8GetDatumFast(Syy);
transdatums[5] = Float8GetDatumFast(Sxy);
result = construct_array(transdatums, 6,
FLOAT8OID,
sizeof(float8), FLOAT8PASSBYVAL, 'd');
PG_RETURN_ARRAYTYPE_P(result);
}
}
/*
* float8_regr_combine
*
* An aggregate combine function used to combine two 6 fields
* aggregate transition data into a single transition data.
* This function is used only in two stage aggregation and
* shouldn't be called outside aggregate context.
*/
Datum
float8_regr_combine(PG_FUNCTION_ARGS)
{
ArrayType *transarray1 = PG_GETARG_ARRAYTYPE_P(0);
ArrayType *transarray2 = PG_GETARG_ARRAYTYPE_P(1);
float8 *transvalues1;
float8 *transvalues2;
float8 N1,
Sx1,
Sxx1,
Sy1,
Syy1,
Sxy1,
N2,
Sx2,
Sxx2,
Sy2,
Syy2,
Sxy2,
tmp1,
tmp2,
N,
Sx,
Sxx,
Sy,
Syy,
Sxy;
transvalues1 = check_float8_array(transarray1, "float8_regr_combine", 6);
transvalues2 = check_float8_array(transarray2, "float8_regr_combine", 6);
N1 = transvalues1[0];
Sx1 = transvalues1[1];
Sxx1 = transvalues1[2];
Sy1 = transvalues1[3];
Syy1 = transvalues1[4];
Sxy1 = transvalues1[5];
N2 = transvalues2[0];
Sx2 = transvalues2[1];
Sxx2 = transvalues2[2];
Sy2 = transvalues2[3];
Syy2 = transvalues2[4];
Sxy2 = transvalues2[5];
/*--------------------
* The transition values combine using a generalization of the
* Youngs-Cramer algorithm as follows:
*
* N = N1 + N2
* Sx = Sx1 + Sx2
* Sxx = Sxx1 + Sxx2 + N1 * N2 * (Sx1/N1 - Sx2/N2)^2 / N
* Sy = Sy1 + Sy2
* Syy = Syy1 + Syy2 + N1 * N2 * (Sy1/N1 - Sy2/N2)^2 / N
* Sxy = Sxy1 + Sxy2 + N1 * N2 * (Sx1/N1 - Sx2/N2) * (Sy1/N1 - Sy2/N2) / N
*
* It's worth handling the special cases N1 = 0 and N2 = 0 separately
* since those cases are trivial, and we then don't need to worry about
* division-by-zero errors in the general case.
*--------------------
*/
if (N1 == 0.0)
{
N = N2;
Sx = Sx2;
Sxx = Sxx2;
Sy = Sy2;
Syy = Syy2;
Sxy = Sxy2;
}
else if (N2 == 0.0)
{
N = N1;
Sx = Sx1;
Sxx = Sxx1;
Sy = Sy1;
Syy = Syy1;
Sxy = Sxy1;
}
else
{
N = N1 + N2;
Sx = float8_pl(Sx1, Sx2);
tmp1 = Sx1 / N1 - Sx2 / N2;
Sxx = Sxx1 + Sxx2 + N1 * N2 * tmp1 * tmp1 / N;
check_float8_val(Sxx, isinf(Sxx1) || isinf(Sxx2), true);
Sy = float8_pl(Sy1, Sy2);
tmp2 = Sy1 / N1 - Sy2 / N2;
Syy = Syy1 + Syy2 + N1 * N2 * tmp2 * tmp2 / N;
check_float8_val(Syy, isinf(Syy1) || isinf(Syy2), true);
Sxy = Sxy1 + Sxy2 + N1 * N2 * tmp1 * tmp2 / N;
check_float8_val(Sxy, isinf(Sxy1) || isinf(Sxy2), true);
}
/*
* If we're invoked as an aggregate, we can cheat and modify our first
* parameter in-place to reduce palloc overhead. Otherwise we construct a
* new array with the updated transition data and return it.
*/
if (AggCheckCallContext(fcinfo, NULL))
{
transvalues1[0] = N;
transvalues1[1] = Sx;
transvalues1[2] = Sxx;
transvalues1[3] = Sy;
transvalues1[4] = Syy;
transvalues1[5] = Sxy;
PG_RETURN_ARRAYTYPE_P(transarray1);
}
else
{
Datum transdatums[6];
ArrayType *result;
transdatums[0] = Float8GetDatumFast(N);
transdatums[1] = Float8GetDatumFast(Sx);
transdatums[2] = Float8GetDatumFast(Sxx);
transdatums[3] = Float8GetDatumFast(Sy);
transdatums[4] = Float8GetDatumFast(Syy);
transdatums[5] = Float8GetDatumFast(Sxy);
result = construct_array(transdatums, 6,
FLOAT8OID,
sizeof(float8), FLOAT8PASSBYVAL, 'd');
PG_RETURN_ARRAYTYPE_P(result);
}
}
Datum
float8_regr_sxx(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
float8 N,
Sxx;
transvalues = check_float8_array(transarray, "float8_regr_sxx", 6);
N = transvalues[0];
Sxx = transvalues[2];
/* if N is 0 we should return NULL */
if (N < 1.0)
PG_RETURN_NULL();
/* Note that Sxx is guaranteed to be non-negative */
PG_RETURN_FLOAT8(Sxx);
}
Datum
float8_regr_syy(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
float8 N,
Syy;
transvalues = check_float8_array(transarray, "float8_regr_syy", 6);
N = transvalues[0];
Syy = transvalues[4];
/* if N is 0 we should return NULL */
if (N < 1.0)
PG_RETURN_NULL();
/* Note that Syy is guaranteed to be non-negative */
PG_RETURN_FLOAT8(Syy);
}
Datum
float8_regr_sxy(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
2006-10-04 02:30:14 +02:00
float8 N,
Sxy;
transvalues = check_float8_array(transarray, "float8_regr_sxy", 6);
N = transvalues[0];
Sxy = transvalues[5];
/* if N is 0 we should return NULL */
if (N < 1.0)
PG_RETURN_NULL();
/* A negative result is valid here */
PG_RETURN_FLOAT8(Sxy);
}
Datum
float8_regr_avgx(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
float8 N,
Sx;
transvalues = check_float8_array(transarray, "float8_regr_avgx", 6);
N = transvalues[0];
Sx = transvalues[1];
/* if N is 0 we should return NULL */
if (N < 1.0)
PG_RETURN_NULL();
PG_RETURN_FLOAT8(Sx / N);
}
Datum
float8_regr_avgy(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
float8 N,
Sy;
transvalues = check_float8_array(transarray, "float8_regr_avgy", 6);
N = transvalues[0];
Sy = transvalues[3];
/* if N is 0 we should return NULL */
if (N < 1.0)
PG_RETURN_NULL();
PG_RETURN_FLOAT8(Sy / N);
}
Datum
float8_covar_pop(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
2006-10-04 02:30:14 +02:00
float8 N,
Sxy;
transvalues = check_float8_array(transarray, "float8_covar_pop", 6);
N = transvalues[0];
Sxy = transvalues[5];
/* if N is 0 we should return NULL */
if (N < 1.0)
PG_RETURN_NULL();
PG_RETURN_FLOAT8(Sxy / N);
}
Datum
float8_covar_samp(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
2006-10-04 02:30:14 +02:00
float8 N,
Sxy;
transvalues = check_float8_array(transarray, "float8_covar_samp", 6);
N = transvalues[0];
Sxy = transvalues[5];
/* if N is <= 1 we should return NULL */
if (N < 2.0)
PG_RETURN_NULL();
PG_RETURN_FLOAT8(Sxy / (N - 1.0));
}
Datum
float8_corr(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
2006-10-04 02:30:14 +02:00
float8 N,
Sxx,
Syy,
Sxy;
transvalues = check_float8_array(transarray, "float8_corr", 6);
N = transvalues[0];
Sxx = transvalues[2];
Syy = transvalues[4];
Sxy = transvalues[5];
/* if N is 0 we should return NULL */
if (N < 1.0)
PG_RETURN_NULL();
/* Note that Sxx and Syy are guaranteed to be non-negative */
/* per spec, return NULL for horizontal and vertical lines */
if (Sxx == 0 || Syy == 0)
PG_RETURN_NULL();
2006-10-04 02:30:14 +02:00
PG_RETURN_FLOAT8(Sxy / sqrt(Sxx * Syy));
}
Datum
float8_regr_r2(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
2006-10-04 02:30:14 +02:00
float8 N,
Sxx,
Syy,
Sxy;
transvalues = check_float8_array(transarray, "float8_regr_r2", 6);
N = transvalues[0];
Sxx = transvalues[2];
Syy = transvalues[4];
Sxy = transvalues[5];
/* if N is 0 we should return NULL */
if (N < 1.0)
PG_RETURN_NULL();
/* Note that Sxx and Syy are guaranteed to be non-negative */
/* per spec, return NULL for a vertical line */
if (Sxx == 0)
PG_RETURN_NULL();
/* per spec, return 1.0 for a horizontal line */
if (Syy == 0)
PG_RETURN_FLOAT8(1.0);
PG_RETURN_FLOAT8((Sxy * Sxy) / (Sxx * Syy));
}
Datum
float8_regr_slope(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
2006-10-04 02:30:14 +02:00
float8 N,
Sxx,
Sxy;
transvalues = check_float8_array(transarray, "float8_regr_slope", 6);
N = transvalues[0];
Sxx = transvalues[2];
Sxy = transvalues[5];
/* if N is 0 we should return NULL */
if (N < 1.0)
PG_RETURN_NULL();
/* Note that Sxx is guaranteed to be non-negative */
/* per spec, return NULL for a vertical line */
if (Sxx == 0)
PG_RETURN_NULL();
2006-10-04 02:30:14 +02:00
PG_RETURN_FLOAT8(Sxy / Sxx);
}
Datum
float8_regr_intercept(PG_FUNCTION_ARGS)
{
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
float8 *transvalues;
2006-10-04 02:30:14 +02:00
float8 N,
Sx,
Sxx,
Sy,
Sxy;
transvalues = check_float8_array(transarray, "float8_regr_intercept", 6);
N = transvalues[0];
Sx = transvalues[1];
Sxx = transvalues[2];
Sy = transvalues[3];
Sxy = transvalues[5];
/* if N is 0 we should return NULL */
if (N < 1.0)
PG_RETURN_NULL();
/* Note that Sxx is guaranteed to be non-negative */
/* per spec, return NULL for a vertical line */
if (Sxx == 0)
PG_RETURN_NULL();
2006-10-04 02:30:14 +02:00
PG_RETURN_FLOAT8((Sy - Sx * Sxy / Sxx) / N);
}
/*
* ====================================
* MIXED-PRECISION ARITHMETIC OPERATORS
* ====================================
*/
/*
* float48pl - returns arg1 + arg2
* float48mi - returns arg1 - arg2
* float48mul - returns arg1 * arg2
* float48div - returns arg1 / arg2
*/
Datum
float48pl(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_FLOAT8(float8_pl((float8) arg1, arg2));
}
Datum
float48mi(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_FLOAT8(float8_mi((float8) arg1, arg2));
}
Datum
float48mul(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_FLOAT8(float8_mul((float8) arg1, arg2));
}
Datum
float48div(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_FLOAT8(float8_div((float8) arg1, arg2));
}
/*
* float84pl - returns arg1 + arg2
* float84mi - returns arg1 - arg2
* float84mul - returns arg1 * arg2
* float84div - returns arg1 / arg2
*/
Datum
float84pl(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_FLOAT8(float8_pl(arg1, (float8) arg2));
}
Datum
float84mi(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_FLOAT8(float8_mi(arg1, (float8) arg2));
}
Datum
float84mul(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_FLOAT8(float8_mul(arg1, (float8) arg2));
}
Datum
float84div(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_FLOAT8(float8_div(arg1, (float8) arg2));
}
/*
* ====================
* COMPARISON OPERATORS
* ====================
*/
/*
* float48{eq,ne,lt,le,gt,ge} - float4/float8 comparison operations
*/
Datum
float48eq(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_eq((float8) arg1, arg2));
}
Datum
float48ne(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_ne((float8) arg1, arg2));
}
Datum
float48lt(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_lt((float8) arg1, arg2));
}
Datum
float48le(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_le((float8) arg1, arg2));
}
Datum
float48gt(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_gt((float8) arg1, arg2));
}
Datum
float48ge(PG_FUNCTION_ARGS)
{
float4 arg1 = PG_GETARG_FLOAT4(0);
float8 arg2 = PG_GETARG_FLOAT8(1);
PG_RETURN_BOOL(float8_ge((float8) arg1, arg2));
}
/*
* float84{eq,ne,lt,le,gt,ge} - float8/float4 comparison operations
*/
Datum
float84eq(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float8_eq(arg1, (float8) arg2));
}
Datum
float84ne(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float8_ne(arg1, (float8) arg2));
}
Datum
float84lt(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float8_lt(arg1, (float8) arg2));
}
Datum
float84le(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float8_le(arg1, (float8) arg2));
}
Datum
float84gt(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float8_gt(arg1, (float8) arg2));
}
Datum
float84ge(PG_FUNCTION_ARGS)
{
float8 arg1 = PG_GETARG_FLOAT8(0);
float4 arg2 = PG_GETARG_FLOAT4(1);
PG_RETURN_BOOL(float8_ge(arg1, (float8) arg2));
}
/*
* Implements the float8 version of the width_bucket() function
* defined by SQL2003. See also width_bucket_numeric().
*
* 'bound1' and 'bound2' are the lower and upper bounds of the
* histogram's range, respectively. 'count' is the number of buckets
* in the histogram. width_bucket() returns an integer indicating the
* bucket number that 'operand' belongs to in an equiwidth histogram
* with the specified characteristics. An operand smaller than the
* lower bound is assigned to bucket 0. An operand greater than the
* upper bound is assigned to an additional bucket (with number
* count+1). We don't allow "NaN" for any of the float8 inputs, and we
* don't allow either of the histogram bounds to be +/- infinity.
*/
Datum
width_bucket_float8(PG_FUNCTION_ARGS)
{
2007-11-15 22:14:46 +01:00
float8 operand = PG_GETARG_FLOAT8(0);
float8 bound1 = PG_GETARG_FLOAT8(1);
float8 bound2 = PG_GETARG_FLOAT8(2);
int32 count = PG_GETARG_INT32(3);
int32 result;
if (count <= 0.0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
errmsg("count must be greater than zero")));
if (isnan(operand) || isnan(bound1) || isnan(bound2))
ereport(ERROR,
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
errmsg("operand, lower bound, and upper bound cannot be NaN")));
/* Note that we allow "operand" to be infinite */
if (isinf(bound1) || isinf(bound2))
ereport(ERROR,
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
errmsg("lower and upper bounds must be finite")));
if (bound1 < bound2)
{
if (operand < bound1)
result = 0;
else if (operand >= bound2)
{
if (pg_add_s32_overflow(count, 1, &result))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("integer out of range")));
}
else
result = ((float8) count * (operand - bound1) / (bound2 - bound1)) + 1;
}
else if (bound1 > bound2)
{
if (operand > bound1)
result = 0;
else if (operand <= bound2)
{
if (pg_add_s32_overflow(count, 1, &result))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("integer out of range")));
}
else
result = ((float8) count * (bound1 - operand) / (bound1 - bound2)) + 1;
}
else
{
ereport(ERROR,
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
errmsg("lower bound cannot equal upper bound")));
2007-11-15 22:14:46 +01:00
result = 0; /* keep the compiler quiet */
}
PG_RETURN_INT32(result);
}
/* ========== PRIVATE ROUTINES ========== */
#ifndef HAVE_CBRT
static double
cbrt(double x)
{
int isneg = (x < 0.0);
double absx = fabs(x);
double tmpres = pow(absx, (double) 1.0 / (double) 3.0);
/*
2006-10-04 02:30:14 +02:00
* The result is somewhat inaccurate --- not really pow()'s fault, as the
* exponent it's handed contains roundoff error. We can improve the
* accuracy by doing one iteration of Newton's formula. Beware of zero
* input however.
*/
if (tmpres > 0.0)
2006-10-04 02:30:14 +02:00
tmpres -= (tmpres - absx / (tmpres * tmpres)) / (double) 3.0;
1998-09-01 05:29:17 +02:00
return isneg ? -tmpres : tmpres;
}
Phase 2 of pgindent updates. Change pg_bsd_indent to follow upstream rules for placement of comments to the right of code, and remove pgindent hack that caused comments following #endif to not obey the general rule. Commit e3860ffa4dd0dad0dd9eea4be9cc1412373a8c89 wasn't actually using the published version of pg_bsd_indent, but a hacked-up version that tried to minimize the amount of movement of comments to the right of code. The situation of interest is where such a comment has to be moved to the right of its default placement at column 33 because there's code there. BSD indent has always moved right in units of tab stops in such cases --- but in the previous incarnation, indent was working in 8-space tab stops, while now it knows we use 4-space tabs. So the net result is that in about half the cases, such comments are placed one tab stop left of before. This is better all around: it leaves more room on the line for comment text, and it means that in such cases the comment uniformly starts at the next 4-space tab stop after the code, rather than sometimes one and sometimes two tabs after. Also, ensure that comments following #endif are indented the same as comments following other preprocessor commands such as #else. That inconsistency turns out to have been self-inflicted damage from a poorly-thought-through post-indent "fixup" in pgindent. This patch is much less interesting than the first round of indent changes, but also bulkier, so I thought it best to separate the effects. Discussion: https://postgr.es/m/E1dAmxK-0006EE-1r@gemulon.postgresql.org Discussion: https://postgr.es/m/30527.1495162840@sss.pgh.pa.us
2017-06-21 21:18:54 +02:00
#endif /* !HAVE_CBRT */