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