LCOV - code coverage report
Current view: top level - src/backend/utils/adt - float.c (source / functions) Hit Total Coverage
Test: PostgreSQL 18devel Lines: 1108 1359 81.5 %
Date: 2025-04-01 14:15:22 Functions: 141 157 89.8 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14