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