Line data Source code
1 : /*-------------------------------------------------------------------------
2 : *
3 : * float.c
4 : * Functions for the built-in floating-point types.
5 : *
6 : * Portions Copyright (c) 1996-2026, PostgreSQL Global Development Group
7 : * Portions Copyright (c) 1994, Regents of the University of California
8 : *
9 : *
10 : * IDENTIFICATION
11 : * src/backend/utils/adt/float.c
12 : *
13 : *-------------------------------------------------------------------------
14 : */
15 : #include "postgres.h"
16 :
17 : #include <ctype.h>
18 : #include <float.h>
19 : #include <math.h>
20 : #include <limits.h>
21 :
22 : #include "catalog/pg_type.h"
23 : #include "common/int.h"
24 : #include "common/shortest_dec.h"
25 : #include "libpq/pqformat.h"
26 : #include "utils/array.h"
27 : #include "utils/float.h"
28 : #include "utils/fmgrprotos.h"
29 : #include "utils/sortsupport.h"
30 :
31 :
32 : /*
33 : * Reject building with gcc's -ffast-math switch. It breaks our handling of
34 : * float Infinity and NaN values (via -ffinite-math-only), causes results to
35 : * be less accurate than expected (via -funsafe-math-optimizations and
36 : * -fexcess-precision=fast), and causes some math error reports to be missed
37 : * (via -fno-math-errno). Unfortunately we can't easily detect cases where
38 : * those options were given individually, but this at least catches the most
39 : * obvious case.
40 : *
41 : * We test this only here, not in any header file, to allow extensions to use
42 : * -ffast-math if they need to. But the inline functions in float.h will
43 : * misbehave in such an extension, so its authors had better be careful.
44 : */
45 : #ifdef __FAST_MATH__
46 : #error -ffast-math is known to break this code
47 : #endif
48 :
49 : /*
50 : * Configurable GUC parameter
51 : *
52 : * If >0, use shortest-decimal format for output; this is both the default and
53 : * allows for compatibility with clients that explicitly set a value here to
54 : * get round-trip-accurate results. If 0 or less, then use the old, slow,
55 : * decimal rounding method.
56 : */
57 : int extra_float_digits = 1;
58 :
59 : /* Cached constants for degree-based trig functions */
60 : static bool degree_consts_set = false;
61 : static float8 sin_30 = 0;
62 : static float8 one_minus_cos_60 = 0;
63 : static float8 asin_0_5 = 0;
64 : static float8 acos_0_5 = 0;
65 : static float8 atan_1_0 = 0;
66 : static float8 tan_45 = 0;
67 : static float8 cot_45 = 0;
68 :
69 : /*
70 : * These are intentionally not static; don't "fix" them. They will never
71 : * be referenced by other files, much less changed; but we don't want the
72 : * compiler to know that, else it might try to precompute expressions
73 : * involving them. See comments for init_degree_constants().
74 : *
75 : * The additional extern declarations are to silence
76 : * -Wmissing-variable-declarations.
77 : */
78 : extern float8 degree_c_thirty;
79 : extern float8 degree_c_forty_five;
80 : extern float8 degree_c_sixty;
81 : extern float8 degree_c_one_half;
82 : extern float8 degree_c_one;
83 : float8 degree_c_thirty = 30.0;
84 : float8 degree_c_forty_five = 45.0;
85 : float8 degree_c_sixty = 60.0;
86 : float8 degree_c_one_half = 0.5;
87 : float8 degree_c_one = 1.0;
88 :
89 : /* Local function prototypes */
90 : static double sind_q1(double x);
91 : static double cosd_q1(double x);
92 : static void init_degree_constants(void);
93 :
94 :
95 : /*
96 : * We use these out-of-line ereport() calls to report float overflow,
97 : * underflow, and zero-divide, because following our usual practice of
98 : * repeating them at each call site would lead to a lot of code bloat.
99 : *
100 : * This does mean that you don't get a useful error location indicator.
101 : */
102 : pg_noinline void
103 56 : float_overflow_error(void)
104 : {
105 56 : ereport(ERROR,
106 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
107 : errmsg("value out of range: overflow")));
108 : }
109 :
110 : pg_noinline void
111 20 : float_underflow_error(void)
112 : {
113 20 : ereport(ERROR,
114 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
115 : errmsg("value out of range: underflow")));
116 : }
117 :
118 : pg_noinline void
119 48 : float_zero_divide_error(void)
120 : {
121 48 : ereport(ERROR,
122 : (errcode(ERRCODE_DIVISION_BY_ZERO),
123 : errmsg("division by zero")));
124 : }
125 :
126 :
127 : /*
128 : * Returns -1 if 'val' represents negative infinity, 1 if 'val'
129 : * represents (positive) infinity, and 0 otherwise. On some platforms,
130 : * this is equivalent to the isinf() macro, but not everywhere: C99
131 : * does not specify that isinf() needs to distinguish between positive
132 : * and negative infinity.
133 : */
134 : int
135 0 : is_infinite(double val)
136 : {
137 0 : int inf = isinf(val);
138 :
139 0 : if (inf == 0)
140 0 : return 0;
141 0 : else if (val > 0)
142 0 : return 1;
143 : else
144 0 : return -1;
145 : }
146 :
147 :
148 : /* ========== USER I/O ROUTINES ========== */
149 :
150 :
151 : /*
152 : * float4in - converts "num" to float4
153 : *
154 : * Note that this code now uses strtof(), where it used to use strtod().
155 : *
156 : * The motivation for using strtof() is to avoid a double-rounding problem:
157 : * for certain decimal inputs, if you round the input correctly to a double,
158 : * and then round the double to a float, the result is incorrect in that it
159 : * does not match the result of rounding the decimal value to float directly.
160 : *
161 : * One of the best examples is 7.038531e-26:
162 : *
163 : * 0xAE43FDp-107 = 7.03853069185120912085...e-26
164 : * midpoint 7.03853100000000022281...e-26
165 : * 0xAE43FEp-107 = 7.03853130814879132477...e-26
166 : *
167 : * making 0xAE43FDp-107 the correct float result, but if you do the conversion
168 : * via a double, you get
169 : *
170 : * 0xAE43FD.7FFFFFF8p-107 = 7.03853099999999907487...e-26
171 : * midpoint 7.03853099999999964884...e-26
172 : * 0xAE43FD.80000000p-107 = 7.03853100000000022281...e-26
173 : * 0xAE43FD.80000008p-107 = 7.03853100000000137076...e-26
174 : *
175 : * so the value rounds to the double exactly on the midpoint between the two
176 : * nearest floats, and then rounding again to a float gives the incorrect
177 : * result of 0xAE43FEp-107.
178 : *
179 : */
180 : Datum
181 372490 : float4in(PG_FUNCTION_ARGS)
182 : {
183 372490 : char *num = PG_GETARG_CSTRING(0);
184 :
185 372490 : PG_RETURN_FLOAT4(float4in_internal(num, NULL, "real", num,
186 : fcinfo->context));
187 : }
188 :
189 : /*
190 : * float4in_internal - guts of float4in()
191 : *
192 : * This is exposed for use by functions that want a reasonably
193 : * platform-independent way of inputting floats. The behavior is
194 : * essentially like strtof + ereturn on error.
195 : *
196 : * Uses the same API as float8in_internal below, so most of its
197 : * comments also apply here, except regarding use in geometric types.
198 : */
199 : float4
200 377668 : float4in_internal(char *num, char **endptr_p,
201 : const char *type_name, const char *orig_string,
202 : struct Node *escontext)
203 : {
204 : float val;
205 : char *endptr;
206 :
207 : /*
208 : * endptr points to the first character _after_ the sequence we recognized
209 : * as a valid floating point number. orig_string points to the original
210 : * input string.
211 : */
212 :
213 : /* skip leading whitespace */
214 377808 : while (*num != '\0' && isspace((unsigned char) *num))
215 140 : num++;
216 :
217 : /*
218 : * Check for an empty-string input to begin with, to avoid the vagaries of
219 : * strtod() on different platforms.
220 : */
221 377668 : if (*num == '\0')
222 8 : ereturn(escontext, 0,
223 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
224 : errmsg("invalid input syntax for type %s: \"%s\"",
225 : type_name, orig_string)));
226 :
227 377660 : errno = 0;
228 377660 : val = strtof(num, &endptr);
229 :
230 : /* did we not see anything that looks like a double? */
231 377660 : if (endptr == num || errno != 0)
232 : {
233 112 : int save_errno = errno;
234 :
235 : /*
236 : * C99 requires that strtof() accept NaN, [+-]Infinity, and [+-]Inf,
237 : * but not all platforms support all of these (and some accept them
238 : * but set ERANGE anyway...) Therefore, we check for these inputs
239 : * ourselves if strtof() fails.
240 : *
241 : * Note: C99 also requires hexadecimal input as well as some extended
242 : * forms of NaN, but we consider these forms unportable and don't try
243 : * to support them. You can use 'em if your strtof() takes 'em.
244 : */
245 112 : if (pg_strncasecmp(num, "NaN", 3) == 0)
246 : {
247 0 : val = get_float4_nan();
248 0 : endptr = num + 3;
249 : }
250 112 : else if (pg_strncasecmp(num, "Infinity", 8) == 0)
251 : {
252 0 : val = get_float4_infinity();
253 0 : endptr = num + 8;
254 : }
255 112 : else if (pg_strncasecmp(num, "+Infinity", 9) == 0)
256 : {
257 0 : val = get_float4_infinity();
258 0 : endptr = num + 9;
259 : }
260 112 : else if (pg_strncasecmp(num, "-Infinity", 9) == 0)
261 : {
262 0 : val = -get_float4_infinity();
263 0 : endptr = num + 9;
264 : }
265 112 : else if (pg_strncasecmp(num, "inf", 3) == 0)
266 : {
267 0 : val = get_float4_infinity();
268 0 : endptr = num + 3;
269 : }
270 112 : else if (pg_strncasecmp(num, "+inf", 4) == 0)
271 : {
272 0 : val = get_float4_infinity();
273 0 : endptr = num + 4;
274 : }
275 112 : else if (pg_strncasecmp(num, "-inf", 4) == 0)
276 : {
277 0 : val = -get_float4_infinity();
278 0 : endptr = num + 4;
279 : }
280 112 : else if (save_errno == ERANGE)
281 : {
282 : /*
283 : * Some platforms return ERANGE for denormalized numbers (those
284 : * that are not zero, but are too close to zero to have full
285 : * precision). We'd prefer not to throw error for that, so try to
286 : * detect whether it's a "real" out-of-range condition by checking
287 : * to see if the result is zero or huge.
288 : */
289 47 : if (val == 0.0 ||
290 : #if !defined(HUGE_VALF)
291 : isinf(val)
292 : #else
293 12 : (val >= HUGE_VALF || val <= -HUGE_VALF)
294 : #endif
295 : )
296 : {
297 : /* see comments in float8in_internal for rationale */
298 43 : char *errnumber = pstrdup(num);
299 :
300 43 : errnumber[endptr - num] = '\0';
301 :
302 43 : ereturn(escontext, 0,
303 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
304 : errmsg("\"%s\" is out of range for type real",
305 : errnumber)));
306 : }
307 : }
308 : else
309 65 : ereturn(escontext, 0,
310 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
311 : errmsg("invalid input syntax for type %s: \"%s\"",
312 : type_name, orig_string)));
313 : }
314 :
315 : /* skip trailing whitespace */
316 377684 : while (*endptr != '\0' && isspace((unsigned char) *endptr))
317 132 : endptr++;
318 :
319 : /* report stopping point if wanted, else complain if not end of string */
320 377552 : if (endptr_p)
321 0 : *endptr_p = endptr;
322 377552 : else if (*endptr != '\0')
323 24 : ereturn(escontext, 0,
324 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
325 : errmsg("invalid input syntax for type %s: \"%s\"",
326 : type_name, orig_string)));
327 :
328 377528 : return val;
329 : }
330 :
331 : /*
332 : * float4out - converts a float4 number to a string
333 : * using a standard output format
334 : */
335 : Datum
336 190459 : float4out(PG_FUNCTION_ARGS)
337 : {
338 190459 : float4 num = PG_GETARG_FLOAT4(0);
339 190459 : char *ascii = (char *) palloc(32);
340 190459 : int ndig = FLT_DIG + extra_float_digits;
341 :
342 190459 : if (extra_float_digits > 0)
343 : {
344 184009 : float_to_shortest_decimal_buf(num, ascii);
345 184009 : PG_RETURN_CSTRING(ascii);
346 : }
347 :
348 6450 : (void) pg_strfromd(ascii, 32, ndig, num);
349 6450 : PG_RETURN_CSTRING(ascii);
350 : }
351 :
352 : /*
353 : * float4recv - converts external binary format to float4
354 : */
355 : Datum
356 0 : float4recv(PG_FUNCTION_ARGS)
357 : {
358 0 : StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
359 :
360 0 : PG_RETURN_FLOAT4(pq_getmsgfloat4(buf));
361 : }
362 :
363 : /*
364 : * float4send - converts float4 to binary format
365 : */
366 : Datum
367 4345 : float4send(PG_FUNCTION_ARGS)
368 : {
369 4345 : float4 num = PG_GETARG_FLOAT4(0);
370 : StringInfoData buf;
371 :
372 4345 : pq_begintypsend(&buf);
373 4345 : pq_sendfloat4(&buf, num);
374 4345 : PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
375 : }
376 :
377 : /*
378 : * float8in - converts "num" to float8
379 : */
380 : Datum
381 455456 : float8in(PG_FUNCTION_ARGS)
382 : {
383 455456 : char *num = PG_GETARG_CSTRING(0);
384 :
385 455456 : PG_RETURN_FLOAT8(float8in_internal(num, NULL, "double precision", num,
386 : fcinfo->context));
387 : }
388 :
389 : /*
390 : * float8in_internal - guts of float8in()
391 : *
392 : * This is exposed for use by functions that want a reasonably
393 : * platform-independent way of inputting doubles. The behavior is
394 : * essentially like strtod + ereturn on error, but note the following
395 : * differences:
396 : * 1. Both leading and trailing whitespace are skipped.
397 : * 2. If endptr_p is NULL, we report error if there's trailing junk.
398 : * Otherwise, it's up to the caller to complain about trailing junk.
399 : * 3. In event of a syntax error, the report mentions the given type_name
400 : * and prints orig_string as the input; this is meant to support use of
401 : * this function with types such as "box" and "point", where what we are
402 : * parsing here is just a substring of orig_string.
403 : *
404 : * If escontext points to an ErrorSaveContext node, that is filled instead
405 : * of throwing an error; the caller must check SOFT_ERROR_OCCURRED()
406 : * to detect errors.
407 : *
408 : * "num" could validly be declared "const char *", but that results in an
409 : * unreasonable amount of extra casting both here and in callers, so we don't.
410 : */
411 : float8
412 651343 : float8in_internal(char *num, char **endptr_p,
413 : const char *type_name, const char *orig_string,
414 : struct Node *escontext)
415 : {
416 : double val;
417 : char *endptr;
418 :
419 : /* skip leading whitespace */
420 652218 : while (*num != '\0' && isspace((unsigned char) *num))
421 875 : num++;
422 :
423 : /*
424 : * Check for an empty-string input to begin with, to avoid the vagaries of
425 : * strtod() on different platforms.
426 : */
427 651343 : if (*num == '\0')
428 12 : ereturn(escontext, 0,
429 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
430 : errmsg("invalid input syntax for type %s: \"%s\"",
431 : type_name, orig_string)));
432 :
433 651331 : errno = 0;
434 651331 : val = strtod(num, &endptr);
435 :
436 : /* did we not see anything that looks like a double? */
437 651331 : if (endptr == num || errno != 0)
438 : {
439 190 : int save_errno = errno;
440 :
441 : /*
442 : * C99 requires that strtod() accept NaN, [+-]Infinity, and [+-]Inf,
443 : * but not all platforms support all of these (and some accept them
444 : * but set ERANGE anyway...) Therefore, we check for these inputs
445 : * ourselves if strtod() fails.
446 : *
447 : * Note: C99 also requires hexadecimal input as well as some extended
448 : * forms of NaN, but we consider these forms unportable and don't try
449 : * to support them. You can use 'em if your strtod() takes 'em.
450 : */
451 190 : if (pg_strncasecmp(num, "NaN", 3) == 0)
452 : {
453 0 : val = get_float8_nan();
454 0 : endptr = num + 3;
455 : }
456 190 : else if (pg_strncasecmp(num, "Infinity", 8) == 0)
457 : {
458 0 : val = get_float8_infinity();
459 0 : endptr = num + 8;
460 : }
461 190 : else if (pg_strncasecmp(num, "+Infinity", 9) == 0)
462 : {
463 0 : val = get_float8_infinity();
464 0 : endptr = num + 9;
465 : }
466 190 : else if (pg_strncasecmp(num, "-Infinity", 9) == 0)
467 : {
468 0 : val = -get_float8_infinity();
469 0 : endptr = num + 9;
470 : }
471 190 : else if (pg_strncasecmp(num, "inf", 3) == 0)
472 : {
473 0 : val = get_float8_infinity();
474 0 : endptr = num + 3;
475 : }
476 190 : else if (pg_strncasecmp(num, "+inf", 4) == 0)
477 : {
478 0 : val = get_float8_infinity();
479 0 : endptr = num + 4;
480 : }
481 190 : else if (pg_strncasecmp(num, "-inf", 4) == 0)
482 : {
483 0 : val = -get_float8_infinity();
484 0 : endptr = num + 4;
485 : }
486 190 : else if (save_errno == ERANGE)
487 : {
488 : /*
489 : * Some platforms return ERANGE for denormalized numbers (those
490 : * that are not zero, but are too close to zero to have full
491 : * precision). We'd prefer not to throw error for that, so try to
492 : * detect whether it's a "real" out-of-range condition by checking
493 : * to see if the result is zero or huge.
494 : *
495 : * On error, we intentionally complain about double precision not
496 : * the given type name, and we print only the part of the string
497 : * that is the current number.
498 : */
499 86 : if (val == 0.0 || val >= HUGE_VAL || val <= -HUGE_VAL)
500 : {
501 71 : char *errnumber = pstrdup(num);
502 :
503 71 : errnumber[endptr - num] = '\0';
504 71 : ereturn(escontext, 0,
505 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
506 : errmsg("\"%s\" is out of range for type double precision",
507 : errnumber)));
508 : }
509 : }
510 : else
511 104 : ereturn(escontext, 0,
512 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
513 : errmsg("invalid input syntax for type %s: \"%s\"",
514 : type_name, orig_string)));
515 : }
516 :
517 : /* skip trailing whitespace */
518 651416 : while (*endptr != '\0' && isspace((unsigned char) *endptr))
519 260 : endptr++;
520 :
521 : /* report stopping point if wanted, else complain if not end of string */
522 651156 : if (endptr_p)
523 195555 : *endptr_p = endptr;
524 455601 : else if (*endptr != '\0')
525 28 : ereturn(escontext, 0,
526 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
527 : errmsg("invalid input syntax for type %s: \"%s\"",
528 : type_name, orig_string)));
529 :
530 651128 : return val;
531 : }
532 :
533 :
534 : /*
535 : * float8out - converts float8 number to a string
536 : * using a standard output format
537 : */
538 : Datum
539 608722 : float8out(PG_FUNCTION_ARGS)
540 : {
541 608722 : float8 num = PG_GETARG_FLOAT8(0);
542 :
543 608722 : PG_RETURN_CSTRING(float8out_internal(num));
544 : }
545 :
546 : /*
547 : * float8out_internal - guts of float8out()
548 : *
549 : * This is exposed for use by functions that want a reasonably
550 : * platform-independent way of outputting doubles.
551 : * The result is always palloc'd.
552 : */
553 : char *
554 2448896 : float8out_internal(double num)
555 : {
556 2448896 : char *ascii = (char *) palloc(32);
557 2448896 : int ndig = DBL_DIG + extra_float_digits;
558 :
559 2448896 : if (extra_float_digits > 0)
560 : {
561 2303916 : double_to_shortest_decimal_buf(num, ascii);
562 2303916 : return ascii;
563 : }
564 :
565 144980 : (void) pg_strfromd(ascii, 32, ndig, num);
566 144980 : return ascii;
567 : }
568 :
569 : /*
570 : * float8recv - converts external binary format to float8
571 : */
572 : Datum
573 16 : float8recv(PG_FUNCTION_ARGS)
574 : {
575 16 : StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
576 :
577 16 : PG_RETURN_FLOAT8(pq_getmsgfloat8(buf));
578 : }
579 :
580 : /*
581 : * float8send - converts float8 to binary format
582 : */
583 : Datum
584 3437 : float8send(PG_FUNCTION_ARGS)
585 : {
586 3437 : float8 num = PG_GETARG_FLOAT8(0);
587 : StringInfoData buf;
588 :
589 3437 : pq_begintypsend(&buf);
590 3437 : pq_sendfloat8(&buf, num);
591 3437 : PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
592 : }
593 :
594 :
595 : /* ========== PUBLIC ROUTINES ========== */
596 :
597 :
598 : /*
599 : * ======================
600 : * FLOAT4 BASE OPERATIONS
601 : * ======================
602 : */
603 :
604 : /*
605 : * float4abs - returns |arg1| (absolute value)
606 : */
607 : Datum
608 20 : float4abs(PG_FUNCTION_ARGS)
609 : {
610 20 : float4 arg1 = PG_GETARG_FLOAT4(0);
611 :
612 20 : PG_RETURN_FLOAT4(fabsf(arg1));
613 : }
614 :
615 : /*
616 : * float4um - returns -arg1 (unary minus)
617 : */
618 : Datum
619 10 : float4um(PG_FUNCTION_ARGS)
620 : {
621 10 : float4 arg1 = PG_GETARG_FLOAT4(0);
622 : float4 result;
623 :
624 10 : result = -arg1;
625 10 : PG_RETURN_FLOAT4(result);
626 : }
627 :
628 : Datum
629 0 : float4up(PG_FUNCTION_ARGS)
630 : {
631 0 : float4 arg = PG_GETARG_FLOAT4(0);
632 :
633 0 : PG_RETURN_FLOAT4(arg);
634 : }
635 :
636 : Datum
637 12 : float4larger(PG_FUNCTION_ARGS)
638 : {
639 12 : float4 arg1 = PG_GETARG_FLOAT4(0);
640 12 : float4 arg2 = PG_GETARG_FLOAT4(1);
641 : float4 result;
642 :
643 12 : if (float4_gt(arg1, arg2))
644 4 : result = arg1;
645 : else
646 8 : result = arg2;
647 12 : PG_RETURN_FLOAT4(result);
648 : }
649 :
650 : Datum
651 0 : float4smaller(PG_FUNCTION_ARGS)
652 : {
653 0 : float4 arg1 = PG_GETARG_FLOAT4(0);
654 0 : float4 arg2 = PG_GETARG_FLOAT4(1);
655 : float4 result;
656 :
657 0 : if (float4_lt(arg1, arg2))
658 0 : result = arg1;
659 : else
660 0 : result = arg2;
661 0 : PG_RETURN_FLOAT4(result);
662 : }
663 :
664 : /*
665 : * ======================
666 : * FLOAT8 BASE OPERATIONS
667 : * ======================
668 : */
669 :
670 : /*
671 : * float8abs - returns |arg1| (absolute value)
672 : */
673 : Datum
674 78241 : float8abs(PG_FUNCTION_ARGS)
675 : {
676 78241 : float8 arg1 = PG_GETARG_FLOAT8(0);
677 :
678 78241 : PG_RETURN_FLOAT8(fabs(arg1));
679 : }
680 :
681 :
682 : /*
683 : * float8um - returns -arg1 (unary minus)
684 : */
685 : Datum
686 211 : float8um(PG_FUNCTION_ARGS)
687 : {
688 211 : float8 arg1 = PG_GETARG_FLOAT8(0);
689 : float8 result;
690 :
691 211 : result = -arg1;
692 211 : PG_RETURN_FLOAT8(result);
693 : }
694 :
695 : Datum
696 0 : float8up(PG_FUNCTION_ARGS)
697 : {
698 0 : float8 arg = PG_GETARG_FLOAT8(0);
699 :
700 0 : PG_RETURN_FLOAT8(arg);
701 : }
702 :
703 : Datum
704 8584 : float8larger(PG_FUNCTION_ARGS)
705 : {
706 8584 : float8 arg1 = PG_GETARG_FLOAT8(0);
707 8584 : float8 arg2 = PG_GETARG_FLOAT8(1);
708 : float8 result;
709 :
710 8584 : if (float8_gt(arg1, arg2))
711 8184 : result = arg1;
712 : else
713 400 : result = arg2;
714 8584 : PG_RETURN_FLOAT8(result);
715 : }
716 :
717 : Datum
718 768 : float8smaller(PG_FUNCTION_ARGS)
719 : {
720 768 : float8 arg1 = PG_GETARG_FLOAT8(0);
721 768 : float8 arg2 = PG_GETARG_FLOAT8(1);
722 : float8 result;
723 :
724 768 : if (float8_lt(arg1, arg2))
725 592 : result = arg1;
726 : else
727 176 : result = arg2;
728 768 : PG_RETURN_FLOAT8(result);
729 : }
730 :
731 :
732 : /*
733 : * ====================
734 : * ARITHMETIC OPERATORS
735 : * ====================
736 : */
737 :
738 : /*
739 : * float4pl - returns arg1 + arg2
740 : * float4mi - returns arg1 - arg2
741 : * float4mul - returns arg1 * arg2
742 : * float4div - returns arg1 / arg2
743 : */
744 : Datum
745 36 : float4pl(PG_FUNCTION_ARGS)
746 : {
747 36 : float4 arg1 = PG_GETARG_FLOAT4(0);
748 36 : float4 arg2 = PG_GETARG_FLOAT4(1);
749 :
750 36 : PG_RETURN_FLOAT4(float4_pl(arg1, arg2));
751 : }
752 :
753 : Datum
754 12 : float4mi(PG_FUNCTION_ARGS)
755 : {
756 12 : float4 arg1 = PG_GETARG_FLOAT4(0);
757 12 : float4 arg2 = PG_GETARG_FLOAT4(1);
758 :
759 12 : PG_RETURN_FLOAT4(float4_mi(arg1, arg2));
760 : }
761 :
762 : Datum
763 24 : float4mul(PG_FUNCTION_ARGS)
764 : {
765 24 : float4 arg1 = PG_GETARG_FLOAT4(0);
766 24 : float4 arg2 = PG_GETARG_FLOAT4(1);
767 :
768 24 : PG_RETURN_FLOAT4(float4_mul(arg1, arg2));
769 : }
770 :
771 : Datum
772 36 : float4div(PG_FUNCTION_ARGS)
773 : {
774 36 : float4 arg1 = PG_GETARG_FLOAT4(0);
775 36 : float4 arg2 = PG_GETARG_FLOAT4(1);
776 :
777 36 : PG_RETURN_FLOAT4(float4_div(arg1, arg2));
778 : }
779 :
780 : /*
781 : * float8pl - returns arg1 + arg2
782 : * float8mi - returns arg1 - arg2
783 : * float8mul - returns arg1 * arg2
784 : * float8div - returns arg1 / arg2
785 : */
786 : Datum
787 65336 : float8pl(PG_FUNCTION_ARGS)
788 : {
789 65336 : float8 arg1 = PG_GETARG_FLOAT8(0);
790 65336 : float8 arg2 = PG_GETARG_FLOAT8(1);
791 :
792 65336 : PG_RETURN_FLOAT8(float8_pl(arg1, arg2));
793 : }
794 :
795 : Datum
796 8308 : float8mi(PG_FUNCTION_ARGS)
797 : {
798 8308 : float8 arg1 = PG_GETARG_FLOAT8(0);
799 8308 : float8 arg2 = PG_GETARG_FLOAT8(1);
800 :
801 8308 : PG_RETURN_FLOAT8(float8_mi(arg1, arg2));
802 : }
803 :
804 : Datum
805 937384 : float8mul(PG_FUNCTION_ARGS)
806 : {
807 937384 : float8 arg1 = PG_GETARG_FLOAT8(0);
808 937384 : float8 arg2 = PG_GETARG_FLOAT8(1);
809 :
810 937384 : PG_RETURN_FLOAT8(float8_mul(arg1, arg2));
811 : }
812 :
813 : Datum
814 10058 : float8div(PG_FUNCTION_ARGS)
815 : {
816 10058 : float8 arg1 = PG_GETARG_FLOAT8(0);
817 10058 : float8 arg2 = PG_GETARG_FLOAT8(1);
818 :
819 10058 : PG_RETURN_FLOAT8(float8_div(arg1, arg2));
820 : }
821 :
822 :
823 : /*
824 : * ====================
825 : * COMPARISON OPERATORS
826 : * ====================
827 : */
828 :
829 : /*
830 : * float4{eq,ne,lt,le,gt,ge} - float4/float4 comparison operations
831 : */
832 : int
833 8550296 : float4_cmp_internal(float4 a, float4 b)
834 : {
835 8550296 : if (float4_gt(a, b))
836 212535 : return 1;
837 8337761 : if (float4_lt(a, b))
838 1067954 : return -1;
839 7269807 : return 0;
840 : }
841 :
842 : Datum
843 24370 : float4eq(PG_FUNCTION_ARGS)
844 : {
845 24370 : float4 arg1 = PG_GETARG_FLOAT4(0);
846 24370 : float4 arg2 = PG_GETARG_FLOAT4(1);
847 :
848 24370 : PG_RETURN_BOOL(float4_eq(arg1, arg2));
849 : }
850 :
851 : Datum
852 20 : float4ne(PG_FUNCTION_ARGS)
853 : {
854 20 : float4 arg1 = PG_GETARG_FLOAT4(0);
855 20 : float4 arg2 = PG_GETARG_FLOAT4(1);
856 :
857 20 : PG_RETURN_BOOL(float4_ne(arg1, arg2));
858 : }
859 :
860 : Datum
861 9665 : float4lt(PG_FUNCTION_ARGS)
862 : {
863 9665 : float4 arg1 = PG_GETARG_FLOAT4(0);
864 9665 : float4 arg2 = PG_GETARG_FLOAT4(1);
865 :
866 9665 : PG_RETURN_BOOL(float4_lt(arg1, arg2));
867 : }
868 :
869 : Datum
870 2552 : float4le(PG_FUNCTION_ARGS)
871 : {
872 2552 : float4 arg1 = PG_GETARG_FLOAT4(0);
873 2552 : float4 arg2 = PG_GETARG_FLOAT4(1);
874 :
875 2552 : PG_RETURN_BOOL(float4_le(arg1, arg2));
876 : }
877 :
878 : Datum
879 3092 : float4gt(PG_FUNCTION_ARGS)
880 : {
881 3092 : float4 arg1 = PG_GETARG_FLOAT4(0);
882 3092 : float4 arg2 = PG_GETARG_FLOAT4(1);
883 :
884 3092 : PG_RETURN_BOOL(float4_gt(arg1, arg2));
885 : }
886 :
887 : Datum
888 2552 : float4ge(PG_FUNCTION_ARGS)
889 : {
890 2552 : float4 arg1 = PG_GETARG_FLOAT4(0);
891 2552 : float4 arg2 = PG_GETARG_FLOAT4(1);
892 :
893 2552 : PG_RETURN_BOOL(float4_ge(arg1, arg2));
894 : }
895 :
896 : Datum
897 953417 : btfloat4cmp(PG_FUNCTION_ARGS)
898 : {
899 953417 : float4 arg1 = PG_GETARG_FLOAT4(0);
900 953417 : float4 arg2 = PG_GETARG_FLOAT4(1);
901 :
902 953417 : PG_RETURN_INT32(float4_cmp_internal(arg1, arg2));
903 : }
904 :
905 : static int
906 7591848 : btfloat4fastcmp(Datum x, Datum y, SortSupport ssup)
907 : {
908 7591848 : float4 arg1 = DatumGetFloat4(x);
909 7591848 : float4 arg2 = DatumGetFloat4(y);
910 :
911 7591848 : return float4_cmp_internal(arg1, arg2);
912 : }
913 :
914 : Datum
915 669 : btfloat4sortsupport(PG_FUNCTION_ARGS)
916 : {
917 669 : SortSupport ssup = (SortSupport) PG_GETARG_POINTER(0);
918 :
919 669 : ssup->comparator = btfloat4fastcmp;
920 669 : PG_RETURN_VOID();
921 : }
922 :
923 : /*
924 : * float8{eq,ne,lt,le,gt,ge} - float8/float8 comparison operations
925 : */
926 : int
927 16405205 : float8_cmp_internal(float8 a, float8 b)
928 : {
929 16405205 : if (float8_gt(a, b))
930 6057675 : return 1;
931 10347530 : if (float8_lt(a, b))
932 10169200 : return -1;
933 178330 : return 0;
934 : }
935 :
936 : Datum
937 381590 : float8eq(PG_FUNCTION_ARGS)
938 : {
939 381590 : float8 arg1 = PG_GETARG_FLOAT8(0);
940 381590 : float8 arg2 = PG_GETARG_FLOAT8(1);
941 :
942 381590 : PG_RETURN_BOOL(float8_eq(arg1, arg2));
943 : }
944 :
945 : Datum
946 8243 : float8ne(PG_FUNCTION_ARGS)
947 : {
948 8243 : float8 arg1 = PG_GETARG_FLOAT8(0);
949 8243 : float8 arg2 = PG_GETARG_FLOAT8(1);
950 :
951 8243 : PG_RETURN_BOOL(float8_ne(arg1, arg2));
952 : }
953 :
954 : Datum
955 31824 : float8lt(PG_FUNCTION_ARGS)
956 : {
957 31824 : float8 arg1 = PG_GETARG_FLOAT8(0);
958 31824 : float8 arg2 = PG_GETARG_FLOAT8(1);
959 :
960 31824 : PG_RETURN_BOOL(float8_lt(arg1, arg2));
961 : }
962 :
963 : Datum
964 3810 : float8le(PG_FUNCTION_ARGS)
965 : {
966 3810 : float8 arg1 = PG_GETARG_FLOAT8(0);
967 3810 : float8 arg2 = PG_GETARG_FLOAT8(1);
968 :
969 3810 : PG_RETURN_BOOL(float8_le(arg1, arg2));
970 : }
971 :
972 : Datum
973 20316 : float8gt(PG_FUNCTION_ARGS)
974 : {
975 20316 : float8 arg1 = PG_GETARG_FLOAT8(0);
976 20316 : float8 arg2 = PG_GETARG_FLOAT8(1);
977 :
978 20316 : PG_RETURN_BOOL(float8_gt(arg1, arg2));
979 : }
980 :
981 : Datum
982 14022 : float8ge(PG_FUNCTION_ARGS)
983 : {
984 14022 : float8 arg1 = PG_GETARG_FLOAT8(0);
985 14022 : float8 arg2 = PG_GETARG_FLOAT8(1);
986 :
987 14022 : PG_RETURN_BOOL(float8_ge(arg1, arg2));
988 : }
989 :
990 : Datum
991 1982 : btfloat8cmp(PG_FUNCTION_ARGS)
992 : {
993 1982 : float8 arg1 = PG_GETARG_FLOAT8(0);
994 1982 : float8 arg2 = PG_GETARG_FLOAT8(1);
995 :
996 1982 : PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
997 : }
998 :
999 : static int
1000 4896067 : btfloat8fastcmp(Datum x, Datum y, SortSupport ssup)
1001 : {
1002 4896067 : float8 arg1 = DatumGetFloat8(x);
1003 4896067 : float8 arg2 = DatumGetFloat8(y);
1004 :
1005 4896067 : return float8_cmp_internal(arg1, arg2);
1006 : }
1007 :
1008 : Datum
1009 714 : btfloat8sortsupport(PG_FUNCTION_ARGS)
1010 : {
1011 714 : SortSupport ssup = (SortSupport) PG_GETARG_POINTER(0);
1012 :
1013 714 : ssup->comparator = btfloat8fastcmp;
1014 714 : PG_RETURN_VOID();
1015 : }
1016 :
1017 : Datum
1018 17 : btfloat48cmp(PG_FUNCTION_ARGS)
1019 : {
1020 17 : float4 arg1 = PG_GETARG_FLOAT4(0);
1021 17 : float8 arg2 = PG_GETARG_FLOAT8(1);
1022 :
1023 : /* widen float4 to float8 and then compare */
1024 17 : PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
1025 : }
1026 :
1027 : Datum
1028 149 : btfloat84cmp(PG_FUNCTION_ARGS)
1029 : {
1030 149 : float8 arg1 = PG_GETARG_FLOAT8(0);
1031 149 : float4 arg2 = PG_GETARG_FLOAT4(1);
1032 :
1033 : /* widen float4 to float8 and then compare */
1034 149 : PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
1035 : }
1036 :
1037 : /*
1038 : * in_range support function for float8.
1039 : *
1040 : * Note: we needn't supply a float8_float4 variant, as implicit coercion
1041 : * of the offset value takes care of that scenario just as well.
1042 : */
1043 : Datum
1044 768 : in_range_float8_float8(PG_FUNCTION_ARGS)
1045 : {
1046 768 : float8 val = PG_GETARG_FLOAT8(0);
1047 768 : float8 base = PG_GETARG_FLOAT8(1);
1048 768 : float8 offset = PG_GETARG_FLOAT8(2);
1049 768 : bool sub = PG_GETARG_BOOL(3);
1050 768 : bool less = PG_GETARG_BOOL(4);
1051 : float8 sum;
1052 :
1053 : /*
1054 : * Reject negative or NaN offset. Negative is per spec, and NaN is
1055 : * because appropriate semantics for that seem non-obvious.
1056 : */
1057 768 : if (isnan(offset) || offset < 0)
1058 4 : ereport(ERROR,
1059 : (errcode(ERRCODE_INVALID_PRECEDING_OR_FOLLOWING_SIZE),
1060 : errmsg("invalid preceding or following size in window function")));
1061 :
1062 : /*
1063 : * Deal with cases where val and/or base is NaN, following the rule that
1064 : * NaN sorts after non-NaN (cf float8_cmp_internal). The offset cannot
1065 : * affect the conclusion.
1066 : */
1067 764 : if (isnan(val))
1068 : {
1069 124 : if (isnan(base))
1070 40 : PG_RETURN_BOOL(true); /* NAN = NAN */
1071 : else
1072 84 : PG_RETURN_BOOL(!less); /* NAN > non-NAN */
1073 : }
1074 640 : else if (isnan(base))
1075 : {
1076 84 : PG_RETURN_BOOL(less); /* non-NAN < NAN */
1077 : }
1078 :
1079 : /*
1080 : * Deal with cases where both base and offset are infinite, and computing
1081 : * base +/- offset would produce NaN. This corresponds to a window frame
1082 : * whose boundary infinitely precedes +inf or infinitely follows -inf,
1083 : * which is not well-defined. For consistency with other cases involving
1084 : * infinities, such as the fact that +inf infinitely follows +inf, we
1085 : * choose to assume that +inf infinitely precedes +inf and -inf infinitely
1086 : * follows -inf, and therefore that all finite and infinite values are in
1087 : * such a window frame.
1088 : *
1089 : * offset is known positive, so we need only check the sign of base in
1090 : * this test.
1091 : */
1092 556 : if (isinf(offset) && isinf(base) &&
1093 : (sub ? base > 0 : base < 0))
1094 116 : PG_RETURN_BOOL(true);
1095 :
1096 : /*
1097 : * Otherwise it should be safe to compute base +/- offset. We trust the
1098 : * FPU to cope if an input is +/-inf or the true sum would overflow, and
1099 : * produce a suitably signed infinity, which will compare properly against
1100 : * val whether or not that's infinity.
1101 : */
1102 440 : if (sub)
1103 240 : sum = base - offset;
1104 : else
1105 200 : sum = base + offset;
1106 :
1107 440 : if (less)
1108 172 : PG_RETURN_BOOL(val <= sum);
1109 : else
1110 268 : PG_RETURN_BOOL(val >= sum);
1111 : }
1112 :
1113 : /*
1114 : * in_range support function for float4.
1115 : *
1116 : * We would need a float4_float8 variant in any case, so we supply that and
1117 : * let implicit coercion take care of the float4_float4 case.
1118 : */
1119 : Datum
1120 768 : in_range_float4_float8(PG_FUNCTION_ARGS)
1121 : {
1122 768 : float4 val = PG_GETARG_FLOAT4(0);
1123 768 : float4 base = PG_GETARG_FLOAT4(1);
1124 768 : float8 offset = PG_GETARG_FLOAT8(2);
1125 768 : bool sub = PG_GETARG_BOOL(3);
1126 768 : bool less = PG_GETARG_BOOL(4);
1127 : float8 sum;
1128 :
1129 : /*
1130 : * Reject negative or NaN offset. Negative is per spec, and NaN is
1131 : * because appropriate semantics for that seem non-obvious.
1132 : */
1133 768 : if (isnan(offset) || offset < 0)
1134 4 : ereport(ERROR,
1135 : (errcode(ERRCODE_INVALID_PRECEDING_OR_FOLLOWING_SIZE),
1136 : errmsg("invalid preceding or following size in window function")));
1137 :
1138 : /*
1139 : * Deal with cases where val and/or base is NaN, following the rule that
1140 : * NaN sorts after non-NaN (cf float8_cmp_internal). The offset cannot
1141 : * affect the conclusion.
1142 : */
1143 764 : if (isnan(val))
1144 : {
1145 124 : if (isnan(base))
1146 40 : PG_RETURN_BOOL(true); /* NAN = NAN */
1147 : else
1148 84 : PG_RETURN_BOOL(!less); /* NAN > non-NAN */
1149 : }
1150 640 : else if (isnan(base))
1151 : {
1152 84 : PG_RETURN_BOOL(less); /* non-NAN < NAN */
1153 : }
1154 :
1155 : /*
1156 : * Deal with cases where both base and offset are infinite, and computing
1157 : * base +/- offset would produce NaN. This corresponds to a window frame
1158 : * whose boundary infinitely precedes +inf or infinitely follows -inf,
1159 : * which is not well-defined. For consistency with other cases involving
1160 : * infinities, such as the fact that +inf infinitely follows +inf, we
1161 : * choose to assume that +inf infinitely precedes +inf and -inf infinitely
1162 : * follows -inf, and therefore that all finite and infinite values are in
1163 : * such a window frame.
1164 : *
1165 : * offset is known positive, so we need only check the sign of base in
1166 : * this test.
1167 : */
1168 556 : if (isinf(offset) && isinf(base) &&
1169 : (sub ? base > 0 : base < 0))
1170 116 : PG_RETURN_BOOL(true);
1171 :
1172 : /*
1173 : * Otherwise it should be safe to compute base +/- offset. We trust the
1174 : * FPU to cope if an input is +/-inf or the true sum would overflow, and
1175 : * produce a suitably signed infinity, which will compare properly against
1176 : * val whether or not that's infinity.
1177 : */
1178 440 : if (sub)
1179 240 : sum = base - offset;
1180 : else
1181 200 : sum = base + offset;
1182 :
1183 440 : if (less)
1184 172 : PG_RETURN_BOOL(val <= sum);
1185 : else
1186 268 : PG_RETURN_BOOL(val >= sum);
1187 : }
1188 :
1189 :
1190 : /*
1191 : * ===================
1192 : * CONVERSION ROUTINES
1193 : * ===================
1194 : */
1195 :
1196 : /*
1197 : * ftod - converts a float4 number to a float8 number
1198 : */
1199 : Datum
1200 201 : ftod(PG_FUNCTION_ARGS)
1201 : {
1202 201 : float4 num = PG_GETARG_FLOAT4(0);
1203 :
1204 201 : PG_RETURN_FLOAT8((float8) num);
1205 : }
1206 :
1207 :
1208 : /*
1209 : * dtof - converts a float8 number to a float4 number
1210 : */
1211 : Datum
1212 36 : dtof(PG_FUNCTION_ARGS)
1213 : {
1214 36 : float8 num = PG_GETARG_FLOAT8(0);
1215 : float4 result;
1216 :
1217 36 : result = (float4) num;
1218 36 : if (unlikely(isinf(result)) && !isinf(num))
1219 8 : float_overflow_error();
1220 28 : if (unlikely(result == 0.0f) && num != 0.0)
1221 8 : float_underflow_error();
1222 :
1223 20 : PG_RETURN_FLOAT4(result);
1224 : }
1225 :
1226 :
1227 : /*
1228 : * dtoi4 - converts a float8 number to an int4 number
1229 : */
1230 : Datum
1231 548175 : dtoi4(PG_FUNCTION_ARGS)
1232 : {
1233 548175 : float8 num = PG_GETARG_FLOAT8(0);
1234 :
1235 : /*
1236 : * Get rid of any fractional part in the input. This is so we don't fail
1237 : * on just-out-of-range values that would round into range. Note
1238 : * assumption that rint() will pass through a NaN or Inf unchanged.
1239 : */
1240 548175 : num = rint(num);
1241 :
1242 : /* Range check */
1243 548175 : if (unlikely(isnan(num) || !FLOAT8_FITS_IN_INT32(num)))
1244 16 : ereport(ERROR,
1245 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1246 : errmsg("integer out of range")));
1247 :
1248 548159 : PG_RETURN_INT32((int32) num);
1249 : }
1250 :
1251 :
1252 : /*
1253 : * dtoi2 - converts a float8 number to an int2 number
1254 : */
1255 : Datum
1256 62 : dtoi2(PG_FUNCTION_ARGS)
1257 : {
1258 62 : float8 num = PG_GETARG_FLOAT8(0);
1259 :
1260 : /*
1261 : * Get rid of any fractional part in the input. This is so we don't fail
1262 : * on just-out-of-range values that would round into range. Note
1263 : * assumption that rint() will pass through a NaN or Inf unchanged.
1264 : */
1265 62 : num = rint(num);
1266 :
1267 : /* Range check */
1268 62 : if (unlikely(isnan(num) || !FLOAT8_FITS_IN_INT16(num)))
1269 8 : ereport(ERROR,
1270 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1271 : errmsg("smallint out of range")));
1272 :
1273 54 : PG_RETURN_INT16((int16) num);
1274 : }
1275 :
1276 :
1277 : /*
1278 : * i4tod - converts an int4 number to a float8 number
1279 : */
1280 : Datum
1281 1533139 : i4tod(PG_FUNCTION_ARGS)
1282 : {
1283 1533139 : int32 num = PG_GETARG_INT32(0);
1284 :
1285 1533139 : PG_RETURN_FLOAT8((float8) num);
1286 : }
1287 :
1288 :
1289 : /*
1290 : * i2tod - converts an int2 number to a float8 number
1291 : */
1292 : Datum
1293 164 : i2tod(PG_FUNCTION_ARGS)
1294 : {
1295 164 : int16 num = PG_GETARG_INT16(0);
1296 :
1297 164 : PG_RETURN_FLOAT8((float8) num);
1298 : }
1299 :
1300 :
1301 : /*
1302 : * ftoi4 - converts a float4 number to an int4 number
1303 : */
1304 : Datum
1305 18 : ftoi4(PG_FUNCTION_ARGS)
1306 : {
1307 18 : float4 num = PG_GETARG_FLOAT4(0);
1308 :
1309 : /*
1310 : * Get rid of any fractional part in the input. This is so we don't fail
1311 : * on just-out-of-range values that would round into range. Note
1312 : * assumption that rint() will pass through a NaN or Inf unchanged.
1313 : */
1314 18 : num = rint(num);
1315 :
1316 : /* Range check */
1317 18 : if (unlikely(isnan(num) || !FLOAT4_FITS_IN_INT32(num)))
1318 8 : ereport(ERROR,
1319 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1320 : errmsg("integer out of range")));
1321 :
1322 10 : PG_RETURN_INT32((int32) num);
1323 : }
1324 :
1325 :
1326 : /*
1327 : * ftoi2 - converts a float4 number to an int2 number
1328 : */
1329 : Datum
1330 18 : ftoi2(PG_FUNCTION_ARGS)
1331 : {
1332 18 : float4 num = PG_GETARG_FLOAT4(0);
1333 :
1334 : /*
1335 : * Get rid of any fractional part in the input. This is so we don't fail
1336 : * on just-out-of-range values that would round into range. Note
1337 : * assumption that rint() will pass through a NaN or Inf unchanged.
1338 : */
1339 18 : num = rint(num);
1340 :
1341 : /* Range check */
1342 18 : if (unlikely(isnan(num) || !FLOAT4_FITS_IN_INT16(num)))
1343 8 : ereport(ERROR,
1344 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1345 : errmsg("smallint out of range")));
1346 :
1347 10 : PG_RETURN_INT16((int16) num);
1348 : }
1349 :
1350 :
1351 : /*
1352 : * i4tof - converts an int4 number to a float4 number
1353 : */
1354 : Datum
1355 367 : i4tof(PG_FUNCTION_ARGS)
1356 : {
1357 367 : int32 num = PG_GETARG_INT32(0);
1358 :
1359 367 : PG_RETURN_FLOAT4((float4) num);
1360 : }
1361 :
1362 :
1363 : /*
1364 : * i2tof - converts an int2 number to a float4 number
1365 : */
1366 : Datum
1367 0 : i2tof(PG_FUNCTION_ARGS)
1368 : {
1369 0 : int16 num = PG_GETARG_INT16(0);
1370 :
1371 0 : PG_RETURN_FLOAT4((float4) num);
1372 : }
1373 :
1374 :
1375 : /*
1376 : * =======================
1377 : * RANDOM FLOAT8 OPERATORS
1378 : * =======================
1379 : */
1380 :
1381 : /*
1382 : * dround - returns ROUND(arg1)
1383 : */
1384 : Datum
1385 592714 : dround(PG_FUNCTION_ARGS)
1386 : {
1387 592714 : float8 arg1 = PG_GETARG_FLOAT8(0);
1388 :
1389 592714 : PG_RETURN_FLOAT8(rint(arg1));
1390 : }
1391 :
1392 : /*
1393 : * dceil - returns the smallest integer greater than or
1394 : * equal to the specified float
1395 : */
1396 : Datum
1397 34240 : dceil(PG_FUNCTION_ARGS)
1398 : {
1399 34240 : float8 arg1 = PG_GETARG_FLOAT8(0);
1400 :
1401 34240 : PG_RETURN_FLOAT8(ceil(arg1));
1402 : }
1403 :
1404 : /*
1405 : * dfloor - returns the largest integer lesser than or
1406 : * equal to the specified float
1407 : */
1408 : Datum
1409 40 : dfloor(PG_FUNCTION_ARGS)
1410 : {
1411 40 : float8 arg1 = PG_GETARG_FLOAT8(0);
1412 :
1413 40 : PG_RETURN_FLOAT8(floor(arg1));
1414 : }
1415 :
1416 : /*
1417 : * dsign - returns -1 if the argument is less than 0, 0
1418 : * if the argument is equal to 0, and 1 if the
1419 : * argument is greater than zero.
1420 : */
1421 : Datum
1422 20 : dsign(PG_FUNCTION_ARGS)
1423 : {
1424 20 : float8 arg1 = PG_GETARG_FLOAT8(0);
1425 : float8 result;
1426 :
1427 20 : if (arg1 > 0)
1428 12 : result = 1.0;
1429 8 : else if (arg1 < 0)
1430 4 : result = -1.0;
1431 : else
1432 4 : result = 0.0;
1433 :
1434 20 : PG_RETURN_FLOAT8(result);
1435 : }
1436 :
1437 : /*
1438 : * dtrunc - returns truncation-towards-zero of arg1,
1439 : * arg1 >= 0 ... the greatest integer less
1440 : * than or equal to arg1
1441 : * arg1 < 0 ... the least integer greater
1442 : * than or equal to arg1
1443 : */
1444 : Datum
1445 23 : dtrunc(PG_FUNCTION_ARGS)
1446 : {
1447 23 : float8 arg1 = PG_GETARG_FLOAT8(0);
1448 : float8 result;
1449 :
1450 23 : if (arg1 >= 0)
1451 19 : result = floor(arg1);
1452 : else
1453 4 : result = -floor(-arg1);
1454 :
1455 23 : PG_RETURN_FLOAT8(result);
1456 : }
1457 :
1458 :
1459 : /*
1460 : * dsqrt - returns square root of arg1
1461 : */
1462 : Datum
1463 2052 : dsqrt(PG_FUNCTION_ARGS)
1464 : {
1465 2052 : float8 arg1 = PG_GETARG_FLOAT8(0);
1466 : float8 result;
1467 :
1468 2052 : if (arg1 < 0)
1469 0 : ereport(ERROR,
1470 : (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
1471 : errmsg("cannot take square root of a negative number")));
1472 :
1473 2052 : result = sqrt(arg1);
1474 2052 : if (unlikely(isinf(result)) && !isinf(arg1))
1475 0 : float_overflow_error();
1476 2052 : if (unlikely(result == 0.0) && arg1 != 0.0)
1477 0 : float_underflow_error();
1478 :
1479 2052 : PG_RETURN_FLOAT8(result);
1480 : }
1481 :
1482 :
1483 : /*
1484 : * dcbrt - returns cube root of arg1
1485 : */
1486 : Datum
1487 25 : dcbrt(PG_FUNCTION_ARGS)
1488 : {
1489 25 : float8 arg1 = PG_GETARG_FLOAT8(0);
1490 : float8 result;
1491 :
1492 25 : result = cbrt(arg1);
1493 25 : if (unlikely(isinf(result)) && !isinf(arg1))
1494 0 : float_overflow_error();
1495 25 : if (unlikely(result == 0.0) && arg1 != 0.0)
1496 0 : float_underflow_error();
1497 :
1498 25 : PG_RETURN_FLOAT8(result);
1499 : }
1500 :
1501 :
1502 : /*
1503 : * dpow - returns pow(arg1,arg2)
1504 : */
1505 : Datum
1506 484 : dpow(PG_FUNCTION_ARGS)
1507 : {
1508 484 : float8 arg1 = PG_GETARG_FLOAT8(0);
1509 484 : float8 arg2 = PG_GETARG_FLOAT8(1);
1510 : float8 result;
1511 :
1512 : /*
1513 : * The POSIX spec says that NaN ^ 0 = 1, and 1 ^ NaN = 1, while all other
1514 : * cases with NaN inputs yield NaN (with no error). Many older platforms
1515 : * get one or more of these cases wrong, so deal with them via explicit
1516 : * logic rather than trusting pow(3).
1517 : */
1518 484 : if (isnan(arg1))
1519 : {
1520 15 : if (isnan(arg2) || arg2 != 0.0)
1521 10 : PG_RETURN_FLOAT8(get_float8_nan());
1522 5 : PG_RETURN_FLOAT8(1.0);
1523 : }
1524 469 : if (isnan(arg2))
1525 : {
1526 15 : if (arg1 != 1.0)
1527 10 : PG_RETURN_FLOAT8(get_float8_nan());
1528 5 : PG_RETURN_FLOAT8(1.0);
1529 : }
1530 :
1531 : /*
1532 : * The SQL spec requires that we emit a particular SQLSTATE error code for
1533 : * certain error conditions. Specifically, we don't return a
1534 : * divide-by-zero error code for 0 ^ -1.
1535 : */
1536 454 : if (arg1 == 0 && arg2 < 0)
1537 4 : ereport(ERROR,
1538 : (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
1539 : errmsg("zero raised to a negative power is undefined")));
1540 450 : if (arg1 < 0 && floor(arg2) != arg2)
1541 4 : ereport(ERROR,
1542 : (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
1543 : errmsg("a negative number raised to a non-integer power yields a complex result")));
1544 :
1545 : /*
1546 : * We don't trust the platform's pow() to handle infinity cases per POSIX
1547 : * spec either, so deal with those explicitly too. It's easier to handle
1548 : * infinite y first, so that it doesn't matter if x is also infinite.
1549 : */
1550 446 : if (isinf(arg2))
1551 : {
1552 85 : float8 absx = fabs(arg1);
1553 :
1554 85 : if (absx == 1.0)
1555 20 : result = 1.0;
1556 65 : else if (arg2 > 0.0) /* y = +Inf */
1557 : {
1558 35 : if (absx > 1.0)
1559 20 : result = arg2;
1560 : else
1561 15 : result = 0.0;
1562 : }
1563 : else /* y = -Inf */
1564 : {
1565 30 : if (absx > 1.0)
1566 20 : result = 0.0;
1567 : else
1568 10 : result = -arg2;
1569 : }
1570 : }
1571 361 : else if (isinf(arg1))
1572 : {
1573 40 : if (arg2 == 0.0)
1574 10 : result = 1.0;
1575 30 : else if (arg1 > 0.0) /* x = +Inf */
1576 : {
1577 10 : if (arg2 > 0.0)
1578 5 : result = arg1;
1579 : else
1580 5 : result = 0.0;
1581 : }
1582 : else /* x = -Inf */
1583 : {
1584 : /*
1585 : * Per POSIX, the sign of the result depends on whether y is an
1586 : * odd integer. Since x < 0, we already know from the previous
1587 : * domain check that y is an integer. It is odd if y/2 is not
1588 : * also an integer.
1589 : */
1590 20 : float8 halfy = arg2 / 2; /* should be computed exactly */
1591 20 : bool yisoddinteger = (floor(halfy) != halfy);
1592 :
1593 20 : if (arg2 > 0.0)
1594 10 : result = yisoddinteger ? arg1 : -arg1;
1595 : else
1596 10 : result = yisoddinteger ? -0.0 : 0.0;
1597 : }
1598 : }
1599 : else
1600 : {
1601 : /*
1602 : * pow() sets errno on only some platforms, depending on whether it
1603 : * follows _IEEE_, _POSIX_, _XOPEN_, or _SVID_, so we must check both
1604 : * errno and invalid output values. (We can't rely on just the
1605 : * latter, either; some old platforms return a large-but-finite
1606 : * HUGE_VAL when reporting overflow.)
1607 : */
1608 321 : errno = 0;
1609 321 : result = pow(arg1, arg2);
1610 321 : if (errno == EDOM || isnan(result))
1611 : {
1612 : /*
1613 : * We handled all possible domain errors above, so this should be
1614 : * impossible. However, old glibc versions on x86 have a bug that
1615 : * causes them to fail this way for abs(y) greater than 2^63:
1616 : *
1617 : * https://sourceware.org/bugzilla/show_bug.cgi?id=3866
1618 : *
1619 : * Hence, if we get here, assume y is finite but large (large
1620 : * enough to be certainly even). The result should be 0 if x == 0,
1621 : * 1.0 if abs(x) == 1.0, otherwise an overflow or underflow error.
1622 : */
1623 0 : if (arg1 == 0.0)
1624 0 : result = 0.0; /* we already verified y is positive */
1625 : else
1626 : {
1627 0 : float8 absx = fabs(arg1);
1628 :
1629 0 : if (absx == 1.0)
1630 0 : result = 1.0;
1631 0 : else if (arg2 >= 0.0 ? (absx > 1.0) : (absx < 1.0))
1632 0 : float_overflow_error();
1633 : else
1634 0 : float_underflow_error();
1635 : }
1636 : }
1637 321 : else if (errno == ERANGE)
1638 : {
1639 4 : if (result != 0.0)
1640 4 : float_overflow_error();
1641 : else
1642 0 : float_underflow_error();
1643 : }
1644 : else
1645 : {
1646 317 : if (unlikely(isinf(result)))
1647 0 : float_overflow_error();
1648 317 : if (unlikely(result == 0.0) && arg1 != 0.0)
1649 0 : float_underflow_error();
1650 : }
1651 : }
1652 :
1653 442 : PG_RETURN_FLOAT8(result);
1654 : }
1655 :
1656 :
1657 : /*
1658 : * dexp - returns the exponential function of arg1
1659 : */
1660 : Datum
1661 39 : dexp(PG_FUNCTION_ARGS)
1662 : {
1663 39 : float8 arg1 = PG_GETARG_FLOAT8(0);
1664 : float8 result;
1665 :
1666 : /*
1667 : * Handle NaN and Inf cases explicitly. This avoids needing to assume
1668 : * that the platform's exp() conforms to POSIX for these cases, and it
1669 : * removes some edge cases for the overflow checks below.
1670 : */
1671 39 : if (isnan(arg1))
1672 5 : result = arg1;
1673 34 : else if (isinf(arg1))
1674 : {
1675 : /* Per POSIX, exp(-Inf) is 0 */
1676 10 : result = (arg1 > 0.0) ? arg1 : 0;
1677 : }
1678 : else
1679 : {
1680 : /*
1681 : * On some platforms, exp() will not set errno but just return Inf or
1682 : * zero to report overflow/underflow; therefore, test both cases.
1683 : */
1684 24 : errno = 0;
1685 24 : result = exp(arg1);
1686 24 : if (unlikely(errno == ERANGE))
1687 : {
1688 4 : if (result != 0.0)
1689 0 : float_overflow_error();
1690 : else
1691 4 : float_underflow_error();
1692 : }
1693 20 : else if (unlikely(isinf(result)))
1694 0 : float_overflow_error();
1695 20 : else if (unlikely(result == 0.0))
1696 0 : float_underflow_error();
1697 : }
1698 :
1699 35 : PG_RETURN_FLOAT8(result);
1700 : }
1701 :
1702 :
1703 : /*
1704 : * dlog1 - returns the natural logarithm of arg1
1705 : */
1706 : Datum
1707 20 : dlog1(PG_FUNCTION_ARGS)
1708 : {
1709 20 : float8 arg1 = PG_GETARG_FLOAT8(0);
1710 : float8 result;
1711 :
1712 : /*
1713 : * Emit particular SQLSTATE error codes for ln(). This is required by the
1714 : * SQL standard.
1715 : */
1716 20 : if (arg1 == 0.0)
1717 4 : ereport(ERROR,
1718 : (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
1719 : errmsg("cannot take logarithm of zero")));
1720 16 : if (arg1 < 0)
1721 4 : ereport(ERROR,
1722 : (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
1723 : errmsg("cannot take logarithm of a negative number")));
1724 :
1725 12 : result = log(arg1);
1726 12 : if (unlikely(isinf(result)) && !isinf(arg1))
1727 0 : float_overflow_error();
1728 12 : if (unlikely(result == 0.0) && arg1 != 1.0)
1729 0 : float_underflow_error();
1730 :
1731 12 : PG_RETURN_FLOAT8(result);
1732 : }
1733 :
1734 :
1735 : /*
1736 : * dlog10 - returns the base 10 logarithm of arg1
1737 : */
1738 : Datum
1739 0 : dlog10(PG_FUNCTION_ARGS)
1740 : {
1741 0 : float8 arg1 = PG_GETARG_FLOAT8(0);
1742 : float8 result;
1743 :
1744 : /*
1745 : * Emit particular SQLSTATE error codes for log(). The SQL spec doesn't
1746 : * define log(), but it does define ln(), so it makes sense to emit the
1747 : * same error code for an analogous error condition.
1748 : */
1749 0 : if (arg1 == 0.0)
1750 0 : ereport(ERROR,
1751 : (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
1752 : errmsg("cannot take logarithm of zero")));
1753 0 : if (arg1 < 0)
1754 0 : ereport(ERROR,
1755 : (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
1756 : errmsg("cannot take logarithm of a negative number")));
1757 :
1758 0 : result = log10(arg1);
1759 0 : if (unlikely(isinf(result)) && !isinf(arg1))
1760 0 : float_overflow_error();
1761 0 : if (unlikely(result == 0.0) && arg1 != 1.0)
1762 0 : float_underflow_error();
1763 :
1764 0 : PG_RETURN_FLOAT8(result);
1765 : }
1766 :
1767 :
1768 : /*
1769 : * dacos - returns the arccos of arg1 (radians)
1770 : */
1771 : Datum
1772 0 : dacos(PG_FUNCTION_ARGS)
1773 : {
1774 0 : float8 arg1 = PG_GETARG_FLOAT8(0);
1775 : float8 result;
1776 :
1777 : /* Per the POSIX spec, return NaN if the input is NaN */
1778 0 : if (isnan(arg1))
1779 0 : PG_RETURN_FLOAT8(get_float8_nan());
1780 :
1781 : /*
1782 : * The principal branch of the inverse cosine function maps values in the
1783 : * range [-1, 1] to values in the range [0, Pi], so we should reject any
1784 : * inputs outside that range and the result will always be finite.
1785 : */
1786 0 : if (arg1 < -1.0 || arg1 > 1.0)
1787 0 : ereport(ERROR,
1788 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1789 : errmsg("input is out of range")));
1790 :
1791 0 : result = acos(arg1);
1792 0 : if (unlikely(isinf(result)))
1793 0 : float_overflow_error();
1794 :
1795 0 : PG_RETURN_FLOAT8(result);
1796 : }
1797 :
1798 :
1799 : /*
1800 : * dasin - returns the arcsin of arg1 (radians)
1801 : */
1802 : Datum
1803 55 : dasin(PG_FUNCTION_ARGS)
1804 : {
1805 55 : float8 arg1 = PG_GETARG_FLOAT8(0);
1806 : float8 result;
1807 :
1808 : /* Per the POSIX spec, return NaN if the input is NaN */
1809 55 : if (isnan(arg1))
1810 0 : PG_RETURN_FLOAT8(get_float8_nan());
1811 :
1812 : /*
1813 : * The principal branch of the inverse sine function maps values in the
1814 : * range [-1, 1] to values in the range [-Pi/2, Pi/2], so we should reject
1815 : * any inputs outside that range and the result will always be finite.
1816 : */
1817 55 : if (arg1 < -1.0 || arg1 > 1.0)
1818 0 : ereport(ERROR,
1819 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1820 : errmsg("input is out of range")));
1821 :
1822 55 : result = asin(arg1);
1823 55 : if (unlikely(isinf(result)))
1824 0 : float_overflow_error();
1825 :
1826 55 : PG_RETURN_FLOAT8(result);
1827 : }
1828 :
1829 :
1830 : /*
1831 : * datan - returns the arctan of arg1 (radians)
1832 : */
1833 : Datum
1834 0 : datan(PG_FUNCTION_ARGS)
1835 : {
1836 0 : float8 arg1 = PG_GETARG_FLOAT8(0);
1837 : float8 result;
1838 :
1839 : /* Per the POSIX spec, return NaN if the input is NaN */
1840 0 : if (isnan(arg1))
1841 0 : PG_RETURN_FLOAT8(get_float8_nan());
1842 :
1843 : /*
1844 : * The principal branch of the inverse tangent function maps all inputs to
1845 : * values in the range [-Pi/2, Pi/2], so the result should always be
1846 : * finite, even if the input is infinite.
1847 : */
1848 0 : result = atan(arg1);
1849 0 : if (unlikely(isinf(result)))
1850 0 : float_overflow_error();
1851 :
1852 0 : PG_RETURN_FLOAT8(result);
1853 : }
1854 :
1855 :
1856 : /*
1857 : * atan2 - returns the arctan of arg1/arg2 (radians)
1858 : */
1859 : Datum
1860 20 : datan2(PG_FUNCTION_ARGS)
1861 : {
1862 20 : float8 arg1 = PG_GETARG_FLOAT8(0);
1863 20 : float8 arg2 = PG_GETARG_FLOAT8(1);
1864 : float8 result;
1865 :
1866 : /* Per the POSIX spec, return NaN if either input is NaN */
1867 20 : if (isnan(arg1) || isnan(arg2))
1868 0 : PG_RETURN_FLOAT8(get_float8_nan());
1869 :
1870 : /*
1871 : * atan2 maps all inputs to values in the range [-Pi, Pi], so the result
1872 : * should always be finite, even if the inputs are infinite.
1873 : */
1874 20 : result = atan2(arg1, arg2);
1875 20 : if (unlikely(isinf(result)))
1876 0 : float_overflow_error();
1877 :
1878 20 : PG_RETURN_FLOAT8(result);
1879 : }
1880 :
1881 :
1882 : /*
1883 : * dcos - returns the cosine of arg1 (radians)
1884 : */
1885 : Datum
1886 589 : dcos(PG_FUNCTION_ARGS)
1887 : {
1888 589 : float8 arg1 = PG_GETARG_FLOAT8(0);
1889 : float8 result;
1890 :
1891 : /* Per the POSIX spec, return NaN if the input is NaN */
1892 589 : if (isnan(arg1))
1893 0 : PG_RETURN_FLOAT8(get_float8_nan());
1894 :
1895 : /*
1896 : * cos() is periodic and so theoretically can work for all finite inputs,
1897 : * but some implementations may choose to throw error if the input is so
1898 : * large that there are no significant digits in the result. So we should
1899 : * check for errors. POSIX allows an error to be reported either via
1900 : * errno or via fetestexcept(), but currently we only support checking
1901 : * errno. (fetestexcept() is rumored to report underflow unreasonably
1902 : * early on some platforms, so it's not clear that believing it would be a
1903 : * net improvement anyway.)
1904 : *
1905 : * For infinite inputs, POSIX specifies that the trigonometric functions
1906 : * should return a domain error; but we won't notice that unless the
1907 : * platform reports via errno, so also explicitly test for infinite
1908 : * inputs.
1909 : */
1910 589 : errno = 0;
1911 589 : result = cos(arg1);
1912 589 : if (errno != 0 || isinf(arg1))
1913 0 : ereport(ERROR,
1914 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1915 : errmsg("input is out of range")));
1916 589 : if (unlikely(isinf(result)))
1917 0 : float_overflow_error();
1918 :
1919 589 : PG_RETURN_FLOAT8(result);
1920 : }
1921 :
1922 :
1923 : /*
1924 : * dcot - returns the cotangent of arg1 (radians)
1925 : */
1926 : Datum
1927 0 : dcot(PG_FUNCTION_ARGS)
1928 : {
1929 0 : float8 arg1 = PG_GETARG_FLOAT8(0);
1930 : float8 result;
1931 :
1932 : /* Per the POSIX spec, return NaN if the input is NaN */
1933 0 : if (isnan(arg1))
1934 0 : PG_RETURN_FLOAT8(get_float8_nan());
1935 :
1936 : /* Be sure to throw an error if the input is infinite --- see dcos() */
1937 0 : errno = 0;
1938 0 : result = tan(arg1);
1939 0 : if (errno != 0 || isinf(arg1))
1940 0 : ereport(ERROR,
1941 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1942 : errmsg("input is out of range")));
1943 :
1944 0 : result = 1.0 / result;
1945 : /* Not checking for overflow because cot(0) == Inf */
1946 :
1947 0 : PG_RETURN_FLOAT8(result);
1948 : }
1949 :
1950 :
1951 : /*
1952 : * dsin - returns the sine of arg1 (radians)
1953 : */
1954 : Datum
1955 518 : dsin(PG_FUNCTION_ARGS)
1956 : {
1957 518 : float8 arg1 = PG_GETARG_FLOAT8(0);
1958 : float8 result;
1959 :
1960 : /* Per the POSIX spec, return NaN if the input is NaN */
1961 518 : if (isnan(arg1))
1962 0 : PG_RETURN_FLOAT8(get_float8_nan());
1963 :
1964 : /* Be sure to throw an error if the input is infinite --- see dcos() */
1965 518 : errno = 0;
1966 518 : result = sin(arg1);
1967 518 : if (errno != 0 || isinf(arg1))
1968 0 : ereport(ERROR,
1969 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1970 : errmsg("input is out of range")));
1971 518 : if (unlikely(isinf(result)))
1972 0 : float_overflow_error();
1973 :
1974 518 : PG_RETURN_FLOAT8(result);
1975 : }
1976 :
1977 :
1978 : /*
1979 : * dtan - returns the tangent of arg1 (radians)
1980 : */
1981 : Datum
1982 0 : dtan(PG_FUNCTION_ARGS)
1983 : {
1984 0 : float8 arg1 = PG_GETARG_FLOAT8(0);
1985 : float8 result;
1986 :
1987 : /* Per the POSIX spec, return NaN if the input is NaN */
1988 0 : if (isnan(arg1))
1989 0 : PG_RETURN_FLOAT8(get_float8_nan());
1990 :
1991 : /* Be sure to throw an error if the input is infinite --- see dcos() */
1992 0 : errno = 0;
1993 0 : result = tan(arg1);
1994 0 : if (errno != 0 || isinf(arg1))
1995 0 : ereport(ERROR,
1996 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1997 : errmsg("input is out of range")));
1998 : /* Not checking for overflow because tan(pi/2) == Inf */
1999 :
2000 0 : PG_RETURN_FLOAT8(result);
2001 : }
2002 :
2003 :
2004 : /* ========== DEGREE-BASED TRIGONOMETRIC FUNCTIONS ========== */
2005 :
2006 :
2007 : /*
2008 : * Initialize the cached constants declared at the head of this file
2009 : * (sin_30 etc). The fact that we need those at all, let alone need this
2010 : * Rube-Goldberg-worthy method of initializing them, is because there are
2011 : * compilers out there that will precompute expressions such as sin(constant)
2012 : * using a sin() function different from what will be used at runtime. If we
2013 : * want exact results, we must ensure that none of the scaling constants used
2014 : * in the degree-based trig functions are computed that way. To do so, we
2015 : * compute them from the variables degree_c_thirty etc, which are also really
2016 : * constants, but the compiler cannot assume that.
2017 : *
2018 : * Other hazards we are trying to forestall with this kluge include the
2019 : * possibility that compilers will rearrange the expressions, or compute
2020 : * some intermediate results in registers wider than a standard double.
2021 : *
2022 : * In the places where we use these constants, the typical pattern is like
2023 : * volatile float8 sin_x = sin(x * RADIANS_PER_DEGREE);
2024 : * return (sin_x / sin_30);
2025 : * where we hope to get a value of exactly 1.0 from the division when x = 30.
2026 : * The volatile temporary variable is needed on machines with wide float
2027 : * registers, to ensure that the result of sin(x) is rounded to double width
2028 : * the same as the value of sin_30 has been. Experimentation with gcc shows
2029 : * that marking the temp variable volatile is necessary to make the store and
2030 : * reload actually happen; hopefully the same trick works for other compilers.
2031 : * (gcc's documentation suggests using the -ffloat-store compiler switch to
2032 : * ensure this, but that is compiler-specific and it also pessimizes code in
2033 : * many places where we don't care about this.)
2034 : */
2035 : static void
2036 4 : init_degree_constants(void)
2037 : {
2038 4 : sin_30 = sin(degree_c_thirty * RADIANS_PER_DEGREE);
2039 4 : one_minus_cos_60 = 1.0 - cos(degree_c_sixty * RADIANS_PER_DEGREE);
2040 4 : asin_0_5 = asin(degree_c_one_half);
2041 4 : acos_0_5 = acos(degree_c_one_half);
2042 4 : atan_1_0 = atan(degree_c_one);
2043 4 : tan_45 = sind_q1(degree_c_forty_five) / cosd_q1(degree_c_forty_five);
2044 4 : cot_45 = cosd_q1(degree_c_forty_five) / sind_q1(degree_c_forty_five);
2045 4 : degree_consts_set = true;
2046 4 : }
2047 :
2048 : #define INIT_DEGREE_CONSTANTS() \
2049 : do { \
2050 : if (!degree_consts_set) \
2051 : init_degree_constants(); \
2052 : } while(0)
2053 :
2054 :
2055 : /*
2056 : * asind_q1 - returns the inverse sine of x in degrees, for x in
2057 : * the range [0, 1]. The result is an angle in the
2058 : * first quadrant --- [0, 90] degrees.
2059 : *
2060 : * For the 3 special case inputs (0, 0.5 and 1), this
2061 : * function will return exact values (0, 30 and 90
2062 : * degrees respectively).
2063 : */
2064 : static double
2065 56 : asind_q1(double x)
2066 : {
2067 : /*
2068 : * Stitch together inverse sine and cosine functions for the ranges [0,
2069 : * 0.5] and (0.5, 1]. Each expression below is guaranteed to return
2070 : * exactly 30 for x=0.5, so the result is a continuous monotonic function
2071 : * over the full range.
2072 : */
2073 56 : if (x <= 0.5)
2074 : {
2075 32 : volatile float8 asin_x = asin(x);
2076 :
2077 32 : return (asin_x / asin_0_5) * 30.0;
2078 : }
2079 : else
2080 : {
2081 24 : volatile float8 acos_x = acos(x);
2082 :
2083 24 : return 90.0 - (acos_x / acos_0_5) * 60.0;
2084 : }
2085 : }
2086 :
2087 :
2088 : /*
2089 : * acosd_q1 - returns the inverse cosine of x in degrees, for x in
2090 : * the range [0, 1]. The result is an angle in the
2091 : * first quadrant --- [0, 90] degrees.
2092 : *
2093 : * For the 3 special case inputs (0, 0.5 and 1), this
2094 : * function will return exact values (0, 60 and 90
2095 : * degrees respectively).
2096 : */
2097 : static double
2098 24 : acosd_q1(double x)
2099 : {
2100 : /*
2101 : * Stitch together inverse sine and cosine functions for the ranges [0,
2102 : * 0.5] and (0.5, 1]. Each expression below is guaranteed to return
2103 : * exactly 60 for x=0.5, so the result is a continuous monotonic function
2104 : * over the full range.
2105 : */
2106 24 : if (x <= 0.5)
2107 : {
2108 16 : volatile float8 asin_x = asin(x);
2109 :
2110 16 : return 90.0 - (asin_x / asin_0_5) * 30.0;
2111 : }
2112 : else
2113 : {
2114 8 : volatile float8 acos_x = acos(x);
2115 :
2116 8 : return (acos_x / acos_0_5) * 60.0;
2117 : }
2118 : }
2119 :
2120 :
2121 : /*
2122 : * dacosd - returns the arccos of arg1 (degrees)
2123 : */
2124 : Datum
2125 40 : dacosd(PG_FUNCTION_ARGS)
2126 : {
2127 40 : float8 arg1 = PG_GETARG_FLOAT8(0);
2128 : float8 result;
2129 :
2130 : /* Per the POSIX spec, return NaN if the input is NaN */
2131 40 : if (isnan(arg1))
2132 0 : PG_RETURN_FLOAT8(get_float8_nan());
2133 :
2134 40 : INIT_DEGREE_CONSTANTS();
2135 :
2136 : /*
2137 : * The principal branch of the inverse cosine function maps values in the
2138 : * range [-1, 1] to values in the range [0, 180], so we should reject any
2139 : * inputs outside that range and the result will always be finite.
2140 : */
2141 40 : if (arg1 < -1.0 || arg1 > 1.0)
2142 0 : ereport(ERROR,
2143 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
2144 : errmsg("input is out of range")));
2145 :
2146 40 : if (arg1 >= 0.0)
2147 24 : result = acosd_q1(arg1);
2148 : else
2149 16 : result = 90.0 + asind_q1(-arg1);
2150 :
2151 40 : if (unlikely(isinf(result)))
2152 0 : float_overflow_error();
2153 :
2154 40 : PG_RETURN_FLOAT8(result);
2155 : }
2156 :
2157 :
2158 : /*
2159 : * dasind - returns the arcsin of arg1 (degrees)
2160 : */
2161 : Datum
2162 40 : dasind(PG_FUNCTION_ARGS)
2163 : {
2164 40 : float8 arg1 = PG_GETARG_FLOAT8(0);
2165 : float8 result;
2166 :
2167 : /* Per the POSIX spec, return NaN if the input is NaN */
2168 40 : if (isnan(arg1))
2169 0 : PG_RETURN_FLOAT8(get_float8_nan());
2170 :
2171 40 : INIT_DEGREE_CONSTANTS();
2172 :
2173 : /*
2174 : * The principal branch of the inverse sine function maps values in the
2175 : * range [-1, 1] to values in the range [-90, 90], so we should reject any
2176 : * inputs outside that range and the result will always be finite.
2177 : */
2178 40 : if (arg1 < -1.0 || arg1 > 1.0)
2179 0 : ereport(ERROR,
2180 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
2181 : errmsg("input is out of range")));
2182 :
2183 40 : if (arg1 >= 0.0)
2184 24 : result = asind_q1(arg1);
2185 : else
2186 16 : result = -asind_q1(-arg1);
2187 :
2188 40 : if (unlikely(isinf(result)))
2189 0 : float_overflow_error();
2190 :
2191 40 : PG_RETURN_FLOAT8(result);
2192 : }
2193 :
2194 :
2195 : /*
2196 : * datand - returns the arctan of arg1 (degrees)
2197 : */
2198 : Datum
2199 40 : datand(PG_FUNCTION_ARGS)
2200 : {
2201 40 : float8 arg1 = PG_GETARG_FLOAT8(0);
2202 : float8 result;
2203 : volatile float8 atan_arg1;
2204 :
2205 : /* Per the POSIX spec, return NaN if the input is NaN */
2206 40 : if (isnan(arg1))
2207 0 : PG_RETURN_FLOAT8(get_float8_nan());
2208 :
2209 40 : INIT_DEGREE_CONSTANTS();
2210 :
2211 : /*
2212 : * The principal branch of the inverse tangent function maps all inputs to
2213 : * values in the range [-90, 90], so the result should always be finite,
2214 : * even if the input is infinite. Additionally, we take care to ensure
2215 : * than when arg1 is 1, the result is exactly 45.
2216 : */
2217 40 : atan_arg1 = atan(arg1);
2218 40 : result = (atan_arg1 / atan_1_0) * 45.0;
2219 :
2220 40 : if (unlikely(isinf(result)))
2221 0 : float_overflow_error();
2222 :
2223 40 : PG_RETURN_FLOAT8(result);
2224 : }
2225 :
2226 :
2227 : /*
2228 : * atan2d - returns the arctan of arg1/arg2 (degrees)
2229 : */
2230 : Datum
2231 40 : datan2d(PG_FUNCTION_ARGS)
2232 : {
2233 40 : float8 arg1 = PG_GETARG_FLOAT8(0);
2234 40 : float8 arg2 = PG_GETARG_FLOAT8(1);
2235 : float8 result;
2236 : volatile float8 atan2_arg1_arg2;
2237 :
2238 : /* Per the POSIX spec, return NaN if either input is NaN */
2239 40 : if (isnan(arg1) || isnan(arg2))
2240 0 : PG_RETURN_FLOAT8(get_float8_nan());
2241 :
2242 40 : INIT_DEGREE_CONSTANTS();
2243 :
2244 : /*
2245 : * atan2d maps all inputs to values in the range [-180, 180], so the
2246 : * result should always be finite, even if the inputs are infinite.
2247 : *
2248 : * Note: this coding assumes that atan(1.0) is a suitable scaling constant
2249 : * to get an exact result from atan2(). This might well fail on us at
2250 : * some point, requiring us to decide exactly what inputs we think we're
2251 : * going to guarantee an exact result for.
2252 : */
2253 40 : atan2_arg1_arg2 = atan2(arg1, arg2);
2254 40 : result = (atan2_arg1_arg2 / atan_1_0) * 45.0;
2255 :
2256 40 : if (unlikely(isinf(result)))
2257 0 : float_overflow_error();
2258 :
2259 40 : PG_RETURN_FLOAT8(result);
2260 : }
2261 :
2262 :
2263 : /*
2264 : * sind_0_to_30 - returns the sine of an angle that lies between 0 and
2265 : * 30 degrees. This will return exactly 0 when x is 0,
2266 : * and exactly 0.5 when x is 30 degrees.
2267 : */
2268 : static double
2269 212 : sind_0_to_30(double x)
2270 : {
2271 212 : volatile float8 sin_x = sin(x * RADIANS_PER_DEGREE);
2272 :
2273 212 : return (sin_x / sin_30) / 2.0;
2274 : }
2275 :
2276 :
2277 : /*
2278 : * cosd_0_to_60 - returns the cosine of an angle that lies between 0
2279 : * and 60 degrees. This will return exactly 1 when x
2280 : * is 0, and exactly 0.5 when x is 60 degrees.
2281 : */
2282 : static double
2283 356 : cosd_0_to_60(double x)
2284 : {
2285 356 : volatile float8 one_minus_cos_x = 1.0 - cos(x * RADIANS_PER_DEGREE);
2286 :
2287 356 : return 1.0 - (one_minus_cos_x / one_minus_cos_60) / 2.0;
2288 : }
2289 :
2290 :
2291 : /*
2292 : * sind_q1 - returns the sine of an angle in the first quadrant
2293 : * (0 to 90 degrees).
2294 : */
2295 : static double
2296 284 : sind_q1(double x)
2297 : {
2298 : /*
2299 : * Stitch together the sine and cosine functions for the ranges [0, 30]
2300 : * and (30, 90]. These guarantee to return exact answers at their
2301 : * endpoints, so the overall result is a continuous monotonic function
2302 : * that gives exact results when x = 0, 30 and 90 degrees.
2303 : */
2304 284 : if (x <= 30.0)
2305 140 : return sind_0_to_30(x);
2306 : else
2307 144 : return cosd_0_to_60(90.0 - x);
2308 : }
2309 :
2310 :
2311 : /*
2312 : * cosd_q1 - returns the cosine of an angle in the first quadrant
2313 : * (0 to 90 degrees).
2314 : */
2315 : static double
2316 284 : cosd_q1(double x)
2317 : {
2318 : /*
2319 : * Stitch together the sine and cosine functions for the ranges [0, 60]
2320 : * and (60, 90]. These guarantee to return exact answers at their
2321 : * endpoints, so the overall result is a continuous monotonic function
2322 : * that gives exact results when x = 0, 60 and 90 degrees.
2323 : */
2324 284 : if (x <= 60.0)
2325 212 : return cosd_0_to_60(x);
2326 : else
2327 72 : return sind_0_to_30(90.0 - x);
2328 : }
2329 :
2330 :
2331 : /*
2332 : * dcosd - returns the cosine of arg1 (degrees)
2333 : */
2334 : Datum
2335 132 : dcosd(PG_FUNCTION_ARGS)
2336 : {
2337 132 : float8 arg1 = PG_GETARG_FLOAT8(0);
2338 : float8 result;
2339 132 : int sign = 1;
2340 :
2341 : /*
2342 : * Per the POSIX spec, return NaN if the input is NaN and throw an error
2343 : * if the input is infinite.
2344 : */
2345 132 : if (isnan(arg1))
2346 0 : PG_RETURN_FLOAT8(get_float8_nan());
2347 :
2348 132 : if (isinf(arg1))
2349 0 : ereport(ERROR,
2350 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
2351 : errmsg("input is out of range")));
2352 :
2353 132 : INIT_DEGREE_CONSTANTS();
2354 :
2355 : /* Reduce the range of the input to [0,90] degrees */
2356 132 : arg1 = fmod(arg1, 360.0);
2357 :
2358 132 : if (arg1 < 0.0)
2359 : {
2360 : /* cosd(-x) = cosd(x) */
2361 0 : arg1 = -arg1;
2362 : }
2363 :
2364 132 : if (arg1 > 180.0)
2365 : {
2366 : /* cosd(360-x) = cosd(x) */
2367 36 : arg1 = 360.0 - arg1;
2368 : }
2369 :
2370 132 : if (arg1 > 90.0)
2371 : {
2372 : /* cosd(180-x) = -cosd(x) */
2373 36 : arg1 = 180.0 - arg1;
2374 36 : sign = -sign;
2375 : }
2376 :
2377 132 : result = sign * cosd_q1(arg1);
2378 :
2379 132 : if (unlikely(isinf(result)))
2380 0 : float_overflow_error();
2381 :
2382 132 : PG_RETURN_FLOAT8(result);
2383 : }
2384 :
2385 :
2386 : /*
2387 : * dcotd - returns the cotangent of arg1 (degrees)
2388 : */
2389 : Datum
2390 72 : dcotd(PG_FUNCTION_ARGS)
2391 : {
2392 72 : float8 arg1 = PG_GETARG_FLOAT8(0);
2393 : float8 result;
2394 : volatile float8 cot_arg1;
2395 72 : int sign = 1;
2396 :
2397 : /*
2398 : * Per the POSIX spec, return NaN if the input is NaN and throw an error
2399 : * if the input is infinite.
2400 : */
2401 72 : if (isnan(arg1))
2402 0 : PG_RETURN_FLOAT8(get_float8_nan());
2403 :
2404 72 : if (isinf(arg1))
2405 0 : ereport(ERROR,
2406 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
2407 : errmsg("input is out of range")));
2408 :
2409 72 : INIT_DEGREE_CONSTANTS();
2410 :
2411 : /* Reduce the range of the input to [0,90] degrees */
2412 72 : arg1 = fmod(arg1, 360.0);
2413 :
2414 72 : if (arg1 < 0.0)
2415 : {
2416 : /* cotd(-x) = -cotd(x) */
2417 0 : arg1 = -arg1;
2418 0 : sign = -sign;
2419 : }
2420 :
2421 72 : if (arg1 > 180.0)
2422 : {
2423 : /* cotd(360-x) = -cotd(x) */
2424 24 : arg1 = 360.0 - arg1;
2425 24 : sign = -sign;
2426 : }
2427 :
2428 72 : if (arg1 > 90.0)
2429 : {
2430 : /* cotd(180-x) = -cotd(x) */
2431 24 : arg1 = 180.0 - arg1;
2432 24 : sign = -sign;
2433 : }
2434 :
2435 72 : cot_arg1 = cosd_q1(arg1) / sind_q1(arg1);
2436 72 : result = sign * (cot_arg1 / cot_45);
2437 :
2438 : /*
2439 : * On some machines we get cotd(270) = minus zero, but this isn't always
2440 : * true. For portability, and because the user constituency for this
2441 : * function probably doesn't want minus zero, force it to plain zero.
2442 : */
2443 72 : if (result == 0.0)
2444 16 : result = 0.0;
2445 :
2446 : /* Not checking for overflow because cotd(0) == Inf */
2447 :
2448 72 : PG_RETURN_FLOAT8(result);
2449 : }
2450 :
2451 :
2452 : /*
2453 : * dsind - returns the sine of arg1 (degrees)
2454 : */
2455 : Datum
2456 132 : dsind(PG_FUNCTION_ARGS)
2457 : {
2458 132 : float8 arg1 = PG_GETARG_FLOAT8(0);
2459 : float8 result;
2460 132 : int sign = 1;
2461 :
2462 : /*
2463 : * Per the POSIX spec, return NaN if the input is NaN and throw an error
2464 : * if the input is infinite.
2465 : */
2466 132 : if (isnan(arg1))
2467 0 : PG_RETURN_FLOAT8(get_float8_nan());
2468 :
2469 132 : if (isinf(arg1))
2470 0 : ereport(ERROR,
2471 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
2472 : errmsg("input is out of range")));
2473 :
2474 132 : INIT_DEGREE_CONSTANTS();
2475 :
2476 : /* Reduce the range of the input to [0,90] degrees */
2477 132 : arg1 = fmod(arg1, 360.0);
2478 :
2479 132 : if (arg1 < 0.0)
2480 : {
2481 : /* sind(-x) = -sind(x) */
2482 0 : arg1 = -arg1;
2483 0 : sign = -sign;
2484 : }
2485 :
2486 132 : if (arg1 > 180.0)
2487 : {
2488 : /* sind(360-x) = -sind(x) */
2489 36 : arg1 = 360.0 - arg1;
2490 36 : sign = -sign;
2491 : }
2492 :
2493 132 : if (arg1 > 90.0)
2494 : {
2495 : /* sind(180-x) = sind(x) */
2496 36 : arg1 = 180.0 - arg1;
2497 : }
2498 :
2499 132 : result = sign * sind_q1(arg1);
2500 :
2501 132 : if (unlikely(isinf(result)))
2502 0 : float_overflow_error();
2503 :
2504 132 : PG_RETURN_FLOAT8(result);
2505 : }
2506 :
2507 :
2508 : /*
2509 : * dtand - returns the tangent of arg1 (degrees)
2510 : */
2511 : Datum
2512 72 : dtand(PG_FUNCTION_ARGS)
2513 : {
2514 72 : float8 arg1 = PG_GETARG_FLOAT8(0);
2515 : float8 result;
2516 : volatile float8 tan_arg1;
2517 72 : int sign = 1;
2518 :
2519 : /*
2520 : * Per the POSIX spec, return NaN if the input is NaN and throw an error
2521 : * if the input is infinite.
2522 : */
2523 72 : if (isnan(arg1))
2524 0 : PG_RETURN_FLOAT8(get_float8_nan());
2525 :
2526 72 : if (isinf(arg1))
2527 0 : ereport(ERROR,
2528 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
2529 : errmsg("input is out of range")));
2530 :
2531 72 : INIT_DEGREE_CONSTANTS();
2532 :
2533 : /* Reduce the range of the input to [0,90] degrees */
2534 72 : arg1 = fmod(arg1, 360.0);
2535 :
2536 72 : if (arg1 < 0.0)
2537 : {
2538 : /* tand(-x) = -tand(x) */
2539 0 : arg1 = -arg1;
2540 0 : sign = -sign;
2541 : }
2542 :
2543 72 : if (arg1 > 180.0)
2544 : {
2545 : /* tand(360-x) = -tand(x) */
2546 24 : arg1 = 360.0 - arg1;
2547 24 : sign = -sign;
2548 : }
2549 :
2550 72 : if (arg1 > 90.0)
2551 : {
2552 : /* tand(180-x) = -tand(x) */
2553 24 : arg1 = 180.0 - arg1;
2554 24 : sign = -sign;
2555 : }
2556 :
2557 72 : tan_arg1 = sind_q1(arg1) / cosd_q1(arg1);
2558 72 : result = sign * (tan_arg1 / tan_45);
2559 :
2560 : /*
2561 : * On some machines we get tand(180) = minus zero, but this isn't always
2562 : * true. For portability, and because the user constituency for this
2563 : * function probably doesn't want minus zero, force it to plain zero.
2564 : */
2565 72 : if (result == 0.0)
2566 24 : result = 0.0;
2567 :
2568 : /* Not checking for overflow because tand(90) == Inf */
2569 :
2570 72 : PG_RETURN_FLOAT8(result);
2571 : }
2572 :
2573 :
2574 : /*
2575 : * degrees - returns degrees converted from radians
2576 : */
2577 : Datum
2578 40 : degrees(PG_FUNCTION_ARGS)
2579 : {
2580 40 : float8 arg1 = PG_GETARG_FLOAT8(0);
2581 :
2582 40 : PG_RETURN_FLOAT8(float8_div(arg1, RADIANS_PER_DEGREE));
2583 : }
2584 :
2585 :
2586 : /*
2587 : * dpi - returns the constant PI
2588 : */
2589 : Datum
2590 36 : dpi(PG_FUNCTION_ARGS)
2591 : {
2592 36 : PG_RETURN_FLOAT8(M_PI);
2593 : }
2594 :
2595 :
2596 : /*
2597 : * radians - returns radians converted from degrees
2598 : */
2599 : Datum
2600 955 : radians(PG_FUNCTION_ARGS)
2601 : {
2602 955 : float8 arg1 = PG_GETARG_FLOAT8(0);
2603 :
2604 955 : PG_RETURN_FLOAT8(float8_mul(arg1, RADIANS_PER_DEGREE));
2605 : }
2606 :
2607 :
2608 : /* ========== HYPERBOLIC FUNCTIONS ========== */
2609 :
2610 :
2611 : /*
2612 : * dsinh - returns the hyperbolic sine of arg1
2613 : */
2614 : Datum
2615 20 : dsinh(PG_FUNCTION_ARGS)
2616 : {
2617 20 : float8 arg1 = PG_GETARG_FLOAT8(0);
2618 : float8 result;
2619 :
2620 20 : errno = 0;
2621 20 : result = sinh(arg1);
2622 :
2623 : /*
2624 : * if an ERANGE error occurs, it means there is an overflow. For sinh,
2625 : * the result should be either -infinity or infinity, depending on the
2626 : * sign of arg1.
2627 : */
2628 20 : if (errno == ERANGE)
2629 : {
2630 0 : if (arg1 < 0)
2631 0 : result = -get_float8_infinity();
2632 : else
2633 0 : result = get_float8_infinity();
2634 : }
2635 :
2636 20 : PG_RETURN_FLOAT8(result);
2637 : }
2638 :
2639 :
2640 : /*
2641 : * dcosh - returns the hyperbolic cosine of arg1
2642 : */
2643 : Datum
2644 20 : dcosh(PG_FUNCTION_ARGS)
2645 : {
2646 20 : float8 arg1 = PG_GETARG_FLOAT8(0);
2647 : float8 result;
2648 :
2649 20 : errno = 0;
2650 20 : result = cosh(arg1);
2651 :
2652 : /*
2653 : * if an ERANGE error occurs, it means there is an overflow. As cosh is
2654 : * always positive, it always means the result is positive infinity.
2655 : */
2656 20 : if (errno == ERANGE)
2657 0 : result = get_float8_infinity();
2658 :
2659 20 : if (unlikely(result == 0.0))
2660 0 : float_underflow_error();
2661 :
2662 20 : PG_RETURN_FLOAT8(result);
2663 : }
2664 :
2665 : /*
2666 : * dtanh - returns the hyperbolic tangent of arg1
2667 : */
2668 : Datum
2669 20 : dtanh(PG_FUNCTION_ARGS)
2670 : {
2671 20 : float8 arg1 = PG_GETARG_FLOAT8(0);
2672 : float8 result;
2673 :
2674 : /*
2675 : * For tanh, we don't need an errno check because it never overflows.
2676 : */
2677 20 : result = tanh(arg1);
2678 :
2679 20 : if (unlikely(isinf(result)))
2680 0 : float_overflow_error();
2681 :
2682 20 : PG_RETURN_FLOAT8(result);
2683 : }
2684 :
2685 : /*
2686 : * dasinh - returns the inverse hyperbolic sine of arg1
2687 : */
2688 : Datum
2689 20 : dasinh(PG_FUNCTION_ARGS)
2690 : {
2691 20 : float8 arg1 = PG_GETARG_FLOAT8(0);
2692 : float8 result;
2693 :
2694 : /*
2695 : * For asinh, we don't need an errno check because it never overflows.
2696 : */
2697 20 : result = asinh(arg1);
2698 :
2699 20 : PG_RETURN_FLOAT8(result);
2700 : }
2701 :
2702 : /*
2703 : * dacosh - returns the inverse hyperbolic cosine of arg1
2704 : */
2705 : Datum
2706 14 : dacosh(PG_FUNCTION_ARGS)
2707 : {
2708 14 : float8 arg1 = PG_GETARG_FLOAT8(0);
2709 : float8 result;
2710 :
2711 : /*
2712 : * acosh is only defined for inputs >= 1.0. By checking this ourselves,
2713 : * we need not worry about checking for an EDOM error, which is a good
2714 : * thing because some implementations will report that for NaN. Otherwise,
2715 : * no error is possible.
2716 : */
2717 14 : if (arg1 < 1.0)
2718 4 : ereport(ERROR,
2719 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
2720 : errmsg("input is out of range")));
2721 :
2722 10 : result = acosh(arg1);
2723 :
2724 10 : PG_RETURN_FLOAT8(result);
2725 : }
2726 :
2727 : /*
2728 : * datanh - returns the inverse hyperbolic tangent of arg1
2729 : */
2730 : Datum
2731 18 : datanh(PG_FUNCTION_ARGS)
2732 : {
2733 18 : float8 arg1 = PG_GETARG_FLOAT8(0);
2734 : float8 result;
2735 :
2736 : /*
2737 : * atanh is only defined for inputs between -1 and 1. By checking this
2738 : * ourselves, we need not worry about checking for an EDOM error, which is
2739 : * a good thing because some implementations will report that for NaN.
2740 : */
2741 18 : if (arg1 < -1.0 || arg1 > 1.0)
2742 8 : ereport(ERROR,
2743 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
2744 : errmsg("input is out of range")));
2745 :
2746 : /*
2747 : * Also handle the infinity cases ourselves; this is helpful because old
2748 : * glibc versions may produce the wrong errno for this. All other inputs
2749 : * cannot produce an error.
2750 : */
2751 10 : if (arg1 == -1.0)
2752 0 : result = -get_float8_infinity();
2753 10 : else if (arg1 == 1.0)
2754 0 : result = get_float8_infinity();
2755 : else
2756 10 : result = atanh(arg1);
2757 :
2758 10 : PG_RETURN_FLOAT8(result);
2759 : }
2760 :
2761 :
2762 : /* ========== ERROR FUNCTIONS ========== */
2763 :
2764 :
2765 : /*
2766 : * derf - returns the error function: erf(arg1)
2767 : */
2768 : Datum
2769 4088 : derf(PG_FUNCTION_ARGS)
2770 : {
2771 4088 : float8 arg1 = PG_GETARG_FLOAT8(0);
2772 : float8 result;
2773 :
2774 : /*
2775 : * For erf, we don't need an errno check because it never overflows.
2776 : */
2777 4088 : result = erf(arg1);
2778 :
2779 4088 : if (unlikely(isinf(result)))
2780 0 : float_overflow_error();
2781 :
2782 4088 : PG_RETURN_FLOAT8(result);
2783 : }
2784 :
2785 : /*
2786 : * derfc - returns the complementary error function: 1 - erf(arg1)
2787 : */
2788 : Datum
2789 88 : derfc(PG_FUNCTION_ARGS)
2790 : {
2791 88 : float8 arg1 = PG_GETARG_FLOAT8(0);
2792 : float8 result;
2793 :
2794 : /*
2795 : * For erfc, we don't need an errno check because it never overflows.
2796 : */
2797 88 : result = erfc(arg1);
2798 :
2799 88 : if (unlikely(isinf(result)))
2800 0 : float_overflow_error();
2801 :
2802 88 : PG_RETURN_FLOAT8(result);
2803 : }
2804 :
2805 :
2806 : /* ========== GAMMA FUNCTIONS ========== */
2807 :
2808 :
2809 : /*
2810 : * dgamma - returns the gamma function of arg1
2811 : */
2812 : Datum
2813 52 : dgamma(PG_FUNCTION_ARGS)
2814 : {
2815 52 : float8 arg1 = PG_GETARG_FLOAT8(0);
2816 : float8 result;
2817 :
2818 : /*
2819 : * Handle NaN and Inf cases explicitly. This simplifies the overflow
2820 : * checks on platforms that do not set errno.
2821 : */
2822 52 : if (isnan(arg1))
2823 4 : result = arg1;
2824 48 : else if (isinf(arg1))
2825 : {
2826 : /* Per POSIX, an input of -Inf causes a domain error */
2827 8 : if (arg1 < 0)
2828 : {
2829 4 : float_overflow_error();
2830 : result = get_float8_nan(); /* keep compiler quiet */
2831 : }
2832 : else
2833 4 : result = arg1;
2834 : }
2835 : else
2836 : {
2837 : /*
2838 : * Note: the POSIX/C99 gamma function is called "tgamma", not "gamma".
2839 : *
2840 : * On some platforms, tgamma() will not set errno but just return Inf,
2841 : * NaN, or zero to report overflow/underflow; therefore, test those
2842 : * cases explicitly (note that, like the exponential function, the
2843 : * gamma function has no zeros).
2844 : */
2845 40 : errno = 0;
2846 40 : result = tgamma(arg1);
2847 :
2848 40 : if (errno != 0 || isinf(result) || isnan(result))
2849 : {
2850 16 : if (result != 0.0)
2851 12 : float_overflow_error();
2852 : else
2853 4 : float_underflow_error();
2854 : }
2855 24 : else if (result == 0.0)
2856 0 : float_underflow_error();
2857 : }
2858 :
2859 32 : PG_RETURN_FLOAT8(result);
2860 : }
2861 :
2862 :
2863 : /*
2864 : * dlgamma - natural logarithm of absolute value of gamma of arg1
2865 : */
2866 : Datum
2867 59 : dlgamma(PG_FUNCTION_ARGS)
2868 : {
2869 59 : float8 arg1 = PG_GETARG_FLOAT8(0);
2870 : float8 result;
2871 :
2872 : /* On some versions of AIX, lgamma(NaN) fails with ERANGE */
2873 : #if defined(_AIX)
2874 : if (isnan(arg1))
2875 : PG_RETURN_FLOAT8(arg1);
2876 : #endif
2877 :
2878 : /*
2879 : * Note: lgamma may not be thread-safe because it may write to a global
2880 : * variable signgam, which may not be thread-local. However, this doesn't
2881 : * matter to us, since we don't use signgam.
2882 : */
2883 59 : errno = 0;
2884 59 : result = lgamma(arg1);
2885 :
2886 : /*
2887 : * If an ERANGE error occurs, it means there was an overflow or a pole
2888 : * error (which happens for zero and negative integer inputs).
2889 : *
2890 : * On some platforms, lgamma() will not set errno but just return infinity
2891 : * to report overflow, but it should never underflow.
2892 : */
2893 59 : if (errno == ERANGE || (isinf(result) && !isinf(arg1)))
2894 12 : float_overflow_error();
2895 :
2896 47 : PG_RETURN_FLOAT8(result);
2897 : }
2898 :
2899 :
2900 :
2901 : /*
2902 : * =========================
2903 : * FLOAT AGGREGATE OPERATORS
2904 : * =========================
2905 : *
2906 : * float8_accum - accumulate for AVG(), variance aggregates, etc.
2907 : * float4_accum - same, but input data is float4
2908 : * float8_avg - produce final result for float AVG()
2909 : * float8_var_samp - produce final result for float VAR_SAMP()
2910 : * float8_var_pop - produce final result for float VAR_POP()
2911 : * float8_stddev_samp - produce final result for float STDDEV_SAMP()
2912 : * float8_stddev_pop - produce final result for float STDDEV_POP()
2913 : *
2914 : * The naive schoolbook implementation of these aggregates works by
2915 : * accumulating sum(X) and sum(X^2). However, this approach suffers from
2916 : * large rounding errors in the final computation of quantities like the
2917 : * population variance (N*sum(X^2) - sum(X)^2) / N^2, since each of the
2918 : * intermediate terms is potentially very large, while the difference is often
2919 : * quite small.
2920 : *
2921 : * Instead we use the Youngs-Cramer algorithm [1] which works by accumulating
2922 : * Sx=sum(X) and Sxx=sum((X-Sx/N)^2), using a numerically stable algorithm to
2923 : * incrementally update those quantities. The final computations of each of
2924 : * the aggregate values is then trivial and gives more accurate results (for
2925 : * example, the population variance is just Sxx/N). This algorithm is also
2926 : * fairly easy to generalize to allow parallel execution without loss of
2927 : * precision (see, for example, [2]). For more details, and a comparison of
2928 : * this with other algorithms, see [3].
2929 : *
2930 : * The transition datatype for all these aggregates is a 3-element array
2931 : * of float8, holding the values N, Sx, Sxx in that order.
2932 : *
2933 : * Note that we represent N as a float to avoid having to build a special
2934 : * datatype. Given a reasonable floating-point implementation, there should
2935 : * be no accuracy loss unless N exceeds 2 ^ 52 or so (by which time the
2936 : * user will have doubtless lost interest anyway...)
2937 : *
2938 : * [1] Some Results Relevant to Choice of Sum and Sum-of-Product Algorithms,
2939 : * E. A. Youngs and E. M. Cramer, Technometrics Vol 13, No 3, August 1971.
2940 : *
2941 : * [2] Updating Formulae and a Pairwise Algorithm for Computing Sample
2942 : * Variances, T. F. Chan, G. H. Golub & R. J. LeVeque, COMPSTAT 1982.
2943 : *
2944 : * [3] Numerically Stable Parallel Computation of (Co-)Variance, Erich
2945 : * Schubert and Michael Gertz, Proceedings of the 30th International
2946 : * Conference on Scientific and Statistical Database Management, 2018.
2947 : */
2948 :
2949 : static float8 *
2950 33375 : check_float8_array(ArrayType *transarray, const char *caller, int n)
2951 : {
2952 : /*
2953 : * We expect the input to be an N-element float array; verify that. We
2954 : * don't need to use deconstruct_array() since the array data is just
2955 : * going to look like a C array of N float8 values.
2956 : */
2957 33375 : if (ARR_NDIM(transarray) != 1 ||
2958 33375 : ARR_DIMS(transarray)[0] != n ||
2959 33375 : ARR_HASNULL(transarray) ||
2960 33375 : ARR_ELEMTYPE(transarray) != FLOAT8OID)
2961 0 : elog(ERROR, "%s: expected %d-element float8 array", caller, n);
2962 33375 : return (float8 *) ARR_DATA_PTR(transarray);
2963 : }
2964 :
2965 : /*
2966 : * float8_combine
2967 : *
2968 : * An aggregate combine function used to combine two 3 fields
2969 : * aggregate transition data into a single transition data.
2970 : * This function is used only in two stage aggregation and
2971 : * shouldn't be called outside aggregate context.
2972 : */
2973 : Datum
2974 271 : float8_combine(PG_FUNCTION_ARGS)
2975 : {
2976 271 : ArrayType *transarray1 = PG_GETARG_ARRAYTYPE_P(0);
2977 271 : ArrayType *transarray2 = PG_GETARG_ARRAYTYPE_P(1);
2978 : float8 *transvalues1;
2979 : float8 *transvalues2;
2980 : float8 N1,
2981 : Sx1,
2982 : Sxx1,
2983 : N2,
2984 : Sx2,
2985 : Sxx2,
2986 : tmp,
2987 : N,
2988 : Sx,
2989 : Sxx;
2990 :
2991 271 : transvalues1 = check_float8_array(transarray1, "float8_combine", 3);
2992 271 : transvalues2 = check_float8_array(transarray2, "float8_combine", 3);
2993 :
2994 271 : N1 = transvalues1[0];
2995 271 : Sx1 = transvalues1[1];
2996 271 : Sxx1 = transvalues1[2];
2997 :
2998 271 : N2 = transvalues2[0];
2999 271 : Sx2 = transvalues2[1];
3000 271 : Sxx2 = transvalues2[2];
3001 :
3002 : /*--------------------
3003 : * The transition values combine using a generalization of the
3004 : * Youngs-Cramer algorithm as follows:
3005 : *
3006 : * N = N1 + N2
3007 : * Sx = Sx1 + Sx2
3008 : * Sxx = Sxx1 + Sxx2 + N1 * N2 * (Sx1/N1 - Sx2/N2)^2 / N;
3009 : *
3010 : * It's worth handling the special cases N1 = 0 and N2 = 0 separately
3011 : * since those cases are trivial, and we then don't need to worry about
3012 : * division-by-zero errors in the general case.
3013 : *--------------------
3014 : */
3015 271 : if (N1 == 0.0)
3016 : {
3017 261 : N = N2;
3018 261 : Sx = Sx2;
3019 261 : Sxx = Sxx2;
3020 : }
3021 10 : else if (N2 == 0.0)
3022 : {
3023 5 : N = N1;
3024 5 : Sx = Sx1;
3025 5 : Sxx = Sxx1;
3026 : }
3027 : else
3028 : {
3029 5 : N = N1 + N2;
3030 5 : Sx = float8_pl(Sx1, Sx2);
3031 5 : tmp = Sx1 / N1 - Sx2 / N2;
3032 5 : Sxx = Sxx1 + Sxx2 + N1 * N2 * tmp * tmp / N;
3033 5 : if (unlikely(isinf(Sxx)) && !isinf(Sxx1) && !isinf(Sxx2))
3034 0 : float_overflow_error();
3035 : }
3036 :
3037 : /*
3038 : * If we're invoked as an aggregate, we can cheat and modify our first
3039 : * parameter in-place to reduce palloc overhead. Otherwise we construct a
3040 : * new array with the updated transition data and return it.
3041 : */
3042 271 : if (AggCheckCallContext(fcinfo, NULL))
3043 : {
3044 256 : transvalues1[0] = N;
3045 256 : transvalues1[1] = Sx;
3046 256 : transvalues1[2] = Sxx;
3047 :
3048 256 : PG_RETURN_ARRAYTYPE_P(transarray1);
3049 : }
3050 : else
3051 : {
3052 : Datum transdatums[3];
3053 : ArrayType *result;
3054 :
3055 15 : transdatums[0] = Float8GetDatumFast(N);
3056 15 : transdatums[1] = Float8GetDatumFast(Sx);
3057 15 : transdatums[2] = Float8GetDatumFast(Sxx);
3058 :
3059 15 : result = construct_array_builtin(transdatums, 3, FLOAT8OID);
3060 :
3061 15 : PG_RETURN_ARRAYTYPE_P(result);
3062 : }
3063 : }
3064 :
3065 : Datum
3066 31009 : float8_accum(PG_FUNCTION_ARGS)
3067 : {
3068 31009 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3069 31009 : float8 newval = PG_GETARG_FLOAT8(1);
3070 : float8 *transvalues;
3071 : float8 N,
3072 : Sx,
3073 : Sxx,
3074 : tmp;
3075 :
3076 31009 : transvalues = check_float8_array(transarray, "float8_accum", 3);
3077 31009 : N = transvalues[0];
3078 31009 : Sx = transvalues[1];
3079 31009 : Sxx = transvalues[2];
3080 :
3081 : /*
3082 : * Use the Youngs-Cramer algorithm to incorporate the new value into the
3083 : * transition values.
3084 : */
3085 31009 : N += 1.0;
3086 31009 : Sx += newval;
3087 31009 : if (transvalues[0] > 0.0)
3088 : {
3089 30636 : tmp = newval * N - Sx;
3090 30636 : Sxx += tmp * tmp / (N * transvalues[0]);
3091 :
3092 : /*
3093 : * Overflow check. We only report an overflow error when finite
3094 : * inputs lead to infinite results. Note also that Sxx should be NaN
3095 : * if any of the inputs are infinite, so we intentionally prevent Sxx
3096 : * from becoming infinite.
3097 : */
3098 30636 : if (isinf(Sx) || isinf(Sxx))
3099 : {
3100 16 : if (!isinf(transvalues[1]) && !isinf(newval))
3101 0 : float_overflow_error();
3102 :
3103 16 : Sxx = get_float8_nan();
3104 : }
3105 : }
3106 : else
3107 : {
3108 : /*
3109 : * At the first input, we normally can leave Sxx as 0. However, if
3110 : * the first input is Inf or NaN, we'd better force Sxx to NaN;
3111 : * otherwise we will falsely report variance zero when there are no
3112 : * more inputs.
3113 : */
3114 373 : if (isnan(newval) || isinf(newval))
3115 32 : Sxx = get_float8_nan();
3116 : }
3117 :
3118 : /*
3119 : * If we're invoked as an aggregate, we can cheat and modify our first
3120 : * parameter in-place to reduce palloc overhead. Otherwise we construct a
3121 : * new array with the updated transition data and return it.
3122 : */
3123 31009 : if (AggCheckCallContext(fcinfo, NULL))
3124 : {
3125 31004 : transvalues[0] = N;
3126 31004 : transvalues[1] = Sx;
3127 31004 : transvalues[2] = Sxx;
3128 :
3129 31004 : PG_RETURN_ARRAYTYPE_P(transarray);
3130 : }
3131 : else
3132 : {
3133 : Datum transdatums[3];
3134 : ArrayType *result;
3135 :
3136 5 : transdatums[0] = Float8GetDatumFast(N);
3137 5 : transdatums[1] = Float8GetDatumFast(Sx);
3138 5 : transdatums[2] = Float8GetDatumFast(Sxx);
3139 :
3140 5 : result = construct_array_builtin(transdatums, 3, FLOAT8OID);
3141 :
3142 5 : PG_RETURN_ARRAYTYPE_P(result);
3143 : }
3144 : }
3145 :
3146 : Datum
3147 192 : float4_accum(PG_FUNCTION_ARGS)
3148 : {
3149 192 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3150 :
3151 : /* do computations as float8 */
3152 192 : float8 newval = PG_GETARG_FLOAT4(1);
3153 : float8 *transvalues;
3154 : float8 N,
3155 : Sx,
3156 : Sxx,
3157 : tmp;
3158 :
3159 192 : transvalues = check_float8_array(transarray, "float4_accum", 3);
3160 192 : N = transvalues[0];
3161 192 : Sx = transvalues[1];
3162 192 : Sxx = transvalues[2];
3163 :
3164 : /*
3165 : * Use the Youngs-Cramer algorithm to incorporate the new value into the
3166 : * transition values.
3167 : */
3168 192 : N += 1.0;
3169 192 : Sx += newval;
3170 192 : if (transvalues[0] > 0.0)
3171 : {
3172 136 : tmp = newval * N - Sx;
3173 136 : Sxx += tmp * tmp / (N * transvalues[0]);
3174 :
3175 : /*
3176 : * Overflow check. We only report an overflow error when finite
3177 : * inputs lead to infinite results. Note also that Sxx should be NaN
3178 : * if any of the inputs are infinite, so we intentionally prevent Sxx
3179 : * from becoming infinite.
3180 : */
3181 136 : if (isinf(Sx) || isinf(Sxx))
3182 : {
3183 0 : if (!isinf(transvalues[1]) && !isinf(newval))
3184 0 : float_overflow_error();
3185 :
3186 0 : Sxx = get_float8_nan();
3187 : }
3188 : }
3189 : else
3190 : {
3191 : /*
3192 : * At the first input, we normally can leave Sxx as 0. However, if
3193 : * the first input is Inf or NaN, we'd better force Sxx to NaN;
3194 : * otherwise we will falsely report variance zero when there are no
3195 : * more inputs.
3196 : */
3197 56 : if (isnan(newval) || isinf(newval))
3198 16 : Sxx = get_float8_nan();
3199 : }
3200 :
3201 : /*
3202 : * If we're invoked as an aggregate, we can cheat and modify our first
3203 : * parameter in-place to reduce palloc overhead. Otherwise we construct a
3204 : * new array with the updated transition data and return it.
3205 : */
3206 192 : if (AggCheckCallContext(fcinfo, NULL))
3207 : {
3208 192 : transvalues[0] = N;
3209 192 : transvalues[1] = Sx;
3210 192 : transvalues[2] = Sxx;
3211 :
3212 192 : PG_RETURN_ARRAYTYPE_P(transarray);
3213 : }
3214 : else
3215 : {
3216 : Datum transdatums[3];
3217 : ArrayType *result;
3218 :
3219 0 : transdatums[0] = Float8GetDatumFast(N);
3220 0 : transdatums[1] = Float8GetDatumFast(Sx);
3221 0 : transdatums[2] = Float8GetDatumFast(Sxx);
3222 :
3223 0 : result = construct_array_builtin(transdatums, 3, FLOAT8OID);
3224 :
3225 0 : PG_RETURN_ARRAYTYPE_P(result);
3226 : }
3227 : }
3228 :
3229 : Datum
3230 337 : float8_avg(PG_FUNCTION_ARGS)
3231 : {
3232 337 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3233 : float8 *transvalues;
3234 : float8 N,
3235 : Sx;
3236 :
3237 337 : transvalues = check_float8_array(transarray, "float8_avg", 3);
3238 337 : N = transvalues[0];
3239 337 : Sx = transvalues[1];
3240 : /* ignore Sxx */
3241 :
3242 : /* SQL defines AVG of no values to be NULL */
3243 337 : if (N == 0.0)
3244 8 : PG_RETURN_NULL();
3245 :
3246 329 : PG_RETURN_FLOAT8(Sx / N);
3247 : }
3248 :
3249 : Datum
3250 56 : float8_var_pop(PG_FUNCTION_ARGS)
3251 : {
3252 56 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3253 : float8 *transvalues;
3254 : float8 N,
3255 : Sxx;
3256 :
3257 56 : transvalues = check_float8_array(transarray, "float8_var_pop", 3);
3258 56 : N = transvalues[0];
3259 : /* ignore Sx */
3260 56 : Sxx = transvalues[2];
3261 :
3262 : /* Population variance is undefined when N is 0, so return NULL */
3263 56 : if (N == 0.0)
3264 0 : PG_RETURN_NULL();
3265 :
3266 : /* Note that Sxx is guaranteed to be non-negative */
3267 :
3268 56 : PG_RETURN_FLOAT8(Sxx / N);
3269 : }
3270 :
3271 : Datum
3272 28 : float8_var_samp(PG_FUNCTION_ARGS)
3273 : {
3274 28 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3275 : float8 *transvalues;
3276 : float8 N,
3277 : Sxx;
3278 :
3279 28 : transvalues = check_float8_array(transarray, "float8_var_samp", 3);
3280 28 : N = transvalues[0];
3281 : /* ignore Sx */
3282 28 : Sxx = transvalues[2];
3283 :
3284 : /* Sample variance is undefined when N is 0 or 1, so return NULL */
3285 28 : if (N <= 1.0)
3286 24 : PG_RETURN_NULL();
3287 :
3288 : /* Note that Sxx is guaranteed to be non-negative */
3289 :
3290 4 : PG_RETURN_FLOAT8(Sxx / (N - 1.0));
3291 : }
3292 :
3293 : Datum
3294 28 : float8_stddev_pop(PG_FUNCTION_ARGS)
3295 : {
3296 28 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3297 : float8 *transvalues;
3298 : float8 N,
3299 : Sxx;
3300 :
3301 28 : transvalues = check_float8_array(transarray, "float8_stddev_pop", 3);
3302 28 : N = transvalues[0];
3303 : /* ignore Sx */
3304 28 : Sxx = transvalues[2];
3305 :
3306 : /* Population stddev is undefined when N is 0, so return NULL */
3307 28 : if (N == 0.0)
3308 0 : PG_RETURN_NULL();
3309 :
3310 : /* Note that Sxx is guaranteed to be non-negative */
3311 :
3312 28 : PG_RETURN_FLOAT8(sqrt(Sxx / N));
3313 : }
3314 :
3315 : Datum
3316 32 : float8_stddev_samp(PG_FUNCTION_ARGS)
3317 : {
3318 32 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3319 : float8 *transvalues;
3320 : float8 N,
3321 : Sxx;
3322 :
3323 32 : transvalues = check_float8_array(transarray, "float8_stddev_samp", 3);
3324 32 : N = transvalues[0];
3325 : /* ignore Sx */
3326 32 : Sxx = transvalues[2];
3327 :
3328 : /* Sample stddev is undefined when N is 0 or 1, so return NULL */
3329 32 : if (N <= 1.0)
3330 24 : PG_RETURN_NULL();
3331 :
3332 : /* Note that Sxx is guaranteed to be non-negative */
3333 :
3334 8 : PG_RETURN_FLOAT8(sqrt(Sxx / (N - 1.0)));
3335 : }
3336 :
3337 : /*
3338 : * =========================
3339 : * SQL2003 BINARY AGGREGATES
3340 : * =========================
3341 : *
3342 : * As with the preceding aggregates, we use the Youngs-Cramer algorithm to
3343 : * reduce rounding errors in the aggregate final functions.
3344 : *
3345 : * The transition datatype for all these aggregates is an 8-element array of
3346 : * float8, holding the values N, Sx=sum(X), Sxx=sum((X-Sx/N)^2), Sy=sum(Y),
3347 : * Syy=sum((Y-Sy/N)^2), Sxy=sum((X-Sx/N)*(Y-Sy/N)), commonX, and commonY
3348 : * in that order.
3349 : *
3350 : * commonX is defined as the common X value if all the X values were the same,
3351 : * else NaN; likewise for commonY. This is useful for deciding whether corr()
3352 : * and related functions should return NULL. This representation cannot
3353 : * distinguish the-values-were-all-NaN from the-values-were-not-all-the-same,
3354 : * but that's okay because for this purpose we use the IEEE float arithmetic
3355 : * principle that two NaNs are never equal. The SQL standard doesn't mention
3356 : * NaNs, but it says that NULL is to be returned when N*sum(X*X) equals
3357 : * sum(X)*sum(X) (etc), and that shouldn't be considered true for NaNs.
3358 : * Testing this as written in the spec would be highly subject to roundoff
3359 : * error, so instead we directly track whether all the inputs are equal.
3360 : *
3361 : * Note that Y is the first argument to all these aggregates!
3362 : *
3363 : * It might seem attractive to optimize this by having multiple accumulator
3364 : * functions that only calculate the sums actually needed. But on most
3365 : * modern machines, a couple of extra floating-point multiplies will be
3366 : * insignificant compared to the other per-tuple overhead, so I've chosen
3367 : * to minimize code space instead.
3368 : */
3369 :
3370 : Datum
3371 957 : float8_regr_accum(PG_FUNCTION_ARGS)
3372 : {
3373 957 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3374 957 : float8 newvalY = PG_GETARG_FLOAT8(1);
3375 957 : float8 newvalX = PG_GETARG_FLOAT8(2);
3376 : float8 *transvalues;
3377 : float8 N,
3378 : Sx,
3379 : Sxx,
3380 : Sy,
3381 : Syy,
3382 : Sxy,
3383 : commonX,
3384 : commonY,
3385 : tmpX,
3386 : tmpY,
3387 : scale;
3388 :
3389 957 : transvalues = check_float8_array(transarray, "float8_regr_accum", 8);
3390 957 : N = transvalues[0];
3391 957 : Sx = transvalues[1];
3392 957 : Sxx = transvalues[2];
3393 957 : Sy = transvalues[3];
3394 957 : Syy = transvalues[4];
3395 957 : Sxy = transvalues[5];
3396 957 : commonX = transvalues[6];
3397 957 : commonY = transvalues[7];
3398 :
3399 : /*
3400 : * Use the Youngs-Cramer algorithm to incorporate the new values into the
3401 : * transition values.
3402 : */
3403 957 : N += 1.0;
3404 957 : Sx += newvalX;
3405 957 : Sy += newvalY;
3406 957 : if (transvalues[0] > 0.0)
3407 : {
3408 : /*
3409 : * Check to see if we have seen distinct inputs. We can use a test
3410 : * that's a bit cheaper than float8_ne() because if commonX is already
3411 : * NaN, it does not matter whether the != test returns true or not.
3412 : */
3413 853 : if (newvalX != commonX || isnan(newvalX))
3414 732 : commonX = get_float8_nan();
3415 853 : if (newvalY != commonY || isnan(newvalY))
3416 621 : commonY = get_float8_nan();
3417 :
3418 853 : tmpX = newvalX * N - Sx;
3419 853 : tmpY = newvalY * N - Sy;
3420 853 : scale = 1.0 / (N * transvalues[0]);
3421 :
3422 : /*
3423 : * If we have not seen distinct inputs, then Sxx, Syy, and/or Sxy
3424 : * should remain zero (since Sx's exact value would be N * commonX,
3425 : * etc). Updating them would just create the possibility of injecting
3426 : * roundoff error, and we need exact zero results so that the final
3427 : * functions will return NULL in the right cases.
3428 : */
3429 853 : if (isnan(commonX))
3430 732 : Sxx += tmpX * tmpX * scale;
3431 853 : if (isnan(commonY))
3432 621 : Syy += tmpY * tmpY * scale;
3433 853 : if (isnan(commonX) && isnan(commonY))
3434 500 : Sxy += tmpX * tmpY * scale;
3435 :
3436 : /*
3437 : * Overflow check. We only report an overflow error when finite
3438 : * inputs lead to infinite results. Note also that Sxx, Syy and Sxy
3439 : * should be NaN if any of the relevant inputs are infinite, so we
3440 : * intentionally prevent them from becoming infinite.
3441 : */
3442 853 : if (isinf(Sx) || isinf(Sxx) || isinf(Sy) || isinf(Syy) || isinf(Sxy))
3443 : {
3444 0 : if (((isinf(Sx) || isinf(Sxx)) &&
3445 0 : !isinf(transvalues[1]) && !isinf(newvalX)) ||
3446 0 : ((isinf(Sy) || isinf(Syy)) &&
3447 0 : !isinf(transvalues[3]) && !isinf(newvalY)) ||
3448 0 : (isinf(Sxy) &&
3449 0 : !isinf(transvalues[1]) && !isinf(newvalX) &&
3450 0 : !isinf(transvalues[3]) && !isinf(newvalY)))
3451 0 : float_overflow_error();
3452 :
3453 0 : if (isinf(Sxx))
3454 0 : Sxx = get_float8_nan();
3455 0 : if (isinf(Syy))
3456 0 : Syy = get_float8_nan();
3457 0 : if (isinf(Sxy))
3458 0 : Sxy = get_float8_nan();
3459 : }
3460 : }
3461 : else
3462 : {
3463 : /*
3464 : * At the first input, we normally can leave Sxx et al as 0. However,
3465 : * if the first input is Inf or NaN, we'd better force the dependent
3466 : * sums to NaN; otherwise we will falsely report variance zero when
3467 : * there are no more inputs.
3468 : */
3469 104 : if (isnan(newvalX) || isinf(newvalX))
3470 28 : Sxx = Sxy = get_float8_nan();
3471 104 : if (isnan(newvalY) || isinf(newvalY))
3472 4 : Syy = Sxy = get_float8_nan();
3473 :
3474 104 : commonX = newvalX;
3475 104 : commonY = newvalY;
3476 : }
3477 :
3478 : /*
3479 : * If we're invoked as an aggregate, we can cheat and modify our first
3480 : * parameter in-place to reduce palloc overhead. Otherwise we construct a
3481 : * new array with the updated transition data and return it.
3482 : */
3483 957 : if (AggCheckCallContext(fcinfo, NULL))
3484 : {
3485 952 : transvalues[0] = N;
3486 952 : transvalues[1] = Sx;
3487 952 : transvalues[2] = Sxx;
3488 952 : transvalues[3] = Sy;
3489 952 : transvalues[4] = Syy;
3490 952 : transvalues[5] = Sxy;
3491 952 : transvalues[6] = commonX;
3492 952 : transvalues[7] = commonY;
3493 :
3494 952 : PG_RETURN_ARRAYTYPE_P(transarray);
3495 : }
3496 : else
3497 : {
3498 : Datum transdatums[8];
3499 : ArrayType *result;
3500 :
3501 5 : transdatums[0] = Float8GetDatumFast(N);
3502 5 : transdatums[1] = Float8GetDatumFast(Sx);
3503 5 : transdatums[2] = Float8GetDatumFast(Sxx);
3504 5 : transdatums[3] = Float8GetDatumFast(Sy);
3505 5 : transdatums[4] = Float8GetDatumFast(Syy);
3506 5 : transdatums[5] = Float8GetDatumFast(Sxy);
3507 5 : transdatums[6] = Float8GetDatumFast(commonX);
3508 5 : transdatums[7] = Float8GetDatumFast(commonY);
3509 :
3510 5 : result = construct_array_builtin(transdatums, 8, FLOAT8OID);
3511 :
3512 5 : PG_RETURN_ARRAYTYPE_P(result);
3513 : }
3514 : }
3515 :
3516 : /*
3517 : * float8_regr_combine
3518 : *
3519 : * An aggregate combine function used to combine two 8-fields
3520 : * aggregate transition data into a single transition data.
3521 : * This function is used only in two stage aggregation and
3522 : * shouldn't be called outside aggregate context.
3523 : */
3524 : Datum
3525 15 : float8_regr_combine(PG_FUNCTION_ARGS)
3526 : {
3527 15 : ArrayType *transarray1 = PG_GETARG_ARRAYTYPE_P(0);
3528 15 : ArrayType *transarray2 = PG_GETARG_ARRAYTYPE_P(1);
3529 : float8 *transvalues1;
3530 : float8 *transvalues2;
3531 : float8 N1,
3532 : Sx1,
3533 : Sxx1,
3534 : Sy1,
3535 : Syy1,
3536 : Sxy1,
3537 : Cx1,
3538 : Cy1,
3539 : N2,
3540 : Sx2,
3541 : Sxx2,
3542 : Sy2,
3543 : Syy2,
3544 : Sxy2,
3545 : Cx2,
3546 : Cy2,
3547 : tmp1,
3548 : tmp2,
3549 : N,
3550 : Sx,
3551 : Sxx,
3552 : Sy,
3553 : Syy,
3554 : Sxy,
3555 : Cx,
3556 : Cy;
3557 :
3558 15 : transvalues1 = check_float8_array(transarray1, "float8_regr_combine", 8);
3559 15 : transvalues2 = check_float8_array(transarray2, "float8_regr_combine", 8);
3560 :
3561 15 : N1 = transvalues1[0];
3562 15 : Sx1 = transvalues1[1];
3563 15 : Sxx1 = transvalues1[2];
3564 15 : Sy1 = transvalues1[3];
3565 15 : Syy1 = transvalues1[4];
3566 15 : Sxy1 = transvalues1[5];
3567 15 : Cx1 = transvalues1[6];
3568 15 : Cy1 = transvalues1[7];
3569 :
3570 15 : N2 = transvalues2[0];
3571 15 : Sx2 = transvalues2[1];
3572 15 : Sxx2 = transvalues2[2];
3573 15 : Sy2 = transvalues2[3];
3574 15 : Syy2 = transvalues2[4];
3575 15 : Sxy2 = transvalues2[5];
3576 15 : Cx2 = transvalues2[6];
3577 15 : Cy2 = transvalues2[7];
3578 :
3579 : /*--------------------
3580 : * The transition values combine using a generalization of the
3581 : * Youngs-Cramer algorithm as follows:
3582 : *
3583 : * N = N1 + N2
3584 : * Sx = Sx1 + Sx2
3585 : * Sxx = Sxx1 + Sxx2 + N1 * N2 * (Sx1/N1 - Sx2/N2)^2 / N
3586 : * Sy = Sy1 + Sy2
3587 : * Syy = Syy1 + Syy2 + N1 * N2 * (Sy1/N1 - Sy2/N2)^2 / N
3588 : * Sxy = Sxy1 + Sxy2 + N1 * N2 * (Sx1/N1 - Sx2/N2) * (Sy1/N1 - Sy2/N2) / N
3589 : *
3590 : * It's worth handling the special cases N1 = 0 and N2 = 0 separately
3591 : * since those cases are trivial, and we then don't need to worry about
3592 : * division-by-zero errors in the general case.
3593 : *--------------------
3594 : */
3595 15 : if (N1 == 0.0)
3596 : {
3597 5 : N = N2;
3598 5 : Sx = Sx2;
3599 5 : Sxx = Sxx2;
3600 5 : Sy = Sy2;
3601 5 : Syy = Syy2;
3602 5 : Sxy = Sxy2;
3603 5 : Cx = Cx2;
3604 5 : Cy = Cy2;
3605 : }
3606 10 : else if (N2 == 0.0)
3607 : {
3608 5 : N = N1;
3609 5 : Sx = Sx1;
3610 5 : Sxx = Sxx1;
3611 5 : Sy = Sy1;
3612 5 : Syy = Syy1;
3613 5 : Sxy = Sxy1;
3614 5 : Cx = Cx1;
3615 5 : Cy = Cy1;
3616 : }
3617 : else
3618 : {
3619 5 : N = N1 + N2;
3620 5 : Sx = float8_pl(Sx1, Sx2);
3621 5 : tmp1 = Sx1 / N1 - Sx2 / N2;
3622 5 : Sxx = Sxx1 + Sxx2 + N1 * N2 * tmp1 * tmp1 / N;
3623 5 : if (unlikely(isinf(Sxx)) && !isinf(Sxx1) && !isinf(Sxx2))
3624 0 : float_overflow_error();
3625 5 : Sy = float8_pl(Sy1, Sy2);
3626 5 : tmp2 = Sy1 / N1 - Sy2 / N2;
3627 5 : Syy = Syy1 + Syy2 + N1 * N2 * tmp2 * tmp2 / N;
3628 5 : if (unlikely(isinf(Syy)) && !isinf(Syy1) && !isinf(Syy2))
3629 0 : float_overflow_error();
3630 5 : Sxy = Sxy1 + Sxy2 + N1 * N2 * tmp1 * tmp2 / N;
3631 5 : if (unlikely(isinf(Sxy)) && !isinf(Sxy1) && !isinf(Sxy2))
3632 0 : float_overflow_error();
3633 5 : if (float8_eq(Cx1, Cx2))
3634 5 : Cx = Cx1;
3635 : else
3636 0 : Cx = get_float8_nan();
3637 5 : if (float8_eq(Cy1, Cy2))
3638 0 : Cy = Cy1;
3639 : else
3640 5 : Cy = get_float8_nan();
3641 : }
3642 :
3643 : /*
3644 : * If we're invoked as an aggregate, we can cheat and modify our first
3645 : * parameter in-place to reduce palloc overhead. Otherwise we construct a
3646 : * new array with the updated transition data and return it.
3647 : */
3648 15 : if (AggCheckCallContext(fcinfo, NULL))
3649 : {
3650 0 : transvalues1[0] = N;
3651 0 : transvalues1[1] = Sx;
3652 0 : transvalues1[2] = Sxx;
3653 0 : transvalues1[3] = Sy;
3654 0 : transvalues1[4] = Syy;
3655 0 : transvalues1[5] = Sxy;
3656 0 : transvalues1[6] = Cx;
3657 0 : transvalues1[7] = Cy;
3658 :
3659 0 : PG_RETURN_ARRAYTYPE_P(transarray1);
3660 : }
3661 : else
3662 : {
3663 : Datum transdatums[8];
3664 : ArrayType *result;
3665 :
3666 15 : transdatums[0] = Float8GetDatumFast(N);
3667 15 : transdatums[1] = Float8GetDatumFast(Sx);
3668 15 : transdatums[2] = Float8GetDatumFast(Sxx);
3669 15 : transdatums[3] = Float8GetDatumFast(Sy);
3670 15 : transdatums[4] = Float8GetDatumFast(Syy);
3671 15 : transdatums[5] = Float8GetDatumFast(Sxy);
3672 15 : transdatums[6] = Float8GetDatumFast(Cx);
3673 15 : transdatums[7] = Float8GetDatumFast(Cy);
3674 :
3675 15 : result = construct_array_builtin(transdatums, 8, FLOAT8OID);
3676 :
3677 15 : PG_RETURN_ARRAYTYPE_P(result);
3678 : }
3679 : }
3680 :
3681 :
3682 : Datum
3683 20 : float8_regr_sxx(PG_FUNCTION_ARGS)
3684 : {
3685 20 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3686 : float8 *transvalues;
3687 : float8 N,
3688 : Sxx;
3689 :
3690 20 : transvalues = check_float8_array(transarray, "float8_regr_sxx", 8);
3691 20 : N = transvalues[0];
3692 20 : Sxx = transvalues[2];
3693 :
3694 : /* if N is 0 we should return NULL */
3695 20 : if (N < 1.0)
3696 0 : PG_RETURN_NULL();
3697 :
3698 : /* Note that Sxx is guaranteed to be non-negative */
3699 :
3700 20 : PG_RETURN_FLOAT8(Sxx);
3701 : }
3702 :
3703 : Datum
3704 20 : float8_regr_syy(PG_FUNCTION_ARGS)
3705 : {
3706 20 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3707 : float8 *transvalues;
3708 : float8 N,
3709 : Syy;
3710 :
3711 20 : transvalues = check_float8_array(transarray, "float8_regr_syy", 8);
3712 20 : N = transvalues[0];
3713 20 : Syy = transvalues[4];
3714 :
3715 : /* if N is 0 we should return NULL */
3716 20 : if (N < 1.0)
3717 0 : PG_RETURN_NULL();
3718 :
3719 : /* Note that Syy is guaranteed to be non-negative */
3720 :
3721 20 : PG_RETURN_FLOAT8(Syy);
3722 : }
3723 :
3724 : Datum
3725 20 : float8_regr_sxy(PG_FUNCTION_ARGS)
3726 : {
3727 20 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3728 : float8 *transvalues;
3729 : float8 N,
3730 : Sxy;
3731 :
3732 20 : transvalues = check_float8_array(transarray, "float8_regr_sxy", 8);
3733 20 : N = transvalues[0];
3734 20 : Sxy = transvalues[5];
3735 :
3736 : /* if N is 0 we should return NULL */
3737 20 : if (N < 1.0)
3738 0 : PG_RETURN_NULL();
3739 :
3740 : /* A negative result is valid here */
3741 :
3742 20 : PG_RETURN_FLOAT8(Sxy);
3743 : }
3744 :
3745 : Datum
3746 4 : float8_regr_avgx(PG_FUNCTION_ARGS)
3747 : {
3748 4 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3749 : float8 *transvalues;
3750 : float8 N,
3751 : Sx,
3752 : commonX;
3753 :
3754 4 : transvalues = check_float8_array(transarray, "float8_regr_avgx", 8);
3755 4 : N = transvalues[0];
3756 4 : Sx = transvalues[1];
3757 4 : commonX = transvalues[6];
3758 :
3759 : /* if N is 0 we should return NULL */
3760 4 : if (N < 1.0)
3761 0 : PG_RETURN_NULL();
3762 :
3763 : /* if all inputs were the same just return that, avoiding roundoff error */
3764 4 : if (!isnan(commonX))
3765 0 : PG_RETURN_FLOAT8(commonX);
3766 :
3767 4 : PG_RETURN_FLOAT8(Sx / N);
3768 : }
3769 :
3770 : Datum
3771 4 : float8_regr_avgy(PG_FUNCTION_ARGS)
3772 : {
3773 4 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3774 : float8 *transvalues;
3775 : float8 N,
3776 : Sy,
3777 : commonY;
3778 :
3779 4 : transvalues = check_float8_array(transarray, "float8_regr_avgy", 8);
3780 4 : N = transvalues[0];
3781 4 : Sy = transvalues[3];
3782 4 : commonY = transvalues[7];
3783 :
3784 : /* if N is 0 we should return NULL */
3785 4 : if (N < 1.0)
3786 0 : PG_RETURN_NULL();
3787 :
3788 : /* if all inputs were the same just return that, avoiding roundoff error */
3789 4 : if (!isnan(commonY))
3790 0 : PG_RETURN_FLOAT8(commonY);
3791 :
3792 4 : PG_RETURN_FLOAT8(Sy / N);
3793 : }
3794 :
3795 : Datum
3796 16 : float8_covar_pop(PG_FUNCTION_ARGS)
3797 : {
3798 16 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3799 : float8 *transvalues;
3800 : float8 N,
3801 : Sxy;
3802 :
3803 16 : transvalues = check_float8_array(transarray, "float8_covar_pop", 8);
3804 16 : N = transvalues[0];
3805 16 : Sxy = transvalues[5];
3806 :
3807 : /* if N is 0 we should return NULL */
3808 16 : if (N < 1.0)
3809 0 : PG_RETURN_NULL();
3810 :
3811 16 : PG_RETURN_FLOAT8(Sxy / N);
3812 : }
3813 :
3814 : Datum
3815 16 : float8_covar_samp(PG_FUNCTION_ARGS)
3816 : {
3817 16 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3818 : float8 *transvalues;
3819 : float8 N,
3820 : Sxy;
3821 :
3822 16 : transvalues = check_float8_array(transarray, "float8_covar_samp", 8);
3823 16 : N = transvalues[0];
3824 16 : Sxy = transvalues[5];
3825 :
3826 : /* if N is <= 1 we should return NULL */
3827 16 : if (N < 2.0)
3828 12 : PG_RETURN_NULL();
3829 :
3830 4 : PG_RETURN_FLOAT8(Sxy / (N - 1.0));
3831 : }
3832 :
3833 : Datum
3834 36 : float8_corr(PG_FUNCTION_ARGS)
3835 : {
3836 36 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3837 : float8 *transvalues;
3838 : float8 N,
3839 : Sxx,
3840 : Syy,
3841 : Sxy,
3842 : product,
3843 : sqrtproduct,
3844 : result;
3845 :
3846 36 : transvalues = check_float8_array(transarray, "float8_corr", 8);
3847 36 : N = transvalues[0];
3848 36 : Sxx = transvalues[2];
3849 36 : Syy = transvalues[4];
3850 36 : Sxy = transvalues[5];
3851 :
3852 : /* if N is 0 we should return NULL */
3853 36 : if (N < 1.0)
3854 0 : PG_RETURN_NULL();
3855 :
3856 : /* Note that Sxx and Syy are guaranteed to be non-negative */
3857 :
3858 : /* per spec, return NULL for horizontal and vertical lines */
3859 36 : if (Sxx == 0 || Syy == 0)
3860 16 : PG_RETURN_NULL();
3861 :
3862 : /*
3863 : * The product Sxx * Syy might underflow or overflow. If so, we can
3864 : * recover by computing sqrt(Sxx) * sqrt(Syy) instead of sqrt(Sxx * Syy).
3865 : * However, the double sqrt() calculation is a bit slower and less
3866 : * accurate, so don't do it if we don't have to.
3867 : */
3868 20 : product = Sxx * Syy;
3869 20 : if (product == 0 || isinf(product))
3870 8 : sqrtproduct = sqrt(Sxx) * sqrt(Syy);
3871 : else
3872 12 : sqrtproduct = sqrt(product);
3873 20 : result = Sxy / sqrtproduct;
3874 :
3875 : /*
3876 : * Despite all these precautions, this formula can yield results outside
3877 : * [-1, 1] due to roundoff error. Clamp it to the expected range.
3878 : */
3879 20 : if (result < -1)
3880 0 : result = -1;
3881 20 : else if (result > 1)
3882 4 : result = 1;
3883 :
3884 20 : PG_RETURN_FLOAT8(result);
3885 : }
3886 :
3887 : Datum
3888 12 : float8_regr_r2(PG_FUNCTION_ARGS)
3889 : {
3890 12 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3891 : float8 *transvalues;
3892 : float8 N,
3893 : Sxx,
3894 : Syy,
3895 : Sxy;
3896 :
3897 12 : transvalues = check_float8_array(transarray, "float8_regr_r2", 8);
3898 12 : N = transvalues[0];
3899 12 : Sxx = transvalues[2];
3900 12 : Syy = transvalues[4];
3901 12 : Sxy = transvalues[5];
3902 :
3903 : /* if N is 0 we should return NULL */
3904 12 : if (N < 1.0)
3905 0 : PG_RETURN_NULL();
3906 :
3907 : /* Note that Sxx and Syy are guaranteed to be non-negative */
3908 :
3909 : /* per spec, return NULL for a vertical line */
3910 12 : if (Sxx == 0)
3911 4 : PG_RETURN_NULL();
3912 :
3913 : /* per spec, return 1.0 for a horizontal line */
3914 8 : if (Syy == 0)
3915 4 : PG_RETURN_FLOAT8(1.0);
3916 :
3917 4 : PG_RETURN_FLOAT8((Sxy * Sxy) / (Sxx * Syy));
3918 : }
3919 :
3920 : Datum
3921 8 : float8_regr_slope(PG_FUNCTION_ARGS)
3922 : {
3923 8 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3924 : float8 *transvalues;
3925 : float8 N,
3926 : Sxx,
3927 : Sxy;
3928 :
3929 8 : transvalues = check_float8_array(transarray, "float8_regr_slope", 8);
3930 8 : N = transvalues[0];
3931 8 : Sxx = transvalues[2];
3932 8 : Sxy = transvalues[5];
3933 :
3934 : /* if N is 0 we should return NULL */
3935 8 : if (N < 1.0)
3936 0 : PG_RETURN_NULL();
3937 :
3938 : /* Note that Sxx is guaranteed to be non-negative */
3939 :
3940 : /* per spec, return NULL for a vertical line */
3941 8 : if (Sxx == 0)
3942 4 : PG_RETURN_NULL();
3943 :
3944 4 : PG_RETURN_FLOAT8(Sxy / Sxx);
3945 : }
3946 :
3947 : Datum
3948 8 : float8_regr_intercept(PG_FUNCTION_ARGS)
3949 : {
3950 8 : ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
3951 : float8 *transvalues;
3952 : float8 N,
3953 : Sx,
3954 : Sxx,
3955 : Sy,
3956 : Sxy;
3957 :
3958 8 : transvalues = check_float8_array(transarray, "float8_regr_intercept", 8);
3959 8 : N = transvalues[0];
3960 8 : Sx = transvalues[1];
3961 8 : Sxx = transvalues[2];
3962 8 : Sy = transvalues[3];
3963 8 : Sxy = transvalues[5];
3964 :
3965 : /* if N is 0 we should return NULL */
3966 8 : if (N < 1.0)
3967 0 : PG_RETURN_NULL();
3968 :
3969 : /* Note that Sxx is guaranteed to be non-negative */
3970 :
3971 : /* per spec, return NULL for a vertical line */
3972 8 : if (Sxx == 0)
3973 4 : PG_RETURN_NULL();
3974 :
3975 4 : PG_RETURN_FLOAT8((Sy - Sx * Sxy / Sxx) / N);
3976 : }
3977 :
3978 :
3979 : /*
3980 : * ====================================
3981 : * MIXED-PRECISION ARITHMETIC OPERATORS
3982 : * ====================================
3983 : */
3984 :
3985 : /*
3986 : * float48pl - returns arg1 + arg2
3987 : * float48mi - returns arg1 - arg2
3988 : * float48mul - returns arg1 * arg2
3989 : * float48div - returns arg1 / arg2
3990 : */
3991 : Datum
3992 17 : float48pl(PG_FUNCTION_ARGS)
3993 : {
3994 17 : float4 arg1 = PG_GETARG_FLOAT4(0);
3995 17 : float8 arg2 = PG_GETARG_FLOAT8(1);
3996 :
3997 17 : PG_RETURN_FLOAT8(float8_pl((float8) arg1, arg2));
3998 : }
3999 :
4000 : Datum
4001 4 : float48mi(PG_FUNCTION_ARGS)
4002 : {
4003 4 : float4 arg1 = PG_GETARG_FLOAT4(0);
4004 4 : float8 arg2 = PG_GETARG_FLOAT8(1);
4005 :
4006 4 : PG_RETURN_FLOAT8(float8_mi((float8) arg1, arg2));
4007 : }
4008 :
4009 : Datum
4010 4 : float48mul(PG_FUNCTION_ARGS)
4011 : {
4012 4 : float4 arg1 = PG_GETARG_FLOAT4(0);
4013 4 : float8 arg2 = PG_GETARG_FLOAT8(1);
4014 :
4015 4 : PG_RETURN_FLOAT8(float8_mul((float8) arg1, arg2));
4016 : }
4017 :
4018 : Datum
4019 4 : float48div(PG_FUNCTION_ARGS)
4020 : {
4021 4 : float4 arg1 = PG_GETARG_FLOAT4(0);
4022 4 : float8 arg2 = PG_GETARG_FLOAT8(1);
4023 :
4024 4 : PG_RETURN_FLOAT8(float8_div((float8) arg1, arg2));
4025 : }
4026 :
4027 : /*
4028 : * float84pl - returns arg1 + arg2
4029 : * float84mi - returns arg1 - arg2
4030 : * float84mul - returns arg1 * arg2
4031 : * float84div - returns arg1 / arg2
4032 : */
4033 : Datum
4034 8 : float84pl(PG_FUNCTION_ARGS)
4035 : {
4036 8 : float8 arg1 = PG_GETARG_FLOAT8(0);
4037 8 : float4 arg2 = PG_GETARG_FLOAT4(1);
4038 :
4039 8 : PG_RETURN_FLOAT8(float8_pl(arg1, (float8) arg2));
4040 : }
4041 :
4042 : Datum
4043 0 : float84mi(PG_FUNCTION_ARGS)
4044 : {
4045 0 : float8 arg1 = PG_GETARG_FLOAT8(0);
4046 0 : float4 arg2 = PG_GETARG_FLOAT4(1);
4047 :
4048 0 : PG_RETURN_FLOAT8(float8_mi(arg1, (float8) arg2));
4049 : }
4050 :
4051 : Datum
4052 0 : float84mul(PG_FUNCTION_ARGS)
4053 : {
4054 0 : float8 arg1 = PG_GETARG_FLOAT8(0);
4055 0 : float4 arg2 = PG_GETARG_FLOAT4(1);
4056 :
4057 0 : PG_RETURN_FLOAT8(float8_mul(arg1, (float8) arg2));
4058 : }
4059 :
4060 : Datum
4061 4 : float84div(PG_FUNCTION_ARGS)
4062 : {
4063 4 : float8 arg1 = PG_GETARG_FLOAT8(0);
4064 4 : float4 arg2 = PG_GETARG_FLOAT4(1);
4065 :
4066 4 : PG_RETURN_FLOAT8(float8_div(arg1, (float8) arg2));
4067 : }
4068 :
4069 : /*
4070 : * ====================
4071 : * COMPARISON OPERATORS
4072 : * ====================
4073 : */
4074 :
4075 : /*
4076 : * float48{eq,ne,lt,le,gt,ge} - float4/float8 comparison operations
4077 : */
4078 : Datum
4079 1812 : float48eq(PG_FUNCTION_ARGS)
4080 : {
4081 1812 : float4 arg1 = PG_GETARG_FLOAT4(0);
4082 1812 : float8 arg2 = PG_GETARG_FLOAT8(1);
4083 :
4084 1812 : PG_RETURN_BOOL(float8_eq((float8) arg1, arg2));
4085 : }
4086 :
4087 : Datum
4088 13534 : float48ne(PG_FUNCTION_ARGS)
4089 : {
4090 13534 : float4 arg1 = PG_GETARG_FLOAT4(0);
4091 13534 : float8 arg2 = PG_GETARG_FLOAT8(1);
4092 :
4093 13534 : PG_RETURN_BOOL(float8_ne((float8) arg1, arg2));
4094 : }
4095 :
4096 : Datum
4097 2663 : float48lt(PG_FUNCTION_ARGS)
4098 : {
4099 2663 : float4 arg1 = PG_GETARG_FLOAT4(0);
4100 2663 : float8 arg2 = PG_GETARG_FLOAT8(1);
4101 :
4102 2663 : PG_RETURN_BOOL(float8_lt((float8) arg1, arg2));
4103 : }
4104 :
4105 : Datum
4106 17577 : float48le(PG_FUNCTION_ARGS)
4107 : {
4108 17577 : float4 arg1 = PG_GETARG_FLOAT4(0);
4109 17577 : float8 arg2 = PG_GETARG_FLOAT8(1);
4110 :
4111 17577 : PG_RETURN_BOOL(float8_le((float8) arg1, arg2));
4112 : }
4113 :
4114 : Datum
4115 2811 : float48gt(PG_FUNCTION_ARGS)
4116 : {
4117 2811 : float4 arg1 = PG_GETARG_FLOAT4(0);
4118 2811 : float8 arg2 = PG_GETARG_FLOAT8(1);
4119 :
4120 2811 : PG_RETURN_BOOL(float8_gt((float8) arg1, arg2));
4121 : }
4122 :
4123 : Datum
4124 3079 : float48ge(PG_FUNCTION_ARGS)
4125 : {
4126 3079 : float4 arg1 = PG_GETARG_FLOAT4(0);
4127 3079 : float8 arg2 = PG_GETARG_FLOAT8(1);
4128 :
4129 3079 : PG_RETURN_BOOL(float8_ge((float8) arg1, arg2));
4130 : }
4131 :
4132 : /*
4133 : * float84{eq,ne,lt,le,gt,ge} - float8/float4 comparison operations
4134 : */
4135 : Datum
4136 1208 : float84eq(PG_FUNCTION_ARGS)
4137 : {
4138 1208 : float8 arg1 = PG_GETARG_FLOAT8(0);
4139 1208 : float4 arg2 = PG_GETARG_FLOAT4(1);
4140 :
4141 1208 : PG_RETURN_BOOL(float8_eq(arg1, (float8) arg2));
4142 : }
4143 :
4144 : Datum
4145 0 : float84ne(PG_FUNCTION_ARGS)
4146 : {
4147 0 : float8 arg1 = PG_GETARG_FLOAT8(0);
4148 0 : float4 arg2 = PG_GETARG_FLOAT4(1);
4149 :
4150 0 : PG_RETURN_BOOL(float8_ne(arg1, (float8) arg2));
4151 : }
4152 :
4153 : Datum
4154 2132 : float84lt(PG_FUNCTION_ARGS)
4155 : {
4156 2132 : float8 arg1 = PG_GETARG_FLOAT8(0);
4157 2132 : float4 arg2 = PG_GETARG_FLOAT4(1);
4158 :
4159 2132 : PG_RETURN_BOOL(float8_lt(arg1, (float8) arg2));
4160 : }
4161 :
4162 : Datum
4163 2532 : float84le(PG_FUNCTION_ARGS)
4164 : {
4165 2532 : float8 arg1 = PG_GETARG_FLOAT8(0);
4166 2532 : float4 arg2 = PG_GETARG_FLOAT4(1);
4167 :
4168 2532 : PG_RETURN_BOOL(float8_le(arg1, (float8) arg2));
4169 : }
4170 :
4171 : Datum
4172 2132 : float84gt(PG_FUNCTION_ARGS)
4173 : {
4174 2132 : float8 arg1 = PG_GETARG_FLOAT8(0);
4175 2132 : float4 arg2 = PG_GETARG_FLOAT4(1);
4176 :
4177 2132 : PG_RETURN_BOOL(float8_gt(arg1, (float8) arg2));
4178 : }
4179 :
4180 : Datum
4181 2136 : float84ge(PG_FUNCTION_ARGS)
4182 : {
4183 2136 : float8 arg1 = PG_GETARG_FLOAT8(0);
4184 2136 : float4 arg2 = PG_GETARG_FLOAT4(1);
4185 :
4186 2136 : PG_RETURN_BOOL(float8_ge(arg1, (float8) arg2));
4187 : }
4188 :
4189 : /*
4190 : * Implements the float8 version of the width_bucket() function
4191 : * defined by SQL2003. See also width_bucket_numeric().
4192 : *
4193 : * 'bound1' and 'bound2' are the lower and upper bounds of the
4194 : * histogram's range, respectively. 'count' is the number of buckets
4195 : * in the histogram. width_bucket() returns an integer indicating the
4196 : * bucket number that 'operand' belongs to in an equiwidth histogram
4197 : * with the specified characteristics. An operand smaller than the
4198 : * lower bound is assigned to bucket 0. An operand greater than or equal
4199 : * to the upper bound is assigned to an additional bucket (with number
4200 : * count+1). We don't allow the histogram bounds to be NaN or +/- infinity,
4201 : * but we do allow those values for the operand (taking NaN to be larger
4202 : * than any other value, as we do in comparisons).
4203 : */
4204 : Datum
4205 585 : width_bucket_float8(PG_FUNCTION_ARGS)
4206 : {
4207 585 : float8 operand = PG_GETARG_FLOAT8(0);
4208 585 : float8 bound1 = PG_GETARG_FLOAT8(1);
4209 585 : float8 bound2 = PG_GETARG_FLOAT8(2);
4210 585 : int32 count = PG_GETARG_INT32(3);
4211 : int32 result;
4212 :
4213 585 : if (count <= 0)
4214 8 : ereport(ERROR,
4215 : (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
4216 : errmsg("count must be greater than zero")));
4217 :
4218 577 : if (isnan(bound1) || isnan(bound2))
4219 4 : ereport(ERROR,
4220 : (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
4221 : errmsg("lower and upper bounds cannot be NaN")));
4222 :
4223 573 : if (isinf(bound1) || isinf(bound2))
4224 12 : ereport(ERROR,
4225 : (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
4226 : errmsg("lower and upper bounds must be finite")));
4227 :
4228 561 : if (bound1 < bound2)
4229 : {
4230 404 : if (isnan(operand) || operand >= bound2)
4231 : {
4232 86 : if (pg_add_s32_overflow(count, 1, &result))
4233 4 : ereport(ERROR,
4234 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
4235 : errmsg("integer out of range")));
4236 : }
4237 318 : else if (operand < bound1)
4238 77 : result = 0;
4239 : else
4240 : {
4241 241 : if (!isinf(bound2 - bound1))
4242 : {
4243 : /* The quotient is surely in [0,1], so this can't overflow */
4244 229 : result = count * ((operand - bound1) / (bound2 - bound1));
4245 : }
4246 : else
4247 : {
4248 : /*
4249 : * We get here if bound2 - bound1 overflows DBL_MAX. Since
4250 : * both bounds are finite, their difference can't exceed twice
4251 : * DBL_MAX; so we can perform the computation without overflow
4252 : * by dividing all the inputs by 2. That should be exact too,
4253 : * except in the case where a very small operand underflows to
4254 : * zero, which would have negligible impact on the result
4255 : * given such large bounds.
4256 : */
4257 12 : result = count * ((operand / 2 - bound1 / 2) / (bound2 / 2 - bound1 / 2));
4258 : }
4259 : /* The quotient could round to 1.0, which would be a lie */
4260 241 : if (result >= count)
4261 5 : result = count - 1;
4262 : /* Having done that, we can add 1 without fear of overflow */
4263 241 : result++;
4264 : }
4265 : }
4266 157 : else if (bound1 > bound2)
4267 : {
4268 153 : if (isnan(operand) || operand > bound1)
4269 8 : result = 0;
4270 145 : else if (operand <= bound2)
4271 : {
4272 20 : if (pg_add_s32_overflow(count, 1, &result))
4273 4 : ereport(ERROR,
4274 : (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
4275 : errmsg("integer out of range")));
4276 : }
4277 : else
4278 : {
4279 125 : if (!isinf(bound1 - bound2))
4280 113 : result = count * ((bound1 - operand) / (bound1 - bound2));
4281 : else
4282 12 : result = count * ((bound1 / 2 - operand / 2) / (bound1 / 2 - bound2 / 2));
4283 125 : if (result >= count)
4284 5 : result = count - 1;
4285 125 : result++;
4286 : }
4287 : }
4288 : else
4289 : {
4290 4 : ereport(ERROR,
4291 : (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
4292 : errmsg("lower bound cannot equal upper bound")));
4293 : result = 0; /* keep the compiler quiet */
4294 : }
4295 :
4296 549 : PG_RETURN_INT32(result);
4297 : }
|
