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

Generated by: LCOV version 1.14