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