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