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