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