granicus.if.org Git - postgresql/blob - src/backend/utils/adt/float.c

   1 /*-------------------------------------------------------------------------
   2  *
   3  * float.c
   4  *        Functions for the built-in floating-point types.
   5  *
   6  * Portions Copyright (c) 1996-2014, 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 "libpq/pqformat.h"
  24 #include "utils/array.h"
  25 #include "utils/builtins.h"
  26 #include "utils/sortsupport.h"
  27
  28
  29 #ifndef M_PI
  30 /* from my RH5.2 gcc math.h file - thomas 2000-04-03 */
  31 #define M_PI 3.14159265358979323846
  32 #endif
  33
  34 /* Visual C++ etc lacks NAN, and won't accept 0.0/0.0.  NAN definition from
  35  * http://msdn.microsoft.com/library/default.asp?url=/library/en-us/vclang/html/vclrfNotNumberNANItems.asp
  36  */
  37 #if defined(WIN32) && !defined(NAN)
  38 static const uint32 nan[2] = {0xffffffff, 0x7fffffff};
  39
  40 #define NAN (*(const double *) nan)
  41 #endif
  42
  43 /* not sure what the following should be, but better to make it over-sufficient */
  44 #define MAXFLOATWIDTH   64
  45 #define MAXDOUBLEWIDTH  128
  46
  47 /*
  48  * check to see if a float4/8 val has underflowed or overflowed
  49  */
  50 #define CHECKFLOATVAL(val, inf_is_valid, zero_is_valid)                 \
  51 do {                                                                                                                    \
  52         if (isinf(val) && !(inf_is_valid))                                                      \
  53                 ereport(ERROR,                                                                                  \
  54                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),   \
  55                   errmsg("value out of range: overflow")));                             \
  56                                                                                                                                 \
  57         if ((val) == 0.0 && !(zero_is_valid))                                           \
  58                 ereport(ERROR,                                                                                  \
  59                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),   \
  60                  errmsg("value out of range: underflow")));                             \
  61 } while(0)
  62
  63
  64 /* ========== USER I/O ROUTINES ========== */
  65
  66
  67 /* Configurable GUC parameter */
  68 int                     extra_float_digits = 0;         /* Added to DBL_DIG or FLT_DIG */
  69
  70
  71 static int      float4_cmp_internal(float4 a, float4 b);
  72 static int      float8_cmp_internal(float8 a, float8 b);
  73
  74 #ifndef HAVE_CBRT
  75 /*
  76  * Some machines (in particular, some versions of AIX) have an extern
  77  * declaration for cbrt() in <math.h> but fail to provide the actual
  78  * function, which causes configure to not set HAVE_CBRT.  Furthermore,
  79  * their compilers spit up at the mismatch between extern declaration
  80  * and static definition.  We work around that here by the expedient
  81  * of a #define to make the actual name of the static function different.
  82  */
  83 #define cbrt my_cbrt
  84 static double cbrt(double x);
  85 #endif   /* HAVE_CBRT */
  86
  87
  88 /*
  89  * Routines to provide reasonably platform-independent handling of
  90  * infinity and NaN.  We assume that isinf() and isnan() are available
  91  * and work per spec.  (On some platforms, we have to supply our own;
  92  * see src/port.)  However, generating an Infinity or NaN in the first
  93  * place is less well standardized; pre-C99 systems tend not to have C99's
  94  * INFINITY and NAN macros.  We centralize our workarounds for this here.
  95  */
  96
  97 double
  98 get_float8_infinity(void)
  99 {
 100 #ifdef INFINITY
 101         /* C99 standard way */
 102         return (double) INFINITY;
 103 #else
 104
 105         /*
 106          * On some platforms, HUGE_VAL is an infinity, elsewhere it's just the
 107          * largest normal double.  We assume forcing an overflow will get us a
 108          * true infinity.
 109          */
 110         return (double) (HUGE_VAL * HUGE_VAL);
 111 #endif
 112 }
 113
 114 /*
 115 * The funny placements of the two #pragmas is necessary because of a
 116 * long lived bug in the Microsoft compilers.
 117 * See http://support.microsoft.com/kb/120968/en-us for details
 118 */
 119 #if (_MSC_VER >= 1800)
 120 #pragma warning(disable:4756)
 121 #endif
 122 float
 123 get_float4_infinity(void)
 124 {
 125 #ifdef INFINITY
 126         /* C99 standard way */
 127         return (float) INFINITY;
 128 #else
 129 #if (_MSC_VER >= 1800)
 130 #pragma warning(default:4756)
 131 #endif
 132
 133         /*
 134          * On some platforms, HUGE_VAL is an infinity, elsewhere it's just the
 135          * largest normal double.  We assume forcing an overflow will get us a
 136          * true infinity.
 137          */
 138         return (float) (HUGE_VAL * HUGE_VAL);
 139 #endif
 140 }
 141
 142 double
 143 get_float8_nan(void)
 144 {
 145         /* (double) NAN doesn't work on some NetBSD/MIPS releases */
 146 #if defined(NAN) && !(defined(__NetBSD__) && defined(__mips__))
 147         /* C99 standard way */
 148         return (double) NAN;
 149 #else
 150         /* Assume we can get a NAN via zero divide */
 151         return (double) (0.0 / 0.0);
 152 #endif
 153 }
 154
 155 float
 156 get_float4_nan(void)
 157 {
 158 #ifdef NAN
 159         /* C99 standard way */
 160         return (float) NAN;
 161 #else
 162         /* Assume we can get a NAN via zero divide */
 163         return (float) (0.0 / 0.0);
 164 #endif
 165 }
 166
 167
 168 /*
 169  * Returns -1 if 'val' represents negative infinity, 1 if 'val'
 170  * represents (positive) infinity, and 0 otherwise. On some platforms,
 171  * this is equivalent to the isinf() macro, but not everywhere: C99
 172  * does not specify that isinf() needs to distinguish between positive
 173  * and negative infinity.
 174  */
 175 int
 176 is_infinite(double val)
 177 {
 178         int                     inf = isinf(val);
 179
 180         if (inf == 0)
 181                 return 0;
 182         else if (val > 0)
 183                 return 1;
 184         else
 185                 return -1;
 186 }
 187
 188
 189 /*
 190  *              float4in                - converts "num" to float4
 191  */
 192 Datum
 193 float4in(PG_FUNCTION_ARGS)
 194 {
 195         char       *num = PG_GETARG_CSTRING(0);
 196         char       *orig_num;
 197         double          val;
 198         char       *endptr;
 199
 200         /*
 201          * endptr points to the first character _after_ the sequence we recognized
 202          * as a valid floating point number. orig_num points to the original input
 203          * string.
 204          */
 205         orig_num = num;
 206
 207         /* skip leading whitespace */
 208         while (*num != '\0' && isspace((unsigned char) *num))
 209                 num++;
 210
 211         /*
 212          * Check for an empty-string input to begin with, to avoid the vagaries of
 213          * strtod() on different platforms.
 214          */
 215         if (*num == '\0')
 216                 ereport(ERROR,
 217                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
 218                                  errmsg("invalid input syntax for type real: \"%s\"",
 219                                                 orig_num)));
 220
 221         errno = 0;
 222         val = strtod(num, &endptr);
 223
 224         /* did we not see anything that looks like a double? */
 225         if (endptr == num || errno != 0)
 226         {
 227                 int                     save_errno = errno;
 228
 229                 /*
 230                  * C99 requires that strtod() accept NaN, [+-]Infinity, and [+-]Inf,
 231                  * but not all platforms support all of these (and some accept them
 232                  * but set ERANGE anyway...)  Therefore, we check for these inputs
 233                  * ourselves if strtod() fails.
 234                  *
 235                  * Note: C99 also requires hexadecimal input as well as some extended
 236                  * forms of NaN, but we consider these forms unportable and don't try
 237                  * to support them.  You can use 'em if your strtod() takes 'em.
 238                  */
 239                 if (pg_strncasecmp(num, "NaN", 3) == 0)
 240                 {
 241                         val = get_float4_nan();
 242                         endptr = num + 3;
 243                 }
 244                 else if (pg_strncasecmp(num, "Infinity", 8) == 0)
 245                 {
 246                         val = get_float4_infinity();
 247                         endptr = num + 8;
 248                 }
 249                 else if (pg_strncasecmp(num, "+Infinity", 9) == 0)
 250                 {
 251                         val = get_float4_infinity();
 252                         endptr = num + 9;
 253                 }
 254                 else if (pg_strncasecmp(num, "-Infinity", 9) == 0)
 255                 {
 256                         val = -get_float4_infinity();
 257                         endptr = num + 9;
 258                 }
 259                 else if (pg_strncasecmp(num, "inf", 3) == 0)
 260                 {
 261                         val = get_float4_infinity();
 262                         endptr = num + 3;
 263                 }
 264                 else if (pg_strncasecmp(num, "+inf", 4) == 0)
 265                 {
 266                         val = get_float4_infinity();
 267                         endptr = num + 4;
 268                 }
 269                 else if (pg_strncasecmp(num, "-inf", 4) == 0)
 270                 {
 271                         val = -get_float4_infinity();
 272                         endptr = num + 4;
 273                 }
 274                 else if (save_errno == ERANGE)
 275                 {
 276                         /*
 277                          * Some platforms return ERANGE for denormalized numbers (those
 278                          * that are not zero, but are too close to zero to have full
 279                          * precision).  We'd prefer not to throw error for that, so try to
 280                          * detect whether it's a "real" out-of-range condition by checking
 281                          * to see if the result is zero or huge.
 282                          */
 283                         if (val == 0.0 || val >= HUGE_VAL || val <= -HUGE_VAL)
 284                                 ereport(ERROR,
 285                                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
 286                                                  errmsg("\"%s\" is out of range for type real",
 287                                                                 orig_num)));
 288                 }
 289                 else
 290                         ereport(ERROR,
 291                                         (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
 292                                          errmsg("invalid input syntax for type real: \"%s\"",
 293                                                         orig_num)));
 294         }
 295 #ifdef HAVE_BUGGY_SOLARIS_STRTOD
 296         else
 297         {
 298                 /*
 299                  * Many versions of Solaris have a bug wherein strtod sets endptr to
 300                  * point one byte beyond the end of the string when given "inf" or
 301                  * "infinity".
 302                  */
 303                 if (endptr != num && endptr[-1] == '\0')
 304                         endptr--;
 305         }
 306 #endif   /* HAVE_BUGGY_SOLARIS_STRTOD */
 307
 308         /* skip trailing whitespace */
 309         while (*endptr != '\0' && isspace((unsigned char) *endptr))
 310                 endptr++;
 311
 312         /* if there is any junk left at the end of the string, bail out */
 313         if (*endptr != '\0')
 314                 ereport(ERROR,
 315                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
 316                                  errmsg("invalid input syntax for type real: \"%s\"",
 317                                                 orig_num)));
 318
 319         /*
 320          * if we get here, we have a legal double, still need to check to see if
 321          * it's a legal float4
 322          */
 323         CHECKFLOATVAL((float4) val, isinf(val), val == 0);
 324
 325         PG_RETURN_FLOAT4((float4) val);
 326 }
 327
 328 /*
 329  *              float4out               - converts a float4 number to a string
 330  *                                                using a standard output format
 331  */
 332 Datum
 333 float4out(PG_FUNCTION_ARGS)
 334 {
 335         float4          num = PG_GETARG_FLOAT4(0);
 336         char       *ascii = (char *) palloc(MAXFLOATWIDTH + 1);
 337
 338         if (isnan(num))
 339                 PG_RETURN_CSTRING(strcpy(ascii, "NaN"));
 340
 341         switch (is_infinite(num))
 342         {
 343                 case 1:
 344                         strcpy(ascii, "Infinity");
 345                         break;
 346                 case -1:
 347                         strcpy(ascii, "-Infinity");
 348                         break;
 349                 default:
 350                         {
 351                                 int                     ndig = FLT_DIG + extra_float_digits;
 352
 353                                 if (ndig < 1)
 354                                         ndig = 1;
 355
 356                                 snprintf(ascii, MAXFLOATWIDTH + 1, "%.*g", ndig, num);
 357                         }
 358         }
 359
 360         PG_RETURN_CSTRING(ascii);
 361 }
 362
 363 /*
 364  *              float4recv                      - converts external binary format to float4
 365  */
 366 Datum
 367 float4recv(PG_FUNCTION_ARGS)
 368 {
 369         StringInfo      buf = (StringInfo) PG_GETARG_POINTER(0);
 370
 371         PG_RETURN_FLOAT4(pq_getmsgfloat4(buf));
 372 }
 373
 374 /*
 375  *              float4send                      - converts float4 to binary format
 376  */
 377 Datum
 378 float4send(PG_FUNCTION_ARGS)
 379 {
 380         float4          num = PG_GETARG_FLOAT4(0);
 381         StringInfoData buf;
 382
 383         pq_begintypsend(&buf);
 384         pq_sendfloat4(&buf, num);
 385         PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
 386 }
 387
 388 /*
 389  *              float8in                - converts "num" to float8
 390  */
 391 Datum
 392 float8in(PG_FUNCTION_ARGS)
 393 {
 394         char       *num = PG_GETARG_CSTRING(0);
 395         char       *orig_num;
 396         double          val;
 397         char       *endptr;
 398
 399         /*
 400          * endptr points to the first character _after_ the sequence we recognized
 401          * as a valid floating point number. orig_num points to the original input
 402          * string.
 403          */
 404         orig_num = num;
 405
 406         /* skip leading whitespace */
 407         while (*num != '\0' && isspace((unsigned char) *num))
 408                 num++;
 409
 410         /*
 411          * Check for an empty-string input to begin with, to avoid the vagaries of
 412          * strtod() on different platforms.
 413          */
 414         if (*num == '\0')
 415                 ereport(ERROR,
 416                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
 417                          errmsg("invalid input syntax for type double precision: \"%s\"",
 418                                         orig_num)));
 419
 420         errno = 0;
 421         val = strtod(num, &endptr);
 422
 423         /* did we not see anything that looks like a double? */
 424         if (endptr == num || errno != 0)
 425         {
 426                 int                     save_errno = errno;
 427
 428                 /*
 429                  * C99 requires that strtod() accept NaN, [+-]Infinity, and [+-]Inf,
 430                  * but not all platforms support all of these (and some accept them
 431                  * but set ERANGE anyway...)  Therefore, we check for these inputs
 432                  * ourselves if strtod() fails.
 433                  *
 434                  * Note: C99 also requires hexadecimal input as well as some extended
 435                  * forms of NaN, but we consider these forms unportable and don't try
 436                  * to support them.  You can use 'em if your strtod() takes 'em.
 437                  */
 438                 if (pg_strncasecmp(num, "NaN", 3) == 0)
 439                 {
 440                         val = get_float8_nan();
 441                         endptr = num + 3;
 442                 }
 443                 else if (pg_strncasecmp(num, "Infinity", 8) == 0)
 444                 {
 445                         val = get_float8_infinity();
 446                         endptr = num + 8;
 447                 }
 448                 else if (pg_strncasecmp(num, "+Infinity", 9) == 0)
 449                 {
 450                         val = get_float8_infinity();
 451                         endptr = num + 9;
 452                 }
 453                 else if (pg_strncasecmp(num, "-Infinity", 9) == 0)
 454                 {
 455                         val = -get_float8_infinity();
 456                         endptr = num + 9;
 457                 }
 458                 else if (pg_strncasecmp(num, "inf", 3) == 0)
 459                 {
 460                         val = get_float8_infinity();
 461                         endptr = num + 3;
 462                 }
 463                 else if (pg_strncasecmp(num, "+inf", 4) == 0)
 464                 {
 465                         val = get_float8_infinity();
 466                         endptr = num + 4;
 467                 }
 468                 else if (pg_strncasecmp(num, "-inf", 4) == 0)
 469                 {
 470                         val = -get_float8_infinity();
 471                         endptr = num + 4;
 472                 }
 473                 else if (save_errno == ERANGE)
 474                 {
 475                         /*
 476                          * Some platforms return ERANGE for denormalized numbers (those
 477                          * that are not zero, but are too close to zero to have full
 478                          * precision).  We'd prefer not to throw error for that, so try to
 479                          * detect whether it's a "real" out-of-range condition by checking
 480                          * to see if the result is zero or huge.
 481                          */
 482                         if (val == 0.0 || val >= HUGE_VAL || val <= -HUGE_VAL)
 483                                 ereport(ERROR,
 484                                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
 485                                    errmsg("\"%s\" is out of range for type double precision",
 486                                                   orig_num)));
 487                 }
 488                 else
 489                         ereport(ERROR,
 490                                         (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
 491                          errmsg("invalid input syntax for type double precision: \"%s\"",
 492                                         orig_num)));
 493         }
 494 #ifdef HAVE_BUGGY_SOLARIS_STRTOD
 495         else
 496         {
 497                 /*
 498                  * Many versions of Solaris have a bug wherein strtod sets endptr to
 499                  * point one byte beyond the end of the string when given "inf" or
 500                  * "infinity".
 501                  */
 502                 if (endptr != num && endptr[-1] == '\0')
 503                         endptr--;
 504         }
 505 #endif   /* HAVE_BUGGY_SOLARIS_STRTOD */
 506
 507         /* skip trailing whitespace */
 508         while (*endptr != '\0' && isspace((unsigned char) *endptr))
 509                 endptr++;
 510
 511         /* if there is any junk left at the end of the string, bail out */
 512         if (*endptr != '\0')
 513                 ereport(ERROR,
 514                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
 515                          errmsg("invalid input syntax for type double precision: \"%s\"",
 516                                         orig_num)));
 517
 518         CHECKFLOATVAL(val, true, true);
 519
 520         PG_RETURN_FLOAT8(val);
 521 }
 522
 523 /*
 524  *              float8out               - converts float8 number to a string
 525  *                                                using a standard output format
 526  */
 527 Datum
 528 float8out(PG_FUNCTION_ARGS)
 529 {
 530         float8          num = PG_GETARG_FLOAT8(0);
 531         char       *ascii = (char *) palloc(MAXDOUBLEWIDTH + 1);
 532
 533         if (isnan(num))
 534                 PG_RETURN_CSTRING(strcpy(ascii, "NaN"));
 535
 536         switch (is_infinite(num))
 537         {
 538                 case 1:
 539                         strcpy(ascii, "Infinity");
 540                         break;
 541                 case -1:
 542                         strcpy(ascii, "-Infinity");
 543                         break;
 544                 default:
 545                         {
 546                                 int                     ndig = DBL_DIG + extra_float_digits;
 547
 548                                 if (ndig < 1)
 549                                         ndig = 1;
 550
 551                                 snprintf(ascii, MAXDOUBLEWIDTH + 1, "%.*g", ndig, num);
 552                         }
 553         }
 554
 555         PG_RETURN_CSTRING(ascii);
 556 }
 557
 558 /*
 559  *              float8recv                      - converts external binary format to float8
 560  */
 561 Datum
 562 float8recv(PG_FUNCTION_ARGS)
 563 {
 564         StringInfo      buf = (StringInfo) PG_GETARG_POINTER(0);
 565
 566         PG_RETURN_FLOAT8(pq_getmsgfloat8(buf));
 567 }
 568
 569 /*
 570  *              float8send                      - converts float8 to binary format
 571  */
 572 Datum
 573 float8send(PG_FUNCTION_ARGS)
 574 {
 575         float8          num = PG_GETARG_FLOAT8(0);
 576         StringInfoData buf;
 577
 578         pq_begintypsend(&buf);
 579         pq_sendfloat8(&buf, num);
 580         PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
 581 }
 582
 583
 584 /* ========== PUBLIC ROUTINES ========== */
 585
 586
 587 /*
 588  *              ======================
 589  *              FLOAT4 BASE OPERATIONS
 590  *              ======================
 591  */
 592
 593 /*
 594  *              float4abs               - returns |arg1| (absolute value)
 595  */
 596 Datum
 597 float4abs(PG_FUNCTION_ARGS)
 598 {
 599         float4          arg1 = PG_GETARG_FLOAT4(0);
 600
 601         PG_RETURN_FLOAT4((float4) fabs(arg1));
 602 }
 603
 604 /*
 605  *              float4um                - returns -arg1 (unary minus)
 606  */
 607 Datum
 608 float4um(PG_FUNCTION_ARGS)
 609 {
 610         float4          arg1 = PG_GETARG_FLOAT4(0);
 611         float4          result;
 612
 613         result = -arg1;
 614         PG_RETURN_FLOAT4(result);
 615 }
 616
 617 Datum
 618 float4up(PG_FUNCTION_ARGS)
 619 {
 620         float4          arg = PG_GETARG_FLOAT4(0);
 621
 622         PG_RETURN_FLOAT4(arg);
 623 }
 624
 625 Datum
 626 float4larger(PG_FUNCTION_ARGS)
 627 {
 628         float4          arg1 = PG_GETARG_FLOAT4(0);
 629         float4          arg2 = PG_GETARG_FLOAT4(1);
 630         float4          result;
 631
 632         if (float4_cmp_internal(arg1, arg2) > 0)
 633                 result = arg1;
 634         else
 635                 result = arg2;
 636         PG_RETURN_FLOAT4(result);
 637 }
 638
 639 Datum
 640 float4smaller(PG_FUNCTION_ARGS)
 641 {
 642         float4          arg1 = PG_GETARG_FLOAT4(0);
 643         float4          arg2 = PG_GETARG_FLOAT4(1);
 644         float4          result;
 645
 646         if (float4_cmp_internal(arg1, arg2) < 0)
 647                 result = arg1;
 648         else
 649                 result = arg2;
 650         PG_RETURN_FLOAT4(result);
 651 }
 652
 653 /*
 654  *              ======================
 655  *              FLOAT8 BASE OPERATIONS
 656  *              ======================
 657  */
 658
 659 /*
 660  *              float8abs               - returns |arg1| (absolute value)
 661  */
 662 Datum
 663 float8abs(PG_FUNCTION_ARGS)
 664 {
 665         float8          arg1 = PG_GETARG_FLOAT8(0);
 666
 667         PG_RETURN_FLOAT8(fabs(arg1));
 668 }
 669
 670
 671 /*
 672  *              float8um                - returns -arg1 (unary minus)
 673  */
 674 Datum
 675 float8um(PG_FUNCTION_ARGS)
 676 {
 677         float8          arg1 = PG_GETARG_FLOAT8(0);
 678         float8          result;
 679
 680         result = -arg1;
 681         PG_RETURN_FLOAT8(result);
 682 }
 683
 684 Datum
 685 float8up(PG_FUNCTION_ARGS)
 686 {
 687         float8          arg = PG_GETARG_FLOAT8(0);
 688
 689         PG_RETURN_FLOAT8(arg);
 690 }
 691
 692 Datum
 693 float8larger(PG_FUNCTION_ARGS)
 694 {
 695         float8          arg1 = PG_GETARG_FLOAT8(0);
 696         float8          arg2 = PG_GETARG_FLOAT8(1);
 697         float8          result;
 698
 699         if (float8_cmp_internal(arg1, arg2) > 0)
 700                 result = arg1;
 701         else
 702                 result = arg2;
 703         PG_RETURN_FLOAT8(result);
 704 }
 705
 706 Datum
 707 float8smaller(PG_FUNCTION_ARGS)
 708 {
 709         float8          arg1 = PG_GETARG_FLOAT8(0);
 710         float8          arg2 = PG_GETARG_FLOAT8(1);
 711         float8          result;
 712
 713         if (float8_cmp_internal(arg1, arg2) < 0)
 714                 result = arg1;
 715         else
 716                 result = arg2;
 717         PG_RETURN_FLOAT8(result);
 718 }
 719
 720
 721 /*
 722  *              ====================
 723  *              ARITHMETIC OPERATORS
 724  *              ====================
 725  */
 726
 727 /*
 728  *              float4pl                - returns arg1 + arg2
 729  *              float4mi                - returns arg1 - arg2
 730  *              float4mul               - returns arg1 * arg2
 731  *              float4div               - returns arg1 / arg2
 732  */
 733 Datum
 734 float4pl(PG_FUNCTION_ARGS)
 735 {
 736         float4          arg1 = PG_GETARG_FLOAT4(0);
 737         float4          arg2 = PG_GETARG_FLOAT4(1);
 738         float4          result;
 739
 740         result = arg1 + arg2;
 741
 742         /*
 743          * There isn't any way to check for underflow of addition/subtraction
 744          * because numbers near the underflow value have already been rounded to
 745          * the point where we can't detect that the two values were originally
 746          * different, e.g. on x86, '1e-45'::float4 == '2e-45'::float4 ==
 747          * 1.4013e-45.
 748          */
 749         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
 750         PG_RETURN_FLOAT4(result);
 751 }
 752
 753 Datum
 754 float4mi(PG_FUNCTION_ARGS)
 755 {
 756         float4          arg1 = PG_GETARG_FLOAT4(0);
 757         float4          arg2 = PG_GETARG_FLOAT4(1);
 758         float4          result;
 759
 760         result = arg1 - arg2;
 761         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
 762         PG_RETURN_FLOAT4(result);
 763 }
 764
 765 Datum
 766 float4mul(PG_FUNCTION_ARGS)
 767 {
 768         float4          arg1 = PG_GETARG_FLOAT4(0);
 769         float4          arg2 = PG_GETARG_FLOAT4(1);
 770         float4          result;
 771
 772         result = arg1 * arg2;
 773         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
 774                                   arg1 == 0 || arg2 == 0);
 775         PG_RETURN_FLOAT4(result);
 776 }
 777
 778 Datum
 779 float4div(PG_FUNCTION_ARGS)
 780 {
 781         float4          arg1 = PG_GETARG_FLOAT4(0);
 782         float4          arg2 = PG_GETARG_FLOAT4(1);
 783         float4          result;
 784
 785         if (arg2 == 0.0)
 786                 ereport(ERROR,
 787                                 (errcode(ERRCODE_DIVISION_BY_ZERO),
 788                                  errmsg("division by zero")));
 789
 790         result = arg1 / arg2;
 791
 792         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
 793         PG_RETURN_FLOAT4(result);
 794 }
 795
 796 /*
 797  *              float8pl                - returns arg1 + arg2
 798  *              float8mi                - returns arg1 - arg2
 799  *              float8mul               - returns arg1 * arg2
 800  *              float8div               - returns arg1 / arg2
 801  */
 802 Datum
 803 float8pl(PG_FUNCTION_ARGS)
 804 {
 805         float8          arg1 = PG_GETARG_FLOAT8(0);
 806         float8          arg2 = PG_GETARG_FLOAT8(1);
 807         float8          result;
 808
 809         result = arg1 + arg2;
 810
 811         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
 812         PG_RETURN_FLOAT8(result);
 813 }
 814
 815 Datum
 816 float8mi(PG_FUNCTION_ARGS)
 817 {
 818         float8          arg1 = PG_GETARG_FLOAT8(0);
 819         float8          arg2 = PG_GETARG_FLOAT8(1);
 820         float8          result;
 821
 822         result = arg1 - arg2;
 823
 824         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
 825         PG_RETURN_FLOAT8(result);
 826 }
 827
 828 Datum
 829 float8mul(PG_FUNCTION_ARGS)
 830 {
 831         float8          arg1 = PG_GETARG_FLOAT8(0);
 832         float8          arg2 = PG_GETARG_FLOAT8(1);
 833         float8          result;
 834
 835         result = arg1 * arg2;
 836
 837         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
 838                                   arg1 == 0 || arg2 == 0);
 839         PG_RETURN_FLOAT8(result);
 840 }
 841
 842 Datum
 843 float8div(PG_FUNCTION_ARGS)
 844 {
 845         float8          arg1 = PG_GETARG_FLOAT8(0);
 846         float8          arg2 = PG_GETARG_FLOAT8(1);
 847         float8          result;
 848
 849         if (arg2 == 0.0)
 850                 ereport(ERROR,
 851                                 (errcode(ERRCODE_DIVISION_BY_ZERO),
 852                                  errmsg("division by zero")));
 853
 854         result = arg1 / arg2;
 855
 856         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
 857         PG_RETURN_FLOAT8(result);
 858 }
 859
 860
 861 /*
 862  *              ====================
 863  *              COMPARISON OPERATORS
 864  *              ====================
 865  */
 866
 867 /*
 868  *              float4{eq,ne,lt,le,gt,ge}               - float4/float4 comparison operations
 869  */
 870 static int
 871 float4_cmp_internal(float4 a, float4 b)
 872 {
 873         /*
 874          * We consider all NANs to be equal and larger than any non-NAN. This is
 875          * somewhat arbitrary; the important thing is to have a consistent sort
 876          * order.
 877          */
 878         if (isnan(a))
 879         {
 880                 if (isnan(b))
 881                         return 0;                       /* NAN = NAN */
 882                 else
 883                         return 1;                       /* NAN > non-NAN */
 884         }
 885         else if (isnan(b))
 886         {
 887                 return -1;                              /* non-NAN < NAN */
 888         }
 889         else
 890         {
 891                 if (a > b)
 892                         return 1;
 893                 else if (a < b)
 894                         return -1;
 895                 else
 896                         return 0;
 897         }
 898 }
 899
 900 Datum
 901 float4eq(PG_FUNCTION_ARGS)
 902 {
 903         float4          arg1 = PG_GETARG_FLOAT4(0);
 904         float4          arg2 = PG_GETARG_FLOAT4(1);
 905
 906         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) == 0);
 907 }
 908
 909 Datum
 910 float4ne(PG_FUNCTION_ARGS)
 911 {
 912         float4          arg1 = PG_GETARG_FLOAT4(0);
 913         float4          arg2 = PG_GETARG_FLOAT4(1);
 914
 915         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) != 0);
 916 }
 917
 918 Datum
 919 float4lt(PG_FUNCTION_ARGS)
 920 {
 921         float4          arg1 = PG_GETARG_FLOAT4(0);
 922         float4          arg2 = PG_GETARG_FLOAT4(1);
 923
 924         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) < 0);
 925 }
 926
 927 Datum
 928 float4le(PG_FUNCTION_ARGS)
 929 {
 930         float4          arg1 = PG_GETARG_FLOAT4(0);
 931         float4          arg2 = PG_GETARG_FLOAT4(1);
 932
 933         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) <= 0);
 934 }
 935
 936 Datum
 937 float4gt(PG_FUNCTION_ARGS)
 938 {
 939         float4          arg1 = PG_GETARG_FLOAT4(0);
 940         float4          arg2 = PG_GETARG_FLOAT4(1);
 941
 942         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) > 0);
 943 }
 944
 945 Datum
 946 float4ge(PG_FUNCTION_ARGS)
 947 {
 948         float4          arg1 = PG_GETARG_FLOAT4(0);
 949         float4          arg2 = PG_GETARG_FLOAT4(1);
 950
 951         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) >= 0);
 952 }
 953
 954 Datum
 955 btfloat4cmp(PG_FUNCTION_ARGS)
 956 {
 957         float4          arg1 = PG_GETARG_FLOAT4(0);
 958         float4          arg2 = PG_GETARG_FLOAT4(1);
 959
 960         PG_RETURN_INT32(float4_cmp_internal(arg1, arg2));
 961 }
 962
 963 static int
 964 btfloat4fastcmp(Datum x, Datum y, SortSupport ssup)
 965 {
 966         float4          arg1 = DatumGetFloat4(x);
 967         float4          arg2 = DatumGetFloat4(y);
 968
 969         return float4_cmp_internal(arg1, arg2);
 970 }
 971
 972 Datum
 973 btfloat4sortsupport(PG_FUNCTION_ARGS)
 974 {
 975         SortSupport ssup = (SortSupport) PG_GETARG_POINTER(0);
 976
 977         ssup->comparator = btfloat4fastcmp;
 978         PG_RETURN_VOID();
 979 }
 980
 981 /*
 982  *              float8{eq,ne,lt,le,gt,ge}               - float8/float8 comparison operations
 983  */
 984 static int
 985 float8_cmp_internal(float8 a, float8 b)
 986 {
 987         /*
 988          * We consider all NANs to be equal and larger than any non-NAN. This is
 989          * somewhat arbitrary; the important thing is to have a consistent sort
 990          * order.
 991          */
 992         if (isnan(a))
 993         {
 994                 if (isnan(b))
 995                         return 0;                       /* NAN = NAN */
 996                 else
 997                         return 1;                       /* NAN > non-NAN */
 998         }
 999         else if (isnan(b))
1000         {
1001                 return -1;                              /* non-NAN < NAN */
1002         }
1003         else
1004         {
1005                 if (a > b)
1006                         return 1;
1007                 else if (a < b)
1008                         return -1;
1009                 else
1010                         return 0;
1011         }
1012 }
1013
1014 Datum
1015 float8eq(PG_FUNCTION_ARGS)
1016 {
1017         float8          arg1 = PG_GETARG_FLOAT8(0);
1018         float8          arg2 = PG_GETARG_FLOAT8(1);
1019
1020         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
1021 }
1022
1023 Datum
1024 float8ne(PG_FUNCTION_ARGS)
1025 {
1026         float8          arg1 = PG_GETARG_FLOAT8(0);
1027         float8          arg2 = PG_GETARG_FLOAT8(1);
1028
1029         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
1030 }
1031
1032 Datum
1033 float8lt(PG_FUNCTION_ARGS)
1034 {
1035         float8          arg1 = PG_GETARG_FLOAT8(0);
1036         float8          arg2 = PG_GETARG_FLOAT8(1);
1037
1038         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
1039 }
1040
1041 Datum
1042 float8le(PG_FUNCTION_ARGS)
1043 {
1044         float8          arg1 = PG_GETARG_FLOAT8(0);
1045         float8          arg2 = PG_GETARG_FLOAT8(1);
1046
1047         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
1048 }
1049
1050 Datum
1051 float8gt(PG_FUNCTION_ARGS)
1052 {
1053         float8          arg1 = PG_GETARG_FLOAT8(0);
1054         float8          arg2 = PG_GETARG_FLOAT8(1);
1055
1056         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
1057 }
1058
1059 Datum
1060 float8ge(PG_FUNCTION_ARGS)
1061 {
1062         float8          arg1 = PG_GETARG_FLOAT8(0);
1063         float8          arg2 = PG_GETARG_FLOAT8(1);
1064
1065         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
1066 }
1067
1068 Datum
1069 btfloat8cmp(PG_FUNCTION_ARGS)
1070 {
1071         float8          arg1 = PG_GETARG_FLOAT8(0);
1072         float8          arg2 = PG_GETARG_FLOAT8(1);
1073
1074         PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
1075 }
1076
1077 static int
1078 btfloat8fastcmp(Datum x, Datum y, SortSupport ssup)
1079 {
1080         float8          arg1 = DatumGetFloat8(x);
1081         float8          arg2 = DatumGetFloat8(y);
1082
1083         return float8_cmp_internal(arg1, arg2);
1084 }
1085
1086 Datum
1087 btfloat8sortsupport(PG_FUNCTION_ARGS)
1088 {
1089         SortSupport ssup = (SortSupport) PG_GETARG_POINTER(0);
1090
1091         ssup->comparator = btfloat8fastcmp;
1092         PG_RETURN_VOID();
1093 }
1094
1095 Datum
1096 btfloat48cmp(PG_FUNCTION_ARGS)
1097 {
1098         float4          arg1 = PG_GETARG_FLOAT4(0);
1099         float8          arg2 = PG_GETARG_FLOAT8(1);
1100
1101         /* widen float4 to float8 and then compare */
1102         PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
1103 }
1104
1105 Datum
1106 btfloat84cmp(PG_FUNCTION_ARGS)
1107 {
1108         float8          arg1 = PG_GETARG_FLOAT8(0);
1109         float4          arg2 = PG_GETARG_FLOAT4(1);
1110
1111         /* widen float4 to float8 and then compare */
1112         PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
1113 }
1114
1115
1116 /*
1117  *              ===================
1118  *              CONVERSION ROUTINES
1119  *              ===================
1120  */
1121
1122 /*
1123  *              ftod                    - converts a float4 number to a float8 number
1124  */
1125 Datum
1126 ftod(PG_FUNCTION_ARGS)
1127 {
1128         float4          num = PG_GETARG_FLOAT4(0);
1129
1130         PG_RETURN_FLOAT8((float8) num);
1131 }
1132
1133
1134 /*
1135  *              dtof                    - converts a float8 number to a float4 number
1136  */
1137 Datum
1138 dtof(PG_FUNCTION_ARGS)
1139 {
1140         float8          num = PG_GETARG_FLOAT8(0);
1141
1142         CHECKFLOATVAL((float4) num, isinf(num), num == 0);
1143
1144         PG_RETURN_FLOAT4((float4) num);
1145 }
1146
1147
1148 /*
1149  *              dtoi4                   - converts a float8 number to an int4 number
1150  */
1151 Datum
1152 dtoi4(PG_FUNCTION_ARGS)
1153 {
1154         float8          num = PG_GETARG_FLOAT8(0);
1155         int32           result;
1156
1157         /* 'Inf' is handled by INT_MAX */
1158         if (num < INT_MIN || num > INT_MAX || isnan(num))
1159                 ereport(ERROR,
1160                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1161                                  errmsg("integer out of range")));
1162
1163         result = (int32) rint(num);
1164         PG_RETURN_INT32(result);
1165 }
1166
1167
1168 /*
1169  *              dtoi2                   - converts a float8 number to an int2 number
1170  */
1171 Datum
1172 dtoi2(PG_FUNCTION_ARGS)
1173 {
1174         float8          num = PG_GETARG_FLOAT8(0);
1175
1176         if (num < SHRT_MIN || num > SHRT_MAX || isnan(num))
1177                 ereport(ERROR,
1178                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1179                                  errmsg("smallint out of range")));
1180
1181         PG_RETURN_INT16((int16) rint(num));
1182 }
1183
1184
1185 /*
1186  *              i4tod                   - converts an int4 number to a float8 number
1187  */
1188 Datum
1189 i4tod(PG_FUNCTION_ARGS)
1190 {
1191         int32           num = PG_GETARG_INT32(0);
1192
1193         PG_RETURN_FLOAT8((float8) num);
1194 }
1195
1196
1197 /*
1198  *              i2tod                   - converts an int2 number to a float8 number
1199  */
1200 Datum
1201 i2tod(PG_FUNCTION_ARGS)
1202 {
1203         int16           num = PG_GETARG_INT16(0);
1204
1205         PG_RETURN_FLOAT8((float8) num);
1206 }
1207
1208
1209 /*
1210  *              ftoi4                   - converts a float4 number to an int4 number
1211  */
1212 Datum
1213 ftoi4(PG_FUNCTION_ARGS)
1214 {
1215         float4          num = PG_GETARG_FLOAT4(0);
1216
1217         if (num < INT_MIN || num > INT_MAX || isnan(num))
1218                 ereport(ERROR,
1219                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1220                                  errmsg("integer out of range")));
1221
1222         PG_RETURN_INT32((int32) rint(num));
1223 }
1224
1225
1226 /*
1227  *              ftoi2                   - converts a float4 number to an int2 number
1228  */
1229 Datum
1230 ftoi2(PG_FUNCTION_ARGS)
1231 {
1232         float4          num = PG_GETARG_FLOAT4(0);
1233
1234         if (num < SHRT_MIN || num > SHRT_MAX || isnan(num))
1235                 ereport(ERROR,
1236                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1237                                  errmsg("smallint out of range")));
1238
1239         PG_RETURN_INT16((int16) rint(num));
1240 }
1241
1242
1243 /*
1244  *              i4tof                   - converts an int4 number to a float4 number
1245  */
1246 Datum
1247 i4tof(PG_FUNCTION_ARGS)
1248 {
1249         int32           num = PG_GETARG_INT32(0);
1250
1251         PG_RETURN_FLOAT4((float4) num);
1252 }
1253
1254
1255 /*
1256  *              i2tof                   - converts an int2 number to a float4 number
1257  */
1258 Datum
1259 i2tof(PG_FUNCTION_ARGS)
1260 {
1261         int16           num = PG_GETARG_INT16(0);
1262
1263         PG_RETURN_FLOAT4((float4) num);
1264 }
1265
1266
1267 /*
1268  *              =======================
1269  *              RANDOM FLOAT8 OPERATORS
1270  *              =======================
1271  */
1272
1273 /*
1274  *              dround                  - returns       ROUND(arg1)
1275  */
1276 Datum
1277 dround(PG_FUNCTION_ARGS)
1278 {
1279         float8          arg1 = PG_GETARG_FLOAT8(0);
1280
1281         PG_RETURN_FLOAT8(rint(arg1));
1282 }
1283
1284 /*
1285  *              dceil                   - returns the smallest integer greater than or
1286  *                                                equal to the specified float
1287  */
1288 Datum
1289 dceil(PG_FUNCTION_ARGS)
1290 {
1291         float8          arg1 = PG_GETARG_FLOAT8(0);
1292
1293         PG_RETURN_FLOAT8(ceil(arg1));
1294 }
1295
1296 /*
1297  *              dfloor                  - returns the largest integer lesser than or
1298  *                                                equal to the specified float
1299  */
1300 Datum
1301 dfloor(PG_FUNCTION_ARGS)
1302 {
1303         float8          arg1 = PG_GETARG_FLOAT8(0);
1304
1305         PG_RETURN_FLOAT8(floor(arg1));
1306 }
1307
1308 /*
1309  *              dsign                   - returns -1 if the argument is less than 0, 0
1310  *                                                if the argument is equal to 0, and 1 if the
1311  *                                                argument is greater than zero.
1312  */
1313 Datum
1314 dsign(PG_FUNCTION_ARGS)
1315 {
1316         float8          arg1 = PG_GETARG_FLOAT8(0);
1317         float8          result;
1318
1319         if (arg1 > 0)
1320                 result = 1.0;
1321         else if (arg1 < 0)
1322                 result = -1.0;
1323         else
1324                 result = 0.0;
1325
1326         PG_RETURN_FLOAT8(result);
1327 }
1328
1329 /*
1330  *              dtrunc                  - returns truncation-towards-zero of arg1,
1331  *                                                arg1 >= 0 ... the greatest integer less
1332  *                                                                              than or equal to arg1
1333  *                                                arg1 < 0      ... the least integer greater
1334  *                                                                              than or equal to arg1
1335  */
1336 Datum
1337 dtrunc(PG_FUNCTION_ARGS)
1338 {
1339         float8          arg1 = PG_GETARG_FLOAT8(0);
1340         float8          result;
1341
1342         if (arg1 >= 0)
1343                 result = floor(arg1);
1344         else
1345                 result = -floor(-arg1);
1346
1347         PG_RETURN_FLOAT8(result);
1348 }
1349
1350
1351 /*
1352  *              dsqrt                   - returns square root of arg1
1353  */
1354 Datum
1355 dsqrt(PG_FUNCTION_ARGS)
1356 {
1357         float8          arg1 = PG_GETARG_FLOAT8(0);
1358         float8          result;
1359
1360         if (arg1 < 0)
1361                 ereport(ERROR,
1362                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
1363                                  errmsg("cannot take square root of a negative number")));
1364
1365         result = sqrt(arg1);
1366
1367         CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
1368         PG_RETURN_FLOAT8(result);
1369 }
1370
1371
1372 /*
1373  *              dcbrt                   - returns cube root of arg1
1374  */
1375 Datum
1376 dcbrt(PG_FUNCTION_ARGS)
1377 {
1378         float8          arg1 = PG_GETARG_FLOAT8(0);
1379         float8          result;
1380
1381         result = cbrt(arg1);
1382         CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
1383         PG_RETURN_FLOAT8(result);
1384 }
1385
1386
1387 /*
1388  *              dpow                    - returns pow(arg1,arg2)
1389  */
1390 Datum
1391 dpow(PG_FUNCTION_ARGS)
1392 {
1393         float8          arg1 = PG_GETARG_FLOAT8(0);
1394         float8          arg2 = PG_GETARG_FLOAT8(1);
1395         float8          result;
1396
1397         /*
1398          * The SQL spec requires that we emit a particular SQLSTATE error code for
1399          * certain error conditions.  Specifically, we don't return a
1400          * divide-by-zero error code for 0 ^ -1.
1401          */
1402         if (arg1 == 0 && arg2 < 0)
1403                 ereport(ERROR,
1404                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
1405                                  errmsg("zero raised to a negative power is undefined")));
1406         if (arg1 < 0 && floor(arg2) != arg2)
1407                 ereport(ERROR,
1408                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
1409                                  errmsg("a negative number raised to a non-integer power yields a complex result")));
1410
1411         /*
1412          * pow() sets errno only on some platforms, depending on whether it
1413          * follows _IEEE_, _POSIX_, _XOPEN_, or _SVID_, so we try to avoid using
1414          * errno.  However, some platform/CPU combinations return errno == EDOM
1415          * and result == Nan for negative arg1 and very large arg2 (they must be
1416          * using something different from our floor() test to decide it's
1417          * invalid).  Other platforms (HPPA) return errno == ERANGE and a large
1418          * (HUGE_VAL) but finite result to signal overflow.
1419          */
1420         errno = 0;
1421         result = pow(arg1, arg2);
1422         if (errno == EDOM && isnan(result))
1423         {
1424                 if ((fabs(arg1) > 1 && arg2 >= 0) || (fabs(arg1) < 1 && arg2 < 0))
1425                         /* The sign of Inf is not significant in this case. */
1426                         result = get_float8_infinity();
1427                 else if (fabs(arg1) != 1)
1428                         result = 0;
1429                 else
1430                         result = 1;
1431         }
1432         else if (errno == ERANGE && result != 0 && !isinf(result))
1433                 result = get_float8_infinity();
1434
1435         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
1436         PG_RETURN_FLOAT8(result);
1437 }
1438
1439
1440 /*
1441  *              dexp                    - returns the exponential function of arg1
1442  */
1443 Datum
1444 dexp(PG_FUNCTION_ARGS)
1445 {
1446         float8          arg1 = PG_GETARG_FLOAT8(0);
1447         float8          result;
1448
1449         errno = 0;
1450         result = exp(arg1);
1451         if (errno == ERANGE && result != 0 && !isinf(result))
1452                 result = get_float8_infinity();
1453
1454         CHECKFLOATVAL(result, isinf(arg1), false);
1455         PG_RETURN_FLOAT8(result);
1456 }
1457
1458
1459 /*
1460  *              dlog1                   - returns the natural logarithm of arg1
1461  */
1462 Datum
1463 dlog1(PG_FUNCTION_ARGS)
1464 {
1465         float8          arg1 = PG_GETARG_FLOAT8(0);
1466         float8          result;
1467
1468         /*
1469          * Emit particular SQLSTATE error codes for ln(). This is required by the
1470          * SQL standard.
1471          */
1472         if (arg1 == 0.0)
1473                 ereport(ERROR,
1474                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
1475                                  errmsg("cannot take logarithm of zero")));
1476         if (arg1 < 0)
1477                 ereport(ERROR,
1478                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
1479                                  errmsg("cannot take logarithm of a negative number")));
1480
1481         result = log(arg1);
1482
1483         CHECKFLOATVAL(result, isinf(arg1), arg1 == 1);
1484         PG_RETURN_FLOAT8(result);
1485 }
1486
1487
1488 /*
1489  *              dlog10                  - returns the base 10 logarithm of arg1
1490  */
1491 Datum
1492 dlog10(PG_FUNCTION_ARGS)
1493 {
1494         float8          arg1 = PG_GETARG_FLOAT8(0);
1495         float8          result;
1496
1497         /*
1498          * Emit particular SQLSTATE error codes for log(). The SQL spec doesn't
1499          * define log(), but it does define ln(), so it makes sense to emit the
1500          * same error code for an analogous error condition.
1501          */
1502         if (arg1 == 0.0)
1503                 ereport(ERROR,
1504                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
1505                                  errmsg("cannot take logarithm of zero")));
1506         if (arg1 < 0)
1507                 ereport(ERROR,
1508                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
1509                                  errmsg("cannot take logarithm of a negative number")));
1510
1511         result = log10(arg1);
1512
1513         CHECKFLOATVAL(result, isinf(arg1), arg1 == 1);
1514         PG_RETURN_FLOAT8(result);
1515 }
1516
1517
1518 /*
1519  *              dacos                   - returns the arccos of arg1 (radians)
1520  */
1521 Datum
1522 dacos(PG_FUNCTION_ARGS)
1523 {
1524         float8          arg1 = PG_GETARG_FLOAT8(0);
1525         float8          result;
1526
1527         /*
1528          * We use errno here because the trigonometric functions are cyclic and
1529          * hard to check for underflow.
1530          */
1531         errno = 0;
1532         result = acos(arg1);
1533         if (errno != 0)
1534                 ereport(ERROR,
1535                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1536                                  errmsg("input is out of range")));
1537
1538         CHECKFLOATVAL(result, isinf(arg1), true);
1539         PG_RETURN_FLOAT8(result);
1540 }
1541
1542
1543 /*
1544  *              dasin                   - returns the arcsin of arg1 (radians)
1545  */
1546 Datum
1547 dasin(PG_FUNCTION_ARGS)
1548 {
1549         float8          arg1 = PG_GETARG_FLOAT8(0);
1550         float8          result;
1551
1552         errno = 0;
1553         result = asin(arg1);
1554         if (errno != 0)
1555                 ereport(ERROR,
1556                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1557                                  errmsg("input is out of range")));
1558
1559         CHECKFLOATVAL(result, isinf(arg1), true);
1560         PG_RETURN_FLOAT8(result);
1561 }
1562
1563
1564 /*
1565  *              datan                   - returns the arctan of arg1 (radians)
1566  */
1567 Datum
1568 datan(PG_FUNCTION_ARGS)
1569 {
1570         float8          arg1 = PG_GETARG_FLOAT8(0);
1571         float8          result;
1572
1573         errno = 0;
1574         result = atan(arg1);
1575         if (errno != 0)
1576                 ereport(ERROR,
1577                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1578                                  errmsg("input is out of range")));
1579
1580         CHECKFLOATVAL(result, isinf(arg1), true);
1581         PG_RETURN_FLOAT8(result);
1582 }
1583
1584
1585 /*
1586  *              atan2                   - returns the arctan2 of arg1 (radians)
1587  */
1588 Datum
1589 datan2(PG_FUNCTION_ARGS)
1590 {
1591         float8          arg1 = PG_GETARG_FLOAT8(0);
1592         float8          arg2 = PG_GETARG_FLOAT8(1);
1593         float8          result;
1594
1595         errno = 0;
1596         result = atan2(arg1, arg2);
1597         if (errno != 0)
1598                 ereport(ERROR,
1599                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1600                                  errmsg("input is out of range")));
1601
1602         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
1603         PG_RETURN_FLOAT8(result);
1604 }
1605
1606
1607 /*
1608  *              dcos                    - returns the cosine of arg1 (radians)
1609  */
1610 Datum
1611 dcos(PG_FUNCTION_ARGS)
1612 {
1613         float8          arg1 = PG_GETARG_FLOAT8(0);
1614         float8          result;
1615
1616         errno = 0;
1617         result = cos(arg1);
1618         if (errno != 0)
1619                 ereport(ERROR,
1620                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1621                                  errmsg("input is out of range")));
1622
1623         CHECKFLOATVAL(result, isinf(arg1), true);
1624         PG_RETURN_FLOAT8(result);
1625 }
1626
1627
1628 /*
1629  *              dcot                    - returns the cotangent of arg1 (radians)
1630  */
1631 Datum
1632 dcot(PG_FUNCTION_ARGS)
1633 {
1634         float8          arg1 = PG_GETARG_FLOAT8(0);
1635         float8          result;
1636
1637         errno = 0;
1638         result = tan(arg1);
1639         if (errno != 0)
1640                 ereport(ERROR,
1641                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1642                                  errmsg("input is out of range")));
1643
1644         result = 1.0 / result;
1645         CHECKFLOATVAL(result, true /* cotan(pi/2) == inf */ , true);
1646         PG_RETURN_FLOAT8(result);
1647 }
1648
1649
1650 /*
1651  *              dsin                    - returns the sine of arg1 (radians)
1652  */
1653 Datum
1654 dsin(PG_FUNCTION_ARGS)
1655 {
1656         float8          arg1 = PG_GETARG_FLOAT8(0);
1657         float8          result;
1658
1659         errno = 0;
1660         result = sin(arg1);
1661         if (errno != 0)
1662                 ereport(ERROR,
1663                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1664                                  errmsg("input is out of range")));
1665
1666         CHECKFLOATVAL(result, isinf(arg1), true);
1667         PG_RETURN_FLOAT8(result);
1668 }
1669
1670
1671 /*
1672  *              dtan                    - returns the tangent of arg1 (radians)
1673  */
1674 Datum
1675 dtan(PG_FUNCTION_ARGS)
1676 {
1677         float8          arg1 = PG_GETARG_FLOAT8(0);
1678         float8          result;
1679
1680         errno = 0;
1681         result = tan(arg1);
1682         if (errno != 0)
1683                 ereport(ERROR,
1684                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1685                                  errmsg("input is out of range")));
1686
1687         CHECKFLOATVAL(result, true /* tan(pi/2) == Inf */ , true);
1688         PG_RETURN_FLOAT8(result);
1689 }
1690
1691
1692 /*
1693  *              degrees         - returns degrees converted from radians
1694  */
1695 Datum
1696 degrees(PG_FUNCTION_ARGS)
1697 {
1698         float8          arg1 = PG_GETARG_FLOAT8(0);
1699         float8          result;
1700
1701         result = arg1 * (180.0 / M_PI);
1702
1703         CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
1704         PG_RETURN_FLOAT8(result);
1705 }
1706
1707
1708 /*
1709  *              dpi                             - returns the constant PI
1710  */
1711 Datum
1712 dpi(PG_FUNCTION_ARGS)
1713 {
1714         PG_RETURN_FLOAT8(M_PI);
1715 }
1716
1717
1718 /*
1719  *              radians         - returns radians converted from degrees
1720  */
1721 Datum
1722 radians(PG_FUNCTION_ARGS)
1723 {
1724         float8          arg1 = PG_GETARG_FLOAT8(0);
1725         float8          result;
1726
1727         result = arg1 * (M_PI / 180.0);
1728
1729         CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
1730         PG_RETURN_FLOAT8(result);
1731 }
1732
1733
1734 /*
1735  *              drandom         - returns a random number
1736  */
1737 Datum
1738 drandom(PG_FUNCTION_ARGS)
1739 {
1740         float8          result;
1741
1742         /* result [0.0 - 1.0) */
1743         result = (double) random() / ((double) MAX_RANDOM_VALUE + 1);
1744
1745         PG_RETURN_FLOAT8(result);
1746 }
1747
1748
1749 /*
1750  *              setseed         - set seed for the random number generator
1751  */
1752 Datum
1753 setseed(PG_FUNCTION_ARGS)
1754 {
1755         float8          seed = PG_GETARG_FLOAT8(0);
1756         int                     iseed;
1757
1758         if (seed < -1 || seed > 1)
1759                 elog(ERROR, "setseed parameter %f out of range [-1,1]", seed);
1760
1761         iseed = (int) (seed * MAX_RANDOM_VALUE);
1762         srandom((unsigned int) iseed);
1763
1764         PG_RETURN_VOID();
1765 }
1766
1767
1768
1769 /*
1770  *              =========================
1771  *              FLOAT AGGREGATE OPERATORS
1772  *              =========================
1773  *
1774  *              float8_accum            - accumulate for AVG(), variance aggregates, etc.
1775  *              float4_accum            - same, but input data is float4
1776  *              float8_avg                      - produce final result for float AVG()
1777  *              float8_var_samp         - produce final result for float VAR_SAMP()
1778  *              float8_var_pop          - produce final result for float VAR_POP()
1779  *              float8_stddev_samp      - produce final result for float STDDEV_SAMP()
1780  *              float8_stddev_pop       - produce final result for float STDDEV_POP()
1781  *
1782  * The transition datatype for all these aggregates is a 3-element array
1783  * of float8, holding the values N, sum(X), sum(X*X) in that order.
1784  *
1785  * Note that we represent N as a float to avoid having to build a special
1786  * datatype.  Given a reasonable floating-point implementation, there should
1787  * be no accuracy loss unless N exceeds 2 ^ 52 or so (by which time the
1788  * user will have doubtless lost interest anyway...)
1789  */
1790
1791 static float8 *
1792 check_float8_array(ArrayType *transarray, const char *caller, int n)
1793 {
1794         /*
1795          * We expect the input to be an N-element float array; verify that. We
1796          * don't need to use deconstruct_array() since the array data is just
1797          * going to look like a C array of N float8 values.
1798          */
1799         if (ARR_NDIM(transarray) != 1 ||
1800                 ARR_DIMS(transarray)[0] != n ||
1801                 ARR_HASNULL(transarray) ||
1802                 ARR_ELEMTYPE(transarray) != FLOAT8OID)
1803                 elog(ERROR, "%s: expected %d-element float8 array", caller, n);
1804         return (float8 *) ARR_DATA_PTR(transarray);
1805 }
1806
1807 Datum
1808 float8_accum(PG_FUNCTION_ARGS)
1809 {
1810         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
1811         float8          newval = PG_GETARG_FLOAT8(1);
1812         float8     *transvalues;
1813         float8          N,
1814                                 sumX,
1815                                 sumX2;
1816
1817         transvalues = check_float8_array(transarray, "float8_accum", 3);
1818         N = transvalues[0];
1819         sumX = transvalues[1];
1820         sumX2 = transvalues[2];
1821
1822         N += 1.0;
1823         sumX += newval;
1824         CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newval), true);
1825         sumX2 += newval * newval;
1826         CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newval), true);
1827
1828         /*
1829          * If we're invoked as an aggregate, we can cheat and modify our first
1830          * parameter in-place to reduce palloc overhead. Otherwise we construct a
1831          * new array with the updated transition data and return it.
1832          */
1833         if (AggCheckCallContext(fcinfo, NULL))
1834         {
1835                 transvalues[0] = N;
1836                 transvalues[1] = sumX;
1837                 transvalues[2] = sumX2;
1838
1839                 PG_RETURN_ARRAYTYPE_P(transarray);
1840         }
1841         else
1842         {
1843                 Datum           transdatums[3];
1844                 ArrayType  *result;
1845
1846                 transdatums[0] = Float8GetDatumFast(N);
1847                 transdatums[1] = Float8GetDatumFast(sumX);
1848                 transdatums[2] = Float8GetDatumFast(sumX2);
1849
1850                 result = construct_array(transdatums, 3,
1851                                                                  FLOAT8OID,
1852                                                                  sizeof(float8), FLOAT8PASSBYVAL, 'd');
1853
1854                 PG_RETURN_ARRAYTYPE_P(result);
1855         }
1856 }
1857
1858 Datum
1859 float4_accum(PG_FUNCTION_ARGS)
1860 {
1861         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
1862
1863         /* do computations as float8 */
1864         float8          newval = PG_GETARG_FLOAT4(1);
1865         float8     *transvalues;
1866         float8          N,
1867                                 sumX,
1868                                 sumX2;
1869
1870         transvalues = check_float8_array(transarray, "float4_accum", 3);
1871         N = transvalues[0];
1872         sumX = transvalues[1];
1873         sumX2 = transvalues[2];
1874
1875         N += 1.0;
1876         sumX += newval;
1877         CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newval), true);
1878         sumX2 += newval * newval;
1879         CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newval), true);
1880
1881         /*
1882          * If we're invoked as an aggregate, we can cheat and modify our first
1883          * parameter in-place to reduce palloc overhead. Otherwise we construct a
1884          * new array with the updated transition data and return it.
1885          */
1886         if (AggCheckCallContext(fcinfo, NULL))
1887         {
1888                 transvalues[0] = N;
1889                 transvalues[1] = sumX;
1890                 transvalues[2] = sumX2;
1891
1892                 PG_RETURN_ARRAYTYPE_P(transarray);
1893         }
1894         else
1895         {
1896                 Datum           transdatums[3];
1897                 ArrayType  *result;
1898
1899                 transdatums[0] = Float8GetDatumFast(N);
1900                 transdatums[1] = Float8GetDatumFast(sumX);
1901                 transdatums[2] = Float8GetDatumFast(sumX2);
1902
1903                 result = construct_array(transdatums, 3,
1904                                                                  FLOAT8OID,
1905                                                                  sizeof(float8), FLOAT8PASSBYVAL, 'd');
1906
1907                 PG_RETURN_ARRAYTYPE_P(result);
1908         }
1909 }
1910
1911 Datum
1912 float8_avg(PG_FUNCTION_ARGS)
1913 {
1914         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
1915         float8     *transvalues;
1916         float8          N,
1917                                 sumX;
1918
1919         transvalues = check_float8_array(transarray, "float8_avg", 3);
1920         N = transvalues[0];
1921         sumX = transvalues[1];
1922         /* ignore sumX2 */
1923
1924         /* SQL defines AVG of no values to be NULL */
1925         if (N == 0.0)
1926                 PG_RETURN_NULL();
1927
1928         PG_RETURN_FLOAT8(sumX / N);
1929 }
1930
1931 Datum
1932 float8_var_pop(PG_FUNCTION_ARGS)
1933 {
1934         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
1935         float8     *transvalues;
1936         float8          N,
1937                                 sumX,
1938                                 sumX2,
1939                                 numerator;
1940
1941         transvalues = check_float8_array(transarray, "float8_var_pop", 3);
1942         N = transvalues[0];
1943         sumX = transvalues[1];
1944         sumX2 = transvalues[2];
1945
1946         /* Population variance is undefined when N is 0, so return NULL */
1947         if (N == 0.0)
1948                 PG_RETURN_NULL();
1949
1950         numerator = N * sumX2 - sumX * sumX;
1951         CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
1952
1953         /* Watch out for roundoff error producing a negative numerator */
1954         if (numerator <= 0.0)
1955                 PG_RETURN_FLOAT8(0.0);
1956
1957         PG_RETURN_FLOAT8(numerator / (N * N));
1958 }
1959
1960 Datum
1961 float8_var_samp(PG_FUNCTION_ARGS)
1962 {
1963         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
1964         float8     *transvalues;
1965         float8          N,
1966                                 sumX,
1967                                 sumX2,
1968                                 numerator;
1969
1970         transvalues = check_float8_array(transarray, "float8_var_samp", 3);
1971         N = transvalues[0];
1972         sumX = transvalues[1];
1973         sumX2 = transvalues[2];
1974
1975         /* Sample variance is undefined when N is 0 or 1, so return NULL */
1976         if (N <= 1.0)
1977                 PG_RETURN_NULL();
1978
1979         numerator = N * sumX2 - sumX * sumX;
1980         CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
1981
1982         /* Watch out for roundoff error producing a negative numerator */
1983         if (numerator <= 0.0)
1984                 PG_RETURN_FLOAT8(0.0);
1985
1986         PG_RETURN_FLOAT8(numerator / (N * (N - 1.0)));
1987 }
1988
1989 Datum
1990 float8_stddev_pop(PG_FUNCTION_ARGS)
1991 {
1992         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
1993         float8     *transvalues;
1994         float8          N,
1995                                 sumX,
1996                                 sumX2,
1997                                 numerator;
1998
1999         transvalues = check_float8_array(transarray, "float8_stddev_pop", 3);
2000         N = transvalues[0];
2001         sumX = transvalues[1];
2002         sumX2 = transvalues[2];
2003
2004         /* Population stddev is undefined when N is 0, so return NULL */
2005         if (N == 0.0)
2006                 PG_RETURN_NULL();
2007
2008         numerator = N * sumX2 - sumX * sumX;
2009         CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
2010
2011         /* Watch out for roundoff error producing a negative numerator */
2012         if (numerator <= 0.0)
2013                 PG_RETURN_FLOAT8(0.0);
2014
2015         PG_RETURN_FLOAT8(sqrt(numerator / (N * N)));
2016 }
2017
2018 Datum
2019 float8_stddev_samp(PG_FUNCTION_ARGS)
2020 {
2021         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2022         float8     *transvalues;
2023         float8          N,
2024                                 sumX,
2025                                 sumX2,
2026                                 numerator;
2027
2028         transvalues = check_float8_array(transarray, "float8_stddev_samp", 3);
2029         N = transvalues[0];
2030         sumX = transvalues[1];
2031         sumX2 = transvalues[2];
2032
2033         /* Sample stddev is undefined when N is 0 or 1, so return NULL */
2034         if (N <= 1.0)
2035                 PG_RETURN_NULL();
2036
2037         numerator = N * sumX2 - sumX * sumX;
2038         CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
2039
2040         /* Watch out for roundoff error producing a negative numerator */
2041         if (numerator <= 0.0)
2042                 PG_RETURN_FLOAT8(0.0);
2043
2044         PG_RETURN_FLOAT8(sqrt(numerator / (N * (N - 1.0))));
2045 }
2046
2047 /*
2048  *              =========================
2049  *              SQL2003 BINARY AGGREGATES
2050  *              =========================
2051  *
2052  * The transition datatype for all these aggregates is a 6-element array of
2053  * float8, holding the values N, sum(X), sum(X*X), sum(Y), sum(Y*Y), sum(X*Y)
2054  * in that order.  Note that Y is the first argument to the aggregates!
2055  *
2056  * It might seem attractive to optimize this by having multiple accumulator
2057  * functions that only calculate the sums actually needed.      But on most
2058  * modern machines, a couple of extra floating-point multiplies will be
2059  * insignificant compared to the other per-tuple overhead, so I've chosen
2060  * to minimize code space instead.
2061  */
2062
2063 Datum
2064 float8_regr_accum(PG_FUNCTION_ARGS)
2065 {
2066         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2067         float8          newvalY = PG_GETARG_FLOAT8(1);
2068         float8          newvalX = PG_GETARG_FLOAT8(2);
2069         float8     *transvalues;
2070         float8          N,
2071                                 sumX,
2072                                 sumX2,
2073                                 sumY,
2074                                 sumY2,
2075                                 sumXY;
2076
2077         transvalues = check_float8_array(transarray, "float8_regr_accum", 6);
2078         N = transvalues[0];
2079         sumX = transvalues[1];
2080         sumX2 = transvalues[2];
2081         sumY = transvalues[3];
2082         sumY2 = transvalues[4];
2083         sumXY = transvalues[5];
2084
2085         N += 1.0;
2086         sumX += newvalX;
2087         CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newvalX), true);
2088         sumX2 += newvalX * newvalX;
2089         CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newvalX), true);
2090         sumY += newvalY;
2091         CHECKFLOATVAL(sumY, isinf(transvalues[3]) || isinf(newvalY), true);
2092         sumY2 += newvalY * newvalY;
2093         CHECKFLOATVAL(sumY2, isinf(transvalues[4]) || isinf(newvalY), true);
2094         sumXY += newvalX * newvalY;
2095         CHECKFLOATVAL(sumXY, isinf(transvalues[5]) || isinf(newvalX) ||
2096                                   isinf(newvalY), true);
2097
2098         /*
2099          * If we're invoked as an aggregate, we can cheat and modify our first
2100          * parameter in-place to reduce palloc overhead. Otherwise we construct a
2101          * new array with the updated transition data and return it.
2102          */
2103         if (AggCheckCallContext(fcinfo, NULL))
2104         {
2105                 transvalues[0] = N;
2106                 transvalues[1] = sumX;
2107                 transvalues[2] = sumX2;
2108                 transvalues[3] = sumY;
2109                 transvalues[4] = sumY2;
2110                 transvalues[5] = sumXY;
2111
2112                 PG_RETURN_ARRAYTYPE_P(transarray);
2113         }
2114         else
2115         {
2116                 Datum           transdatums[6];
2117                 ArrayType  *result;
2118
2119                 transdatums[0] = Float8GetDatumFast(N);
2120                 transdatums[1] = Float8GetDatumFast(sumX);
2121                 transdatums[2] = Float8GetDatumFast(sumX2);
2122                 transdatums[3] = Float8GetDatumFast(sumY);
2123                 transdatums[4] = Float8GetDatumFast(sumY2);
2124                 transdatums[5] = Float8GetDatumFast(sumXY);
2125
2126                 result = construct_array(transdatums, 6,
2127                                                                  FLOAT8OID,
2128                                                                  sizeof(float8), FLOAT8PASSBYVAL, 'd');
2129
2130                 PG_RETURN_ARRAYTYPE_P(result);
2131         }
2132 }
2133
2134 Datum
2135 float8_regr_sxx(PG_FUNCTION_ARGS)
2136 {
2137         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2138         float8     *transvalues;
2139         float8          N,
2140                                 sumX,
2141                                 sumX2,
2142                                 numerator;
2143
2144         transvalues = check_float8_array(transarray, "float8_regr_sxx", 6);
2145         N = transvalues[0];
2146         sumX = transvalues[1];
2147         sumX2 = transvalues[2];
2148
2149         /* if N is 0 we should return NULL */
2150         if (N < 1.0)
2151                 PG_RETURN_NULL();
2152
2153         numerator = N * sumX2 - sumX * sumX;
2154         CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
2155
2156         /* Watch out for roundoff error producing a negative numerator */
2157         if (numerator <= 0.0)
2158                 PG_RETURN_FLOAT8(0.0);
2159
2160         PG_RETURN_FLOAT8(numerator / N);
2161 }
2162
2163 Datum
2164 float8_regr_syy(PG_FUNCTION_ARGS)
2165 {
2166         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2167         float8     *transvalues;
2168         float8          N,
2169                                 sumY,
2170                                 sumY2,
2171                                 numerator;
2172
2173         transvalues = check_float8_array(transarray, "float8_regr_syy", 6);
2174         N = transvalues[0];
2175         sumY = transvalues[3];
2176         sumY2 = transvalues[4];
2177
2178         /* if N is 0 we should return NULL */
2179         if (N < 1.0)
2180                 PG_RETURN_NULL();
2181
2182         numerator = N * sumY2 - sumY * sumY;
2183         CHECKFLOATVAL(numerator, isinf(sumY2) || isinf(sumY), true);
2184
2185         /* Watch out for roundoff error producing a negative numerator */
2186         if (numerator <= 0.0)
2187                 PG_RETURN_FLOAT8(0.0);
2188
2189         PG_RETURN_FLOAT8(numerator / N);
2190 }
2191
2192 Datum
2193 float8_regr_sxy(PG_FUNCTION_ARGS)
2194 {
2195         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2196         float8     *transvalues;
2197         float8          N,
2198                                 sumX,
2199                                 sumY,
2200                                 sumXY,
2201                                 numerator;
2202
2203         transvalues = check_float8_array(transarray, "float8_regr_sxy", 6);
2204         N = transvalues[0];
2205         sumX = transvalues[1];
2206         sumY = transvalues[3];
2207         sumXY = transvalues[5];
2208
2209         /* if N is 0 we should return NULL */
2210         if (N < 1.0)
2211                 PG_RETURN_NULL();
2212
2213         numerator = N * sumXY - sumX * sumY;
2214         CHECKFLOATVAL(numerator, isinf(sumXY) || isinf(sumX) ||
2215                                   isinf(sumY), true);
2216
2217         /* A negative result is valid here */
2218
2219         PG_RETURN_FLOAT8(numerator / N);
2220 }
2221
2222 Datum
2223 float8_regr_avgx(PG_FUNCTION_ARGS)
2224 {
2225         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2226         float8     *transvalues;
2227         float8          N,
2228                                 sumX;
2229
2230         transvalues = check_float8_array(transarray, "float8_regr_avgx", 6);
2231         N = transvalues[0];
2232         sumX = transvalues[1];
2233
2234         /* if N is 0 we should return NULL */
2235         if (N < 1.0)
2236                 PG_RETURN_NULL();
2237
2238         PG_RETURN_FLOAT8(sumX / N);
2239 }
2240
2241 Datum
2242 float8_regr_avgy(PG_FUNCTION_ARGS)
2243 {
2244         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2245         float8     *transvalues;
2246         float8          N,
2247                                 sumY;
2248
2249         transvalues = check_float8_array(transarray, "float8_regr_avgy", 6);
2250         N = transvalues[0];
2251         sumY = transvalues[3];
2252
2253         /* if N is 0 we should return NULL */
2254         if (N < 1.0)
2255                 PG_RETURN_NULL();
2256
2257         PG_RETURN_FLOAT8(sumY / N);
2258 }
2259
2260 Datum
2261 float8_covar_pop(PG_FUNCTION_ARGS)
2262 {
2263         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2264         float8     *transvalues;
2265         float8          N,
2266                                 sumX,
2267                                 sumY,
2268                                 sumXY,
2269                                 numerator;
2270
2271         transvalues = check_float8_array(transarray, "float8_covar_pop", 6);
2272         N = transvalues[0];
2273         sumX = transvalues[1];
2274         sumY = transvalues[3];
2275         sumXY = transvalues[5];
2276
2277         /* if N is 0 we should return NULL */
2278         if (N < 1.0)
2279                 PG_RETURN_NULL();
2280
2281         numerator = N * sumXY - sumX * sumY;
2282         CHECKFLOATVAL(numerator, isinf(sumXY) || isinf(sumX) ||
2283                                   isinf(sumY), true);
2284
2285         PG_RETURN_FLOAT8(numerator / (N * N));
2286 }
2287
2288 Datum
2289 float8_covar_samp(PG_FUNCTION_ARGS)
2290 {
2291         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2292         float8     *transvalues;
2293         float8          N,
2294                                 sumX,
2295                                 sumY,
2296                                 sumXY,
2297                                 numerator;
2298
2299         transvalues = check_float8_array(transarray, "float8_covar_samp", 6);
2300         N = transvalues[0];
2301         sumX = transvalues[1];
2302         sumY = transvalues[3];
2303         sumXY = transvalues[5];
2304
2305         /* if N is <= 1 we should return NULL */
2306         if (N < 2.0)
2307                 PG_RETURN_NULL();
2308
2309         numerator = N * sumXY - sumX * sumY;
2310         CHECKFLOATVAL(numerator, isinf(sumXY) || isinf(sumX) ||
2311                                   isinf(sumY), true);
2312
2313         PG_RETURN_FLOAT8(numerator / (N * (N - 1.0)));
2314 }
2315
2316 Datum
2317 float8_corr(PG_FUNCTION_ARGS)
2318 {
2319         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2320         float8     *transvalues;
2321         float8          N,
2322                                 sumX,
2323                                 sumX2,
2324                                 sumY,
2325                                 sumY2,
2326                                 sumXY,
2327                                 numeratorX,
2328                                 numeratorY,
2329                                 numeratorXY;
2330
2331         transvalues = check_float8_array(transarray, "float8_corr", 6);
2332         N = transvalues[0];
2333         sumX = transvalues[1];
2334         sumX2 = transvalues[2];
2335         sumY = transvalues[3];
2336         sumY2 = transvalues[4];
2337         sumXY = transvalues[5];
2338
2339         /* if N is 0 we should return NULL */
2340         if (N < 1.0)
2341                 PG_RETURN_NULL();
2342
2343         numeratorX = N * sumX2 - sumX * sumX;
2344         CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true);
2345         numeratorY = N * sumY2 - sumY * sumY;
2346         CHECKFLOATVAL(numeratorY, isinf(sumY2) || isinf(sumY), true);
2347         numeratorXY = N * sumXY - sumX * sumY;
2348         CHECKFLOATVAL(numeratorXY, isinf(sumXY) || isinf(sumX) ||
2349                                   isinf(sumY), true);
2350         if (numeratorX <= 0 || numeratorY <= 0)
2351                 PG_RETURN_NULL();
2352
2353         PG_RETURN_FLOAT8(numeratorXY / sqrt(numeratorX * numeratorY));
2354 }
2355
2356 Datum
2357 float8_regr_r2(PG_FUNCTION_ARGS)
2358 {
2359         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2360         float8     *transvalues;
2361         float8          N,
2362                                 sumX,
2363                                 sumX2,
2364                                 sumY,
2365                                 sumY2,
2366                                 sumXY,
2367                                 numeratorX,
2368                                 numeratorY,
2369                                 numeratorXY;
2370
2371         transvalues = check_float8_array(transarray, "float8_regr_r2", 6);
2372         N = transvalues[0];
2373         sumX = transvalues[1];
2374         sumX2 = transvalues[2];
2375         sumY = transvalues[3];
2376         sumY2 = transvalues[4];
2377         sumXY = transvalues[5];
2378
2379         /* if N is 0 we should return NULL */
2380         if (N < 1.0)
2381                 PG_RETURN_NULL();
2382
2383         numeratorX = N * sumX2 - sumX * sumX;
2384         CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true);
2385         numeratorY = N * sumY2 - sumY * sumY;
2386         CHECKFLOATVAL(numeratorY, isinf(sumY2) || isinf(sumY), true);
2387         numeratorXY = N * sumXY - sumX * sumY;
2388         CHECKFLOATVAL(numeratorXY, isinf(sumXY) || isinf(sumX) ||
2389                                   isinf(sumY), true);
2390         if (numeratorX <= 0)
2391                 PG_RETURN_NULL();
2392         /* per spec, horizontal line produces 1.0 */
2393         if (numeratorY <= 0)
2394                 PG_RETURN_FLOAT8(1.0);
2395
2396         PG_RETURN_FLOAT8((numeratorXY * numeratorXY) /
2397                                          (numeratorX * numeratorY));
2398 }
2399
2400 Datum
2401 float8_regr_slope(PG_FUNCTION_ARGS)
2402 {
2403         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2404         float8     *transvalues;
2405         float8          N,
2406                                 sumX,
2407                                 sumX2,
2408                                 sumY,
2409                                 sumXY,
2410                                 numeratorX,
2411                                 numeratorXY;
2412
2413         transvalues = check_float8_array(transarray, "float8_regr_slope", 6);
2414         N = transvalues[0];
2415         sumX = transvalues[1];
2416         sumX2 = transvalues[2];
2417         sumY = transvalues[3];
2418         sumXY = transvalues[5];
2419
2420         /* if N is 0 we should return NULL */
2421         if (N < 1.0)
2422                 PG_RETURN_NULL();
2423
2424         numeratorX = N * sumX2 - sumX * sumX;
2425         CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true);
2426         numeratorXY = N * sumXY - sumX * sumY;
2427         CHECKFLOATVAL(numeratorXY, isinf(sumXY) || isinf(sumX) ||
2428                                   isinf(sumY), true);
2429         if (numeratorX <= 0)
2430                 PG_RETURN_NULL();
2431
2432         PG_RETURN_FLOAT8(numeratorXY / numeratorX);
2433 }
2434
2435 Datum
2436 float8_regr_intercept(PG_FUNCTION_ARGS)
2437 {
2438         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2439         float8     *transvalues;
2440         float8          N,
2441                                 sumX,
2442                                 sumX2,
2443                                 sumY,
2444                                 sumXY,
2445                                 numeratorX,
2446                                 numeratorXXY;
2447
2448         transvalues = check_float8_array(transarray, "float8_regr_intercept", 6);
2449         N = transvalues[0];
2450         sumX = transvalues[1];
2451         sumX2 = transvalues[2];
2452         sumY = transvalues[3];
2453         sumXY = transvalues[5];
2454
2455         /* if N is 0 we should return NULL */
2456         if (N < 1.0)
2457                 PG_RETURN_NULL();
2458
2459         numeratorX = N * sumX2 - sumX * sumX;
2460         CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true);
2461         numeratorXXY = sumY * sumX2 - sumX * sumXY;
2462         CHECKFLOATVAL(numeratorXXY, isinf(sumY) || isinf(sumX2) ||
2463                                   isinf(sumX) || isinf(sumXY), true);
2464         if (numeratorX <= 0)
2465                 PG_RETURN_NULL();
2466
2467         PG_RETURN_FLOAT8(numeratorXXY / numeratorX);
2468 }
2469
2470
2471 /*
2472  *              ====================================
2473  *              MIXED-PRECISION ARITHMETIC OPERATORS
2474  *              ====================================
2475  */
2476
2477 /*
2478  *              float48pl               - returns arg1 + arg2
2479  *              float48mi               - returns arg1 - arg2
2480  *              float48mul              - returns arg1 * arg2
2481  *              float48div              - returns arg1 / arg2
2482  */
2483 Datum
2484 float48pl(PG_FUNCTION_ARGS)
2485 {
2486         float4          arg1 = PG_GETARG_FLOAT4(0);
2487         float8          arg2 = PG_GETARG_FLOAT8(1);
2488         float8          result;
2489
2490         result = arg1 + arg2;
2491         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
2492         PG_RETURN_FLOAT8(result);
2493 }
2494
2495 Datum
2496 float48mi(PG_FUNCTION_ARGS)
2497 {
2498         float4          arg1 = PG_GETARG_FLOAT4(0);
2499         float8          arg2 = PG_GETARG_FLOAT8(1);
2500         float8          result;
2501
2502         result = arg1 - arg2;
2503         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
2504         PG_RETURN_FLOAT8(result);
2505 }
2506
2507 Datum
2508 float48mul(PG_FUNCTION_ARGS)
2509 {
2510         float4          arg1 = PG_GETARG_FLOAT4(0);
2511         float8          arg2 = PG_GETARG_FLOAT8(1);
2512         float8          result;
2513
2514         result = arg1 * arg2;
2515         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
2516                                   arg1 == 0 || arg2 == 0);
2517         PG_RETURN_FLOAT8(result);
2518 }
2519
2520 Datum
2521 float48div(PG_FUNCTION_ARGS)
2522 {
2523         float4          arg1 = PG_GETARG_FLOAT4(0);
2524         float8          arg2 = PG_GETARG_FLOAT8(1);
2525         float8          result;
2526
2527         if (arg2 == 0.0)
2528                 ereport(ERROR,
2529                                 (errcode(ERRCODE_DIVISION_BY_ZERO),
2530                                  errmsg("division by zero")));
2531
2532         result = arg1 / arg2;
2533         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
2534         PG_RETURN_FLOAT8(result);
2535 }
2536
2537 /*
2538  *              float84pl               - returns arg1 + arg2
2539  *              float84mi               - returns arg1 - arg2
2540  *              float84mul              - returns arg1 * arg2
2541  *              float84div              - returns arg1 / arg2
2542  */
2543 Datum
2544 float84pl(PG_FUNCTION_ARGS)
2545 {
2546         float8          arg1 = PG_GETARG_FLOAT8(0);
2547         float4          arg2 = PG_GETARG_FLOAT4(1);
2548         float8          result;
2549
2550         result = arg1 + arg2;
2551
2552         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
2553         PG_RETURN_FLOAT8(result);
2554 }
2555
2556 Datum
2557 float84mi(PG_FUNCTION_ARGS)
2558 {
2559         float8          arg1 = PG_GETARG_FLOAT8(0);
2560         float4          arg2 = PG_GETARG_FLOAT4(1);
2561         float8          result;
2562
2563         result = arg1 - arg2;
2564
2565         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
2566         PG_RETURN_FLOAT8(result);
2567 }
2568
2569 Datum
2570 float84mul(PG_FUNCTION_ARGS)
2571 {
2572         float8          arg1 = PG_GETARG_FLOAT8(0);
2573         float4          arg2 = PG_GETARG_FLOAT4(1);
2574         float8          result;
2575
2576         result = arg1 * arg2;
2577
2578         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
2579                                   arg1 == 0 || arg2 == 0);
2580         PG_RETURN_FLOAT8(result);
2581 }
2582
2583 Datum
2584 float84div(PG_FUNCTION_ARGS)
2585 {
2586         float8          arg1 = PG_GETARG_FLOAT8(0);
2587         float4          arg2 = PG_GETARG_FLOAT4(1);
2588         float8          result;
2589
2590         if (arg2 == 0.0)
2591                 ereport(ERROR,
2592                                 (errcode(ERRCODE_DIVISION_BY_ZERO),
2593                                  errmsg("division by zero")));
2594
2595         result = arg1 / arg2;
2596
2597         CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
2598         PG_RETURN_FLOAT8(result);
2599 }
2600
2601 /*
2602  *              ====================
2603  *              COMPARISON OPERATORS
2604  *              ====================
2605  */
2606
2607 /*
2608  *              float48{eq,ne,lt,le,gt,ge}              - float4/float8 comparison operations
2609  */
2610 Datum
2611 float48eq(PG_FUNCTION_ARGS)
2612 {
2613         float4          arg1 = PG_GETARG_FLOAT4(0);
2614         float8          arg2 = PG_GETARG_FLOAT8(1);
2615
2616         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
2617 }
2618
2619 Datum
2620 float48ne(PG_FUNCTION_ARGS)
2621 {
2622         float4          arg1 = PG_GETARG_FLOAT4(0);
2623         float8          arg2 = PG_GETARG_FLOAT8(1);
2624
2625         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
2626 }
2627
2628 Datum
2629 float48lt(PG_FUNCTION_ARGS)
2630 {
2631         float4          arg1 = PG_GETARG_FLOAT4(0);
2632         float8          arg2 = PG_GETARG_FLOAT8(1);
2633
2634         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
2635 }
2636
2637 Datum
2638 float48le(PG_FUNCTION_ARGS)
2639 {
2640         float4          arg1 = PG_GETARG_FLOAT4(0);
2641         float8          arg2 = PG_GETARG_FLOAT8(1);
2642
2643         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
2644 }
2645
2646 Datum
2647 float48gt(PG_FUNCTION_ARGS)
2648 {
2649         float4          arg1 = PG_GETARG_FLOAT4(0);
2650         float8          arg2 = PG_GETARG_FLOAT8(1);
2651
2652         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
2653 }
2654
2655 Datum
2656 float48ge(PG_FUNCTION_ARGS)
2657 {
2658         float4          arg1 = PG_GETARG_FLOAT4(0);
2659         float8          arg2 = PG_GETARG_FLOAT8(1);
2660
2661         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
2662 }
2663
2664 /*
2665  *              float84{eq,ne,lt,le,gt,ge}              - float8/float4 comparison operations
2666  */
2667 Datum
2668 float84eq(PG_FUNCTION_ARGS)
2669 {
2670         float8          arg1 = PG_GETARG_FLOAT8(0);
2671         float4          arg2 = PG_GETARG_FLOAT4(1);
2672
2673         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
2674 }
2675
2676 Datum
2677 float84ne(PG_FUNCTION_ARGS)
2678 {
2679         float8          arg1 = PG_GETARG_FLOAT8(0);
2680         float4          arg2 = PG_GETARG_FLOAT4(1);
2681
2682         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
2683 }
2684
2685 Datum
2686 float84lt(PG_FUNCTION_ARGS)
2687 {
2688         float8          arg1 = PG_GETARG_FLOAT8(0);
2689         float4          arg2 = PG_GETARG_FLOAT4(1);
2690
2691         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
2692 }
2693
2694 Datum
2695 float84le(PG_FUNCTION_ARGS)
2696 {
2697         float8          arg1 = PG_GETARG_FLOAT8(0);
2698         float4          arg2 = PG_GETARG_FLOAT4(1);
2699
2700         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
2701 }
2702
2703 Datum
2704 float84gt(PG_FUNCTION_ARGS)
2705 {
2706         float8          arg1 = PG_GETARG_FLOAT8(0);
2707         float4          arg2 = PG_GETARG_FLOAT4(1);
2708
2709         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
2710 }
2711
2712 Datum
2713 float84ge(PG_FUNCTION_ARGS)
2714 {
2715         float8          arg1 = PG_GETARG_FLOAT8(0);
2716         float4          arg2 = PG_GETARG_FLOAT4(1);
2717
2718         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
2719 }
2720
2721 /*
2722  * Implements the float8 version of the width_bucket() function
2723  * defined by SQL2003. See also width_bucket_numeric().
2724  *
2725  * 'bound1' and 'bound2' are the lower and upper bounds of the
2726  * histogram's range, respectively. 'count' is the number of buckets
2727  * in the histogram. width_bucket() returns an integer indicating the
2728  * bucket number that 'operand' belongs to in an equiwidth histogram
2729  * with the specified characteristics. An operand smaller than the
2730  * lower bound is assigned to bucket 0. An operand greater than the
2731  * upper bound is assigned to an additional bucket (with number
2732  * count+1). We don't allow "NaN" for any of the float8 inputs, and we
2733  * don't allow either of the histogram bounds to be +/- infinity.
2734  */
2735 Datum
2736 width_bucket_float8(PG_FUNCTION_ARGS)
2737 {
2738         float8          operand = PG_GETARG_FLOAT8(0);
2739         float8          bound1 = PG_GETARG_FLOAT8(1);
2740         float8          bound2 = PG_GETARG_FLOAT8(2);
2741         int32           count = PG_GETARG_INT32(3);
2742         int32           result;
2743
2744         if (count <= 0.0)
2745                 ereport(ERROR,
2746                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
2747                                  errmsg("count must be greater than zero")));
2748
2749         if (isnan(operand) || isnan(bound1) || isnan(bound2))
2750                 ereport(ERROR,
2751                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
2752                          errmsg("operand, lower bound, and upper bound cannot be NaN")));
2753
2754         /* Note that we allow "operand" to be infinite */
2755         if (isinf(bound1) || isinf(bound2))
2756                 ereport(ERROR,
2757                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
2758                                  errmsg("lower and upper bounds must be finite")));
2759
2760         if (bound1 < bound2)
2761         {
2762                 if (operand < bound1)
2763                         result = 0;
2764                 else if (operand >= bound2)
2765                 {
2766                         result = count + 1;
2767                         /* check for overflow */
2768                         if (result < count)
2769                                 ereport(ERROR,
2770                                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
2771                                                  errmsg("integer out of range")));
2772                 }
2773                 else
2774                         result = ((float8) count * (operand - bound1) / (bound2 - bound1)) + 1;
2775         }
2776         else if (bound1 > bound2)
2777         {
2778                 if (operand > bound1)
2779                         result = 0;
2780                 else if (operand <= bound2)
2781                 {
2782                         result = count + 1;
2783                         /* check for overflow */
2784                         if (result < count)
2785                                 ereport(ERROR,
2786                                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
2787                                                  errmsg("integer out of range")));
2788                 }
2789                 else
2790                         result = ((float8) count * (bound1 - operand) / (bound1 - bound2)) + 1;
2791         }
2792         else
2793         {
2794                 ereport(ERROR,
2795                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
2796                                  errmsg("lower bound cannot equal upper bound")));
2797                 result = 0;                             /* keep the compiler quiet */
2798         }
2799
2800         PG_RETURN_INT32(result);
2801 }
2802
2803 /* ========== PRIVATE ROUTINES ========== */
2804
2805 #ifndef HAVE_CBRT
2806
2807 static double
2808 cbrt(double x)
2809 {
2810         int                     isneg = (x < 0.0);
2811         double          absx = fabs(x);
2812         double          tmpres = pow(absx, (double) 1.0 / (double) 3.0);
2813
2814         /*
2815          * The result is somewhat inaccurate --- not really pow()'s fault, as the
2816          * exponent it's handed contains roundoff error.  We can improve the
2817          * accuracy by doing one iteration of Newton's formula.  Beware of zero
2818          * input however.
2819          */
2820         if (tmpres > 0.0)
2821                 tmpres -= (tmpres - absx / (tmpres * tmpres)) / (double) 3.0;
2822
2823         return isneg ? -tmpres : tmpres;
2824 }
2825
2826 #endif   /* !HAVE_CBRT */