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