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