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