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-2002, PostgreSQL Global Development Group
   7  * Portions Copyright (c) 1994, Regents of the University of California
   8  *
   9  *
  10  * IDENTIFICATION
  11  *        $Header: /cvsroot/pgsql/src/backend/utils/adt/float.c,v 1.81 2002/09/04 20:31:27 momjian Exp $
  12  *
  13  *-------------------------------------------------------------------------
  14  */
  15 /*----------
  16  * OLD COMMENTS
  17  *              Basic float4 ops:
  18  *               float4in, float4out, float4abs, float4um, float4up
  19  *              Basic float8 ops:
  20  *               float8in, float8out, float8abs, float8um, float8up
  21  *              Arithmetic operators:
  22  *               float4pl, float4mi, float4mul, float4div
  23  *               float8pl, float8mi, float8mul, float8div
  24  *              Comparison operators:
  25  *               float4eq, float4ne, float4lt, float4le, float4gt, float4ge, float4cmp
  26  *               float8eq, float8ne, float8lt, float8le, float8gt, float8ge, float8cmp
  27  *              Conversion routines:
  28  *               ftod, dtof, i4tod, dtoi4, i2tod, dtoi2, itof, ftoi, i2tof, ftoi2
  29  *
  30  *              Random float8 ops:
  31  *               dround, dtrunc, dsqrt, dcbrt, dpow, dexp, dlog1
  32  *              Arithmetic operators:
  33  *               float48pl, float48mi, float48mul, float48div
  34  *               float84pl, float84mi, float84mul, float84div
  35  *              Comparison operators:
  36  *               float48eq, float48ne, float48lt, float48le, float48gt, float48ge
  37  *               float84eq, float84ne, float84lt, float84le, float84gt, float84ge
  38  *
  39  *              (You can do the arithmetic and comparison stuff using conversion
  40  *               routines, but then you pay the overhead of invoking a separate
  41  *               conversion function...)
  42  *
  43  * XXX GLUESOME STUFF. FIX IT! -AY '94
  44  *
  45  *              Added some additional conversion routines and cleaned up
  46  *               a bit of the existing code. Need to change the error checking
  47  *               for calls to pow(), exp() since on some machines (my Linux box
  48  *               included) these routines do not set errno. - tgl 97/05/10
  49  *----------
  50  */
  51 #include "postgres.h"
  52
  53 #include <ctype.h>
  54 #include <errno.h>
  55 #include <float.h>                              /* faked on sunos4 */
  56 #include <math.h>
  57
  58 #include <limits.h>
  59 /* for finite() on Solaris */
  60 #ifdef HAVE_IEEEFP_H
  61 #include <ieeefp.h>
  62 #endif
  63
  64 #include "catalog/pg_type.h"
  65 #include "fmgr.h"
  66 #include "utils/array.h"
  67 #include "utils/builtins.h"
  68
  69
  70 #if !(NeXT && NX_CURRENT_COMPILER_RELEASE > NX_COMPILER_RELEASE_3_2)
  71  /* NS3.3 has conflicting declarations of these in <math.h> */
  72
  73 #ifndef atof
  74 extern double atof(const char *p);
  75 #endif
  76
  77 #ifndef HAVE_CBRT
  78 #define cbrt my_cbrt
  79 static double cbrt(double x);
  80
  81 #else
  82 #if !defined(nextstep)
  83 extern double cbrt(double x);
  84 #endif
  85 #endif   /* HAVE_CBRT */
  86
  87 #ifndef HAVE_RINT
  88 #define rint my_rint
  89 static double rint(double x);
  90
  91 #else
  92 extern double rint(double x);
  93 #endif   /* HAVE_RINT */
  94 #endif   /* NeXT check */
  95
  96
  97 static void CheckFloat4Val(double val);
  98 static void CheckFloat8Val(double val);
  99
 100 #ifndef M_PI
 101 /* from my RH5.2 gcc math.h file - thomas 2000-04-03 */
 102 #define M_PI 3.14159265358979323846
 103 #endif
 104
 105 #ifndef NAN
 106 #define NAN             (0.0/0.0)
 107 #endif
 108
 109 #ifndef SHRT_MAX
 110 #define SHRT_MAX 32767
 111 #endif
 112 #ifndef SHRT_MIN
 113 #define SHRT_MIN (-32768)
 114 #endif
 115
 116 #define FORMAT                  'g'             /* use "g" output format as standard
 117                                                                  * format */
 118 /* not sure what the following should be, but better to make it over-sufficient */
 119 #define MAXFLOATWIDTH   64
 120 #define MAXDOUBLEWIDTH  128
 121
 122 /* ========== USER I/O ROUTINES ========== */
 123
 124
 125 #define FLOAT4_MAX               FLT_MAX
 126 #define FLOAT4_MIN               FLT_MIN
 127 #define FLOAT8_MAX               DBL_MAX
 128 #define FLOAT8_MIN               DBL_MIN
 129
 130
 131 /*
 132    check to see if a float4 val is outside of
 133    the FLOAT4_MIN, FLOAT4_MAX bounds.
 134
 135    raise an elog warning if it is
 136 */
 137 static void
 138 CheckFloat4Val(double val)
 139 {
 140         /*
 141          * defining unsafe floats's will make float4 and float8 ops faster at
 142          * the cost of safety, of course!
 143          */
 144 #ifdef UNSAFE_FLOATS
 145         return;
 146 #else
 147         if (fabs(val) > FLOAT4_MAX)
 148                 elog(ERROR, "Bad float4 input format -- overflow");
 149         if (val != 0.0 && fabs(val) < FLOAT4_MIN)
 150                 elog(ERROR, "Bad float4 input format -- underflow");
 151         return;
 152 #endif   /* UNSAFE_FLOATS */
 153 }
 154
 155 /*
 156    check to see if a float8 val is outside of
 157    the FLOAT8_MIN, FLOAT8_MAX bounds.
 158
 159    raise an elog warning if it is
 160 */
 161 static void
 162 CheckFloat8Val(double val)
 163 {
 164         /*
 165          * defining unsafe floats's will make float4 and float8 ops faster at
 166          * the cost of safety, of course!
 167          */
 168 #ifdef UNSAFE_FLOATS
 169         return;
 170 #else
 171         if (fabs(val) > FLOAT8_MAX)
 172                 elog(ERROR, "Bad float8 input format -- overflow");
 173         if (val != 0.0 && fabs(val) < FLOAT8_MIN)
 174                 elog(ERROR, "Bad float8 input format -- underflow");
 175         return;
 176 #endif   /* UNSAFE_FLOATS */
 177 }
 178
 179 /*
 180  *              float4in                - converts "num" to float
 181  *                                                restricted syntax:
 182  *                                                {<sp>} [+|-] {digit} [.{digit}] [<exp>]
 183  *                                                where <sp> is a space, digit is 0-9,
 184  *                                                <exp> is "e" or "E" followed by an integer.
 185  */
 186 Datum
 187 float4in(PG_FUNCTION_ARGS)
 188 {
 189         char       *num = PG_GETARG_CSTRING(0);
 190         double          val;
 191         char       *endptr;
 192
 193         errno = 0;
 194         val = strtod(num, &endptr);
 195         if (*endptr != '\0')
 196         {
 197                 /*
 198                  * XXX we should accept "Infinity" and "-Infinity" too, but what
 199                  * are the correct values to assign?  HUGE_VAL will provoke an
 200                  * error from CheckFloat4Val.
 201                  */
 202                 if (strcasecmp(num, "NaN") == 0)
 203                         val = NAN;
 204                 else
 205                         elog(ERROR, "Bad float4 input format '%s'", num);
 206         }
 207         else
 208         {
 209                 if (errno == ERANGE)
 210                         elog(ERROR, "Input '%s' is out of range for float4", num);
 211         }
 212
 213         /*
 214          * if we get here, we have a legal double, still need to check to see
 215          * if it's a legal float
 216          */
 217         CheckFloat4Val(val);
 218
 219         PG_RETURN_FLOAT4((float4) val);
 220 }
 221
 222 /*
 223  *              float4out               - converts a float4 number to a string
 224  *                                                using a standard output format
 225  */
 226 Datum
 227 float4out(PG_FUNCTION_ARGS)
 228 {
 229         float4          num = PG_GETARG_FLOAT4(0);
 230         char       *ascii = (char *) palloc(MAXFLOATWIDTH + 1);
 231         int                     infflag;
 232
 233         if (isnan(num))
 234                 PG_RETURN_CSTRING(strcpy(ascii, "NaN"));
 235         infflag = isinf(num);
 236         if (infflag > 0)
 237                 PG_RETURN_CSTRING(strcpy(ascii, "Infinity"));
 238         if (infflag < 0)
 239                 PG_RETURN_CSTRING(strcpy(ascii, "-Infinity"));
 240
 241         sprintf(ascii, "%.*g", FLT_DIG, num);
 242         PG_RETURN_CSTRING(ascii);
 243 }
 244
 245 /*
 246  *              float8in                - converts "num" to float8
 247  *                                                restricted syntax:
 248  *                                                {<sp>} [+|-] {digit} [.{digit}] [<exp>]
 249  *                                                where <sp> is a space, digit is 0-9,
 250  *                                                <exp> is "e" or "E" followed by an integer.
 251  */
 252 Datum
 253 float8in(PG_FUNCTION_ARGS)
 254 {
 255         char       *num = PG_GETARG_CSTRING(0);
 256         double          val;
 257         char       *endptr;
 258
 259         errno = 0;
 260         val = strtod(num, &endptr);
 261         if (*endptr != '\0')
 262         {
 263                 if (strcasecmp(num, "NaN") == 0)
 264                         val = NAN;
 265                 else if (strcasecmp(num, "Infinity") == 0)
 266                         val = HUGE_VAL;
 267                 else if (strcasecmp(num, "-Infinity") == 0)
 268                         val = -HUGE_VAL;
 269                 else
 270                         elog(ERROR, "Bad float8 input format '%s'", num);
 271         }
 272         else
 273         {
 274                 if (errno == ERANGE)
 275                         elog(ERROR, "Input '%s' is out of range for float8", num);
 276         }
 277
 278         CheckFloat8Val(val);
 279
 280         PG_RETURN_FLOAT8(val);
 281 }
 282
 283
 284 /*
 285  *              float8out               - converts float8 number to a string
 286  *                                                using a standard output format
 287  */
 288 Datum
 289 float8out(PG_FUNCTION_ARGS)
 290 {
 291         float8          num = PG_GETARG_FLOAT8(0);
 292         char       *ascii = (char *) palloc(MAXDOUBLEWIDTH + 1);
 293         int                     infflag;
 294
 295         if (isnan(num))
 296                 PG_RETURN_CSTRING(strcpy(ascii, "NaN"));
 297         infflag = isinf(num);
 298         if (infflag > 0)
 299                 PG_RETURN_CSTRING(strcpy(ascii, "Infinity"));
 300         if (infflag < 0)
 301                 PG_RETURN_CSTRING(strcpy(ascii, "-Infinity"));
 302
 303         sprintf(ascii, "%.*g", DBL_DIG, num);
 304         PG_RETURN_CSTRING(ascii);
 305 }
 306
 307 /* ========== PUBLIC ROUTINES ========== */
 308
 309
 310 /*
 311  *              ======================
 312  *              FLOAT4 BASE OPERATIONS
 313  *              ======================
 314  */
 315
 316 /*
 317  *              float4abs               - returns |arg1| (absolute value)
 318  */
 319 Datum
 320 float4abs(PG_FUNCTION_ARGS)
 321 {
 322         float4          arg1 = PG_GETARG_FLOAT4(0);
 323
 324         PG_RETURN_FLOAT4((float4) fabs(arg1));
 325 }
 326
 327 /*
 328  *              float4um                - returns -arg1 (unary minus)
 329  */
 330 Datum
 331 float4um(PG_FUNCTION_ARGS)
 332 {
 333         float4          arg1 = PG_GETARG_FLOAT4(0);
 334
 335         PG_RETURN_FLOAT4((float4) -arg1);
 336 }
 337
 338 Datum
 339 float4up(PG_FUNCTION_ARGS)
 340 {
 341         float4          arg = PG_GETARG_FLOAT4(0);
 342
 343         PG_RETURN_FLOAT4(arg);
 344 }
 345
 346 Datum
 347 float4larger(PG_FUNCTION_ARGS)
 348 {
 349         float4          arg1 = PG_GETARG_FLOAT4(0);
 350         float4          arg2 = PG_GETARG_FLOAT4(1);
 351         float4          result;
 352
 353         result = ((arg1 > arg2) ? arg1 : arg2);
 354         PG_RETURN_FLOAT4(result);
 355 }
 356
 357 Datum
 358 float4smaller(PG_FUNCTION_ARGS)
 359 {
 360         float4          arg1 = PG_GETARG_FLOAT4(0);
 361         float4          arg2 = PG_GETARG_FLOAT4(1);
 362         float4          result;
 363
 364         result = ((arg1 < arg2) ? arg1 : arg2);
 365         PG_RETURN_FLOAT4(result);
 366 }
 367
 368 /*
 369  *              ======================
 370  *              FLOAT8 BASE OPERATIONS
 371  *              ======================
 372  */
 373
 374 /*
 375  *              float8abs               - returns |arg1| (absolute value)
 376  */
 377 Datum
 378 float8abs(PG_FUNCTION_ARGS)
 379 {
 380         float8          arg1 = PG_GETARG_FLOAT8(0);
 381         float8          result;
 382
 383         result = fabs(arg1);
 384
 385         CheckFloat8Val(result);
 386         PG_RETURN_FLOAT8(result);
 387 }
 388
 389
 390 /*
 391  *              float8um                - returns -arg1 (unary minus)
 392  */
 393 Datum
 394 float8um(PG_FUNCTION_ARGS)
 395 {
 396         float8          arg1 = PG_GETARG_FLOAT8(0);
 397         float8          result;
 398
 399         result = ((arg1 != 0) ? -(arg1) : arg1);
 400
 401         CheckFloat8Val(result);
 402         PG_RETURN_FLOAT8(result);
 403 }
 404
 405 Datum
 406 float8up(PG_FUNCTION_ARGS)
 407 {
 408         float8          arg = PG_GETARG_FLOAT8(0);
 409
 410         PG_RETURN_FLOAT8(arg);
 411 }
 412
 413 Datum
 414 float8larger(PG_FUNCTION_ARGS)
 415 {
 416         float8          arg1 = PG_GETARG_FLOAT8(0);
 417         float8          arg2 = PG_GETARG_FLOAT8(1);
 418         float8          result;
 419
 420         result = ((arg1 > arg2) ? arg1 : arg2);
 421
 422         PG_RETURN_FLOAT8(result);
 423 }
 424
 425 Datum
 426 float8smaller(PG_FUNCTION_ARGS)
 427 {
 428         float8          arg1 = PG_GETARG_FLOAT8(0);
 429         float8          arg2 = PG_GETARG_FLOAT8(1);
 430         float8          result;
 431
 432         result = ((arg1 < arg2) ? arg1 : arg2);
 433
 434         PG_RETURN_FLOAT8(result);
 435 }
 436
 437
 438 /*
 439  *              ====================
 440  *              ARITHMETIC OPERATORS
 441  *              ====================
 442  */
 443
 444 /*
 445  *              float4pl                - returns arg1 + arg2
 446  *              float4mi                - returns arg1 - arg2
 447  *              float4mul               - returns arg1 * arg2
 448  *              float4div               - returns arg1 / arg2
 449  */
 450 Datum
 451 float4pl(PG_FUNCTION_ARGS)
 452 {
 453         float4          arg1 = PG_GETARG_FLOAT4(0);
 454         float4          arg2 = PG_GETARG_FLOAT4(1);
 455         double          result;
 456
 457         result = arg1 + arg2;
 458         CheckFloat4Val(result);
 459         PG_RETURN_FLOAT4((float4) result);
 460 }
 461
 462 Datum
 463 float4mi(PG_FUNCTION_ARGS)
 464 {
 465         float4          arg1 = PG_GETARG_FLOAT4(0);
 466         float4          arg2 = PG_GETARG_FLOAT4(1);
 467         double          result;
 468
 469         result = arg1 - arg2;
 470         CheckFloat4Val(result);
 471         PG_RETURN_FLOAT4((float4) result);
 472 }
 473
 474 Datum
 475 float4mul(PG_FUNCTION_ARGS)
 476 {
 477         float4          arg1 = PG_GETARG_FLOAT4(0);
 478         float4          arg2 = PG_GETARG_FLOAT4(1);
 479         double          result;
 480
 481         result = arg1 * arg2;
 482         CheckFloat4Val(result);
 483         PG_RETURN_FLOAT4((float4) result);
 484 }
 485
 486 Datum
 487 float4div(PG_FUNCTION_ARGS)
 488 {
 489         float4          arg1 = PG_GETARG_FLOAT4(0);
 490         float4          arg2 = PG_GETARG_FLOAT4(1);
 491         double          result;
 492
 493         if (arg2 == 0.0)
 494                 elog(ERROR, "float4div: divide by zero error");
 495
 496         /* Do division in float8, then check for overflow */
 497         result = (float8) arg1 / (float8) arg2;
 498
 499         CheckFloat4Val(result);
 500         PG_RETURN_FLOAT4((float4) result);
 501 }
 502
 503 /*
 504  *              float8pl                - returns arg1 + arg2
 505  *              float8mi                - returns arg1 - arg2
 506  *              float8mul               - returns arg1 * arg2
 507  *              float8div               - returns arg1 / arg2
 508  */
 509 Datum
 510 float8pl(PG_FUNCTION_ARGS)
 511 {
 512         float8          arg1 = PG_GETARG_FLOAT8(0);
 513         float8          arg2 = PG_GETARG_FLOAT8(1);
 514         float8          result;
 515
 516         result = arg1 + arg2;
 517
 518         CheckFloat8Val(result);
 519         PG_RETURN_FLOAT8(result);
 520 }
 521
 522 Datum
 523 float8mi(PG_FUNCTION_ARGS)
 524 {
 525         float8          arg1 = PG_GETARG_FLOAT8(0);
 526         float8          arg2 = PG_GETARG_FLOAT8(1);
 527         float8          result;
 528
 529         result = arg1 - arg2;
 530
 531         CheckFloat8Val(result);
 532         PG_RETURN_FLOAT8(result);
 533 }
 534
 535 Datum
 536 float8mul(PG_FUNCTION_ARGS)
 537 {
 538         float8          arg1 = PG_GETARG_FLOAT8(0);
 539         float8          arg2 = PG_GETARG_FLOAT8(1);
 540         float8          result;
 541
 542         result = arg1 * arg2;
 543
 544         CheckFloat8Val(result);
 545         PG_RETURN_FLOAT8(result);
 546 }
 547
 548 Datum
 549 float8div(PG_FUNCTION_ARGS)
 550 {
 551         float8          arg1 = PG_GETARG_FLOAT8(0);
 552         float8          arg2 = PG_GETARG_FLOAT8(1);
 553         float8          result;
 554
 555         if (arg2 == 0.0)
 556                 elog(ERROR, "float8div: divide by zero error");
 557
 558         result = arg1 / arg2;
 559
 560         CheckFloat8Val(result);
 561         PG_RETURN_FLOAT8(result);
 562 }
 563
 564
 565 /*
 566  *              ====================
 567  *              COMPARISON OPERATORS
 568  *              ====================
 569  */
 570
 571 /*
 572  *              float4{eq,ne,lt,le,gt,ge}               - float4/float4 comparison operations
 573  */
 574 static int
 575 float4_cmp_internal(float4 a, float4 b)
 576 {
 577         /*
 578          * We consider all NANs to be equal and larger than any non-NAN. This
 579          * is somewhat arbitrary; the important thing is to have a consistent
 580          * sort order.
 581          */
 582         if (isnan(a))
 583         {
 584                 if (isnan(b))
 585                         return 0;                       /* NAN = NAN */
 586                 else
 587                         return 1;                       /* NAN > non-NAN */
 588         }
 589         else if (isnan(b))
 590         {
 591                 return -1;                              /* non-NAN < NAN */
 592         }
 593         else
 594         {
 595                 if (a > b)
 596                         return 1;
 597                 else if (a < b)
 598                         return -1;
 599                 else
 600                         return 0;
 601         }
 602 }
 603
 604 Datum
 605 float4eq(PG_FUNCTION_ARGS)
 606 {
 607         float4          arg1 = PG_GETARG_FLOAT4(0);
 608         float4          arg2 = PG_GETARG_FLOAT4(1);
 609
 610         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) == 0);
 611 }
 612
 613 Datum
 614 float4ne(PG_FUNCTION_ARGS)
 615 {
 616         float4          arg1 = PG_GETARG_FLOAT4(0);
 617         float4          arg2 = PG_GETARG_FLOAT4(1);
 618
 619         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) != 0);
 620 }
 621
 622 Datum
 623 float4lt(PG_FUNCTION_ARGS)
 624 {
 625         float4          arg1 = PG_GETARG_FLOAT4(0);
 626         float4          arg2 = PG_GETARG_FLOAT4(1);
 627
 628         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) < 0);
 629 }
 630
 631 Datum
 632 float4le(PG_FUNCTION_ARGS)
 633 {
 634         float4          arg1 = PG_GETARG_FLOAT4(0);
 635         float4          arg2 = PG_GETARG_FLOAT4(1);
 636
 637         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) <= 0);
 638 }
 639
 640 Datum
 641 float4gt(PG_FUNCTION_ARGS)
 642 {
 643         float4          arg1 = PG_GETARG_FLOAT4(0);
 644         float4          arg2 = PG_GETARG_FLOAT4(1);
 645
 646         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) > 0);
 647 }
 648
 649 Datum
 650 float4ge(PG_FUNCTION_ARGS)
 651 {
 652         float4          arg1 = PG_GETARG_FLOAT4(0);
 653         float4          arg2 = PG_GETARG_FLOAT4(1);
 654
 655         PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) >= 0);
 656 }
 657
 658 Datum
 659 btfloat4cmp(PG_FUNCTION_ARGS)
 660 {
 661         float4          arg1 = PG_GETARG_FLOAT4(0);
 662         float4          arg2 = PG_GETARG_FLOAT4(1);
 663
 664         PG_RETURN_INT32(float4_cmp_internal(arg1, arg2));
 665 }
 666
 667 /*
 668  *              float8{eq,ne,lt,le,gt,ge}               - float8/float8 comparison operations
 669  */
 670 static int
 671 float8_cmp_internal(float8 a, float8 b)
 672 {
 673         /*
 674          * We consider all NANs to be equal and larger than any non-NAN. This
 675          * is somewhat arbitrary; the important thing is to have a consistent
 676          * sort order.
 677          */
 678         if (isnan(a))
 679         {
 680                 if (isnan(b))
 681                         return 0;                       /* NAN = NAN */
 682                 else
 683                         return 1;                       /* NAN > non-NAN */
 684         }
 685         else if (isnan(b))
 686         {
 687                 return -1;                              /* non-NAN < NAN */
 688         }
 689         else
 690         {
 691                 if (a > b)
 692                         return 1;
 693                 else if (a < b)
 694                         return -1;
 695                 else
 696                         return 0;
 697         }
 698 }
 699
 700 Datum
 701 float8eq(PG_FUNCTION_ARGS)
 702 {
 703         float8          arg1 = PG_GETARG_FLOAT8(0);
 704         float8          arg2 = PG_GETARG_FLOAT8(1);
 705
 706         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
 707 }
 708
 709 Datum
 710 float8ne(PG_FUNCTION_ARGS)
 711 {
 712         float8          arg1 = PG_GETARG_FLOAT8(0);
 713         float8          arg2 = PG_GETARG_FLOAT8(1);
 714
 715         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
 716 }
 717
 718 Datum
 719 float8lt(PG_FUNCTION_ARGS)
 720 {
 721         float8          arg1 = PG_GETARG_FLOAT8(0);
 722         float8          arg2 = PG_GETARG_FLOAT8(1);
 723
 724         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
 725 }
 726
 727 Datum
 728 float8le(PG_FUNCTION_ARGS)
 729 {
 730         float8          arg1 = PG_GETARG_FLOAT8(0);
 731         float8          arg2 = PG_GETARG_FLOAT8(1);
 732
 733         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
 734 }
 735
 736 Datum
 737 float8gt(PG_FUNCTION_ARGS)
 738 {
 739         float8          arg1 = PG_GETARG_FLOAT8(0);
 740         float8          arg2 = PG_GETARG_FLOAT8(1);
 741
 742         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
 743 }
 744
 745 Datum
 746 float8ge(PG_FUNCTION_ARGS)
 747 {
 748         float8          arg1 = PG_GETARG_FLOAT8(0);
 749         float8          arg2 = PG_GETARG_FLOAT8(1);
 750
 751         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
 752 }
 753
 754 Datum
 755 btfloat8cmp(PG_FUNCTION_ARGS)
 756 {
 757         float8          arg1 = PG_GETARG_FLOAT8(0);
 758         float8          arg2 = PG_GETARG_FLOAT8(1);
 759
 760         PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
 761 }
 762
 763
 764 /*
 765  *              ===================
 766  *              CONVERSION ROUTINES
 767  *              ===================
 768  */
 769
 770 /*
 771  *              ftod                    - converts a float4 number to a float8 number
 772  */
 773 Datum
 774 ftod(PG_FUNCTION_ARGS)
 775 {
 776         float4          num = PG_GETARG_FLOAT4(0);
 777
 778         PG_RETURN_FLOAT8((float8) num);
 779 }
 780
 781
 782 /*
 783  *              dtof                    - converts a float8 number to a float4 number
 784  */
 785 Datum
 786 dtof(PG_FUNCTION_ARGS)
 787 {
 788         float8          num = PG_GETARG_FLOAT8(0);
 789
 790         CheckFloat4Val(num);
 791
 792         PG_RETURN_FLOAT4((float4) num);
 793 }
 794
 795
 796 /*
 797  *              dtoi4                   - converts a float8 number to an int4 number
 798  */
 799 Datum
 800 dtoi4(PG_FUNCTION_ARGS)
 801 {
 802         float8          num = PG_GETARG_FLOAT8(0);
 803         int32           result;
 804
 805         if ((num < INT_MIN) || (num > INT_MAX))
 806                 elog(ERROR, "dtoi4: integer out of range");
 807
 808         result = (int32) rint(num);
 809         PG_RETURN_INT32(result);
 810 }
 811
 812
 813 /*
 814  *              dtoi2                   - converts a float8 number to an int2 number
 815  */
 816 Datum
 817 dtoi2(PG_FUNCTION_ARGS)
 818 {
 819         float8          num = PG_GETARG_FLOAT8(0);
 820         int16           result;
 821
 822         if ((num < SHRT_MIN) || (num > SHRT_MAX))
 823                 elog(ERROR, "dtoi2: integer out of range");
 824
 825         result = (int16) rint(num);
 826         PG_RETURN_INT16(result);
 827 }
 828
 829
 830 /*
 831  *              i4tod                   - converts an int4 number to a float8 number
 832  */
 833 Datum
 834 i4tod(PG_FUNCTION_ARGS)
 835 {
 836         int32           num = PG_GETARG_INT32(0);
 837         float8          result;
 838
 839         result = num;
 840         PG_RETURN_FLOAT8(result);
 841 }
 842
 843
 844 /*
 845  *              i2tod                   - converts an int2 number to a float8 number
 846  */
 847 Datum
 848 i2tod(PG_FUNCTION_ARGS)
 849 {
 850         int16           num = PG_GETARG_INT16(0);
 851         float8          result;
 852
 853         result = num;
 854         PG_RETURN_FLOAT8(result);
 855 }
 856
 857
 858 /*
 859  *              ftoi4                   - converts a float4 number to an int4 number
 860  */
 861 Datum
 862 ftoi4(PG_FUNCTION_ARGS)
 863 {
 864         float4          num = PG_GETARG_FLOAT4(0);
 865         int32           result;
 866
 867         if ((num < INT_MIN) || (num > INT_MAX))
 868                 elog(ERROR, "ftoi4: integer out of range");
 869
 870         result = (int32) rint(num);
 871         PG_RETURN_INT32(result);
 872 }
 873
 874
 875 /*
 876  *              ftoi2                   - converts a float4 number to an int2 number
 877  */
 878 Datum
 879 ftoi2(PG_FUNCTION_ARGS)
 880 {
 881         float4          num = PG_GETARG_FLOAT4(0);
 882         int16           result;
 883
 884         if ((num < SHRT_MIN) || (num > SHRT_MAX))
 885                 elog(ERROR, "ftoi2: integer out of range");
 886
 887         result = (int16) rint(num);
 888         PG_RETURN_INT16(result);
 889 }
 890
 891
 892 /*
 893  *              i4tof                   - converts an int4 number to a float8 number
 894  */
 895 Datum
 896 i4tof(PG_FUNCTION_ARGS)
 897 {
 898         int32           num = PG_GETARG_INT32(0);
 899         float4          result;
 900
 901         result = num;
 902         PG_RETURN_FLOAT4(result);
 903 }
 904
 905
 906 /*
 907  *              i2tof                   - converts an int2 number to a float4 number
 908  */
 909 Datum
 910 i2tof(PG_FUNCTION_ARGS)
 911 {
 912         int16           num = PG_GETARG_INT16(0);
 913         float4          result;
 914
 915         result = num;
 916         PG_RETURN_FLOAT4(result);
 917 }
 918
 919
 920 /*
 921  *              float8_text             - converts a float8 number to a text string
 922  */
 923 Datum
 924 float8_text(PG_FUNCTION_ARGS)
 925 {
 926         float8          num = PG_GETARG_FLOAT8(0);
 927         text       *result;
 928         int                     len;
 929         char       *str;
 930
 931         str = DatumGetCString(DirectFunctionCall1(float8out,
 932                                                                                           Float8GetDatum(num)));
 933
 934         len = strlen(str) + VARHDRSZ;
 935
 936         result = (text *) palloc(len);
 937
 938         VARATT_SIZEP(result) = len;
 939         memcpy(VARDATA(result), str, (len - VARHDRSZ));
 940
 941         pfree(str);
 942
 943         PG_RETURN_TEXT_P(result);
 944 }
 945
 946
 947 /*
 948  *              text_float8             - converts a text string to a float8 number
 949  */
 950 Datum
 951 text_float8(PG_FUNCTION_ARGS)
 952 {
 953         text       *string = PG_GETARG_TEXT_P(0);
 954         Datum           result;
 955         int                     len;
 956         char       *str;
 957
 958         len = (VARSIZE(string) - VARHDRSZ);
 959         str = palloc(len + 1);
 960         memcpy(str, VARDATA(string), len);
 961         *(str + len) = '\0';
 962
 963         result = DirectFunctionCall1(float8in, CStringGetDatum(str));
 964
 965         pfree(str);
 966
 967         PG_RETURN_DATUM(result);
 968 }
 969
 970
 971 /*
 972  *              float4_text             - converts a float4 number to a text string
 973  */
 974 Datum
 975 float4_text(PG_FUNCTION_ARGS)
 976 {
 977         float4          num = PG_GETARG_FLOAT4(0);
 978         text       *result;
 979         int                     len;
 980         char       *str;
 981
 982         str = DatumGetCString(DirectFunctionCall1(float4out,
 983                                                                                           Float4GetDatum(num)));
 984
 985         len = strlen(str) + VARHDRSZ;
 986
 987         result = (text *) palloc(len);
 988
 989         VARATT_SIZEP(result) = len;
 990         memcpy(VARDATA(result), str, (len - VARHDRSZ));
 991
 992         pfree(str);
 993
 994         PG_RETURN_TEXT_P(result);
 995 }
 996
 997
 998 /*
 999  *              text_float4             - converts a text string to a float4 number
1000  */
1001 Datum
1002 text_float4(PG_FUNCTION_ARGS)
1003 {
1004         text       *string = PG_GETARG_TEXT_P(0);
1005         Datum           result;
1006         int                     len;
1007         char       *str;
1008
1009         len = (VARSIZE(string) - VARHDRSZ);
1010         str = palloc(len + 1);
1011         memcpy(str, VARDATA(string), len);
1012         *(str + len) = '\0';
1013
1014         result = DirectFunctionCall1(float4in, CStringGetDatum(str));
1015
1016         pfree(str);
1017
1018         PG_RETURN_DATUM(result);
1019 }
1020
1021
1022 /*
1023  *              =======================
1024  *              RANDOM FLOAT8 OPERATORS
1025  *              =======================
1026  */
1027
1028 /*
1029  *              dround                  - returns       ROUND(arg1)
1030  */
1031 Datum
1032 dround(PG_FUNCTION_ARGS)
1033 {
1034         float8          arg1 = PG_GETARG_FLOAT8(0);
1035         float8          result;
1036
1037         result = rint(arg1);
1038
1039         PG_RETURN_FLOAT8(result);
1040 }
1041
1042
1043 /*
1044  *              dtrunc                  - returns truncation-towards-zero of arg1,
1045  *                                                arg1 >= 0 ... the greatest integer less
1046  *                                                                              than or equal to arg1
1047  *                                                arg1 < 0      ... the least integer greater
1048  *                                                                              than or equal to arg1
1049  */
1050 Datum
1051 dtrunc(PG_FUNCTION_ARGS)
1052 {
1053         float8          arg1 = PG_GETARG_FLOAT8(0);
1054         float8          result;
1055
1056         if (arg1 >= 0)
1057                 result = floor(arg1);
1058         else
1059                 result = -floor(-arg1);
1060
1061         PG_RETURN_FLOAT8(result);
1062 }
1063
1064
1065 /*
1066  *              dsqrt                   - returns square root of arg1
1067  */
1068 Datum
1069 dsqrt(PG_FUNCTION_ARGS)
1070 {
1071         float8          arg1 = PG_GETARG_FLOAT8(0);
1072         float8          result;
1073
1074         if (arg1 < 0)
1075                 elog(ERROR, "can't take sqrt of a negative number");
1076
1077         result = sqrt(arg1);
1078
1079         CheckFloat8Val(result);
1080         PG_RETURN_FLOAT8(result);
1081 }
1082
1083
1084 /*
1085  *              dcbrt                   - returns cube root of arg1
1086  */
1087 Datum
1088 dcbrt(PG_FUNCTION_ARGS)
1089 {
1090         float8          arg1 = PG_GETARG_FLOAT8(0);
1091         float8          result;
1092
1093         result = cbrt(arg1);
1094         PG_RETURN_FLOAT8(result);
1095 }
1096
1097
1098 /*
1099  *              dpow                    - returns pow(arg1,arg2)
1100  */
1101 Datum
1102 dpow(PG_FUNCTION_ARGS)
1103 {
1104         float8          arg1 = PG_GETARG_FLOAT8(0);
1105         float8          arg2 = PG_GETARG_FLOAT8(1);
1106         float8          result;
1107
1108         /*
1109          * We must check both for errno getting set and for a NaN result, in
1110          * order to deal with the vagaries of different platforms...
1111          */
1112         errno = 0;
1113         result = pow(arg1, arg2);
1114         if (errno != 0
1115 #ifdef HAVE_FINITE
1116                 || !finite(result)
1117 #endif
1118                 )
1119                 elog(ERROR, "pow() result is out of range");
1120
1121         CheckFloat8Val(result);
1122         PG_RETURN_FLOAT8(result);
1123 }
1124
1125
1126 /*
1127  *              dexp                    - returns the exponential function of arg1
1128  */
1129 Datum
1130 dexp(PG_FUNCTION_ARGS)
1131 {
1132         float8          arg1 = PG_GETARG_FLOAT8(0);
1133         float8          result;
1134
1135         /*
1136          * We must check both for errno getting set and for a NaN result, in
1137          * order to deal with the vagaries of different platforms. Also, a
1138          * zero result implies unreported underflow.
1139          */
1140         errno = 0;
1141         result = exp(arg1);
1142         if (errno != 0 || result == 0.0
1143 #ifdef HAVE_FINITE
1144                 || !finite(result)
1145 #endif
1146                 )
1147                 elog(ERROR, "exp() result is out of range");
1148
1149         CheckFloat8Val(result);
1150         PG_RETURN_FLOAT8(result);
1151 }
1152
1153
1154 /*
1155  *              dlog1                   - returns the natural logarithm of arg1
1156  *                                                ("dlog" is already a logging routine...)
1157  */
1158 Datum
1159 dlog1(PG_FUNCTION_ARGS)
1160 {
1161         float8          arg1 = PG_GETARG_FLOAT8(0);
1162         float8          result;
1163
1164         if (arg1 == 0.0)
1165                 elog(ERROR, "can't take log of zero");
1166         if (arg1 < 0)
1167                 elog(ERROR, "can't take log of a negative number");
1168
1169         result = log(arg1);
1170
1171         CheckFloat8Val(result);
1172         PG_RETURN_FLOAT8(result);
1173 }
1174
1175
1176 /*
1177  *              dlog10                  - returns the base 10 logarithm of arg1
1178  */
1179 Datum
1180 dlog10(PG_FUNCTION_ARGS)
1181 {
1182         float8          arg1 = PG_GETARG_FLOAT8(0);
1183         float8          result;
1184
1185         if (arg1 == 0.0)
1186                 elog(ERROR, "can't take log of zero");
1187         if (arg1 < 0)
1188                 elog(ERROR, "can't take log of a negative number");
1189
1190         result = log10(arg1);
1191
1192         CheckFloat8Val(result);
1193         PG_RETURN_FLOAT8(result);
1194 }
1195
1196
1197 /*
1198  *              dacos                   - returns the arccos of arg1 (radians)
1199  */
1200 Datum
1201 dacos(PG_FUNCTION_ARGS)
1202 {
1203         float8          arg1 = PG_GETARG_FLOAT8(0);
1204         float8          result;
1205
1206         errno = 0;
1207         result = acos(arg1);
1208         if (errno != 0
1209 #ifdef HAVE_FINITE
1210                 || !finite(result)
1211 #endif
1212                 )
1213                 elog(ERROR, "acos(%f) input is out of range", arg1);
1214
1215         CheckFloat8Val(result);
1216         PG_RETURN_FLOAT8(result);
1217 }
1218
1219
1220 /*
1221  *              dasin                   - returns the arcsin of arg1 (radians)
1222  */
1223 Datum
1224 dasin(PG_FUNCTION_ARGS)
1225 {
1226         float8          arg1 = PG_GETARG_FLOAT8(0);
1227         float8          result;
1228
1229         errno = 0;
1230         result = asin(arg1);
1231         if (errno != 0
1232 #ifdef HAVE_FINITE
1233                 || !finite(result)
1234 #endif
1235                 )
1236                 elog(ERROR, "asin(%f) input is out of range", arg1);
1237
1238         CheckFloat8Val(result);
1239         PG_RETURN_FLOAT8(result);
1240 }
1241
1242
1243 /*
1244  *              datan                   - returns the arctan of arg1 (radians)
1245  */
1246 Datum
1247 datan(PG_FUNCTION_ARGS)
1248 {
1249         float8          arg1 = PG_GETARG_FLOAT8(0);
1250         float8          result;
1251
1252         errno = 0;
1253         result = atan(arg1);
1254         if (errno != 0
1255 #ifdef HAVE_FINITE
1256                 || !finite(result)
1257 #endif
1258                 )
1259                 elog(ERROR, "atan(%f) input is out of range", arg1);
1260
1261         CheckFloat8Val(result);
1262         PG_RETURN_FLOAT8(result);
1263 }
1264
1265
1266 /*
1267  *              atan2                   - returns the arctan2 of arg1 (radians)
1268  */
1269 Datum
1270 datan2(PG_FUNCTION_ARGS)
1271 {
1272         float8          arg1 = PG_GETARG_FLOAT8(0);
1273         float8          arg2 = PG_GETARG_FLOAT8(1);
1274         float8          result;
1275
1276         errno = 0;
1277         result = atan2(arg1, arg2);
1278         if (errno != 0
1279 #ifdef HAVE_FINITE
1280                 || !finite(result)
1281 #endif
1282                 )
1283                 elog(ERROR, "atan2(%f,%f) input is out of range", arg1, arg2);
1284
1285         CheckFloat8Val(result);
1286         PG_RETURN_FLOAT8(result);
1287 }
1288
1289
1290 /*
1291  *              dcos                    - returns the cosine of arg1 (radians)
1292  */
1293 Datum
1294 dcos(PG_FUNCTION_ARGS)
1295 {
1296         float8          arg1 = PG_GETARG_FLOAT8(0);
1297         float8          result;
1298
1299         errno = 0;
1300         result = cos(arg1);
1301         if (errno != 0
1302 #ifdef HAVE_FINITE
1303                 || !finite(result)
1304 #endif
1305                 )
1306                 elog(ERROR, "cos(%f) input is out of range", arg1);
1307
1308         CheckFloat8Val(result);
1309         PG_RETURN_FLOAT8(result);
1310 }
1311
1312
1313 /*
1314  *              dcot                    - returns the cotangent of arg1 (radians)
1315  */
1316 Datum
1317 dcot(PG_FUNCTION_ARGS)
1318 {
1319         float8          arg1 = PG_GETARG_FLOAT8(0);
1320         float8          result;
1321
1322         errno = 0;
1323         result = tan(arg1);
1324         if (errno != 0 || result == 0.0
1325 #ifdef HAVE_FINITE
1326                 || !finite(result)
1327 #endif
1328                 )
1329                 elog(ERROR, "cot(%f) input is out of range", arg1);
1330
1331         result = 1.0 / result;
1332         CheckFloat8Val(result);
1333         PG_RETURN_FLOAT8(result);
1334 }
1335
1336
1337 /*
1338  *              dsin                    - returns the sine of arg1 (radians)
1339  */
1340 Datum
1341 dsin(PG_FUNCTION_ARGS)
1342 {
1343         float8          arg1 = PG_GETARG_FLOAT8(0);
1344         float8          result;
1345
1346         errno = 0;
1347         result = sin(arg1);
1348         if (errno != 0
1349 #ifdef HAVE_FINITE
1350                 || !finite(result)
1351 #endif
1352                 )
1353                 elog(ERROR, "sin(%f) input is out of range", arg1);
1354
1355         CheckFloat8Val(result);
1356         PG_RETURN_FLOAT8(result);
1357 }
1358
1359
1360 /*
1361  *              dtan                    - returns the tangent of arg1 (radians)
1362  */
1363 Datum
1364 dtan(PG_FUNCTION_ARGS)
1365 {
1366         float8          arg1 = PG_GETARG_FLOAT8(0);
1367         float8          result;
1368
1369         errno = 0;
1370         result = tan(arg1);
1371         if (errno != 0
1372 #ifdef HAVE_FINITE
1373                 || !finite(result)
1374 #endif
1375                 )
1376                 elog(ERROR, "tan(%f) input is out of range", arg1);
1377
1378         CheckFloat8Val(result);
1379         PG_RETURN_FLOAT8(result);
1380 }
1381
1382
1383 /*
1384  *              degrees         - returns degrees converted from radians
1385  */
1386 Datum
1387 degrees(PG_FUNCTION_ARGS)
1388 {
1389         float8          arg1 = PG_GETARG_FLOAT8(0);
1390         float8          result;
1391
1392         result = arg1 * (180.0 / M_PI);
1393
1394         CheckFloat8Val(result);
1395         PG_RETURN_FLOAT8(result);
1396 }
1397
1398
1399 /*
1400  *              dpi                             - returns the constant PI
1401  */
1402 Datum
1403 dpi(PG_FUNCTION_ARGS)
1404 {
1405         PG_RETURN_FLOAT8(M_PI);
1406 }
1407
1408
1409 /*
1410  *              radians         - returns radians converted from degrees
1411  */
1412 Datum
1413 radians(PG_FUNCTION_ARGS)
1414 {
1415         float8          arg1 = PG_GETARG_FLOAT8(0);
1416         float8          result;
1417
1418         result = arg1 * (M_PI / 180.0);
1419
1420         CheckFloat8Val(result);
1421         PG_RETURN_FLOAT8(result);
1422 }
1423
1424
1425 /*
1426  *              drandom         - returns a random number
1427  */
1428 Datum
1429 drandom(PG_FUNCTION_ARGS)
1430 {
1431         float8          result;
1432
1433         /* result 0.0-1.0 */
1434         result = ((double) random()) / ((double) MAX_RANDOM_VALUE);
1435
1436         PG_RETURN_FLOAT8(result);
1437 }
1438
1439
1440 /*
1441  *              setseed         - set seed for the random number generator
1442  */
1443 Datum
1444 setseed(PG_FUNCTION_ARGS)
1445 {
1446         float8          seed = PG_GETARG_FLOAT8(0);
1447         int                     iseed = (int) (seed * MAX_RANDOM_VALUE);
1448
1449         srandom((unsigned int) iseed);
1450
1451         PG_RETURN_INT32(iseed);
1452 }
1453
1454
1455
1456 /*
1457  *              =========================
1458  *              FLOAT AGGREGATE OPERATORS
1459  *              =========================
1460  *
1461  *              float8_accum    - accumulate for AVG(), STDDEV(), etc
1462  *              float4_accum    - same, but input data is float4
1463  *              float8_avg              - produce final result for float AVG()
1464  *              float8_variance - produce final result for float VARIANCE()
1465  *              float8_stddev   - produce final result for float STDDEV()
1466  *
1467  * The transition datatype for all these aggregates is a 3-element array
1468  * of float8, holding the values N, sum(X), sum(X*X) in that order.
1469  *
1470  * Note that we represent N as a float to avoid having to build a special
1471  * datatype.  Given a reasonable floating-point implementation, there should
1472  * be no accuracy loss unless N exceeds 2 ^ 52 or so (by which time the
1473  * user will have doubtless lost interest anyway...)
1474  */
1475
1476 static float8 *
1477 check_float8_array(ArrayType *transarray, const char *caller)
1478 {
1479         /*
1480          * We expect the input to be a 3-element float array; verify that. We
1481          * don't need to use deconstruct_array() since the array data is just
1482          * going to look like a C array of 3 float8 values.
1483          */
1484         if (ARR_NDIM(transarray) != 1 ||
1485                 ARR_DIMS(transarray)[0] != 3 ||
1486                 ARR_ELEMTYPE(transarray) != FLOAT8OID)
1487                 elog(ERROR, "%s: expected 3-element float8 array", caller);
1488         return (float8 *) ARR_DATA_PTR(transarray);
1489 }
1490
1491 Datum
1492 float8_accum(PG_FUNCTION_ARGS)
1493 {
1494         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
1495         float8          newval = PG_GETARG_FLOAT8(1);
1496         float8     *transvalues;
1497         float8          N,
1498                                 sumX,
1499                                 sumX2;
1500         Datum           transdatums[3];
1501         ArrayType  *result;
1502
1503         transvalues = check_float8_array(transarray, "float8_accum");
1504         N = transvalues[0];
1505         sumX = transvalues[1];
1506         sumX2 = transvalues[2];
1507
1508         N += 1.0;
1509         sumX += newval;
1510         sumX2 += newval * newval;
1511
1512         transdatums[0] = Float8GetDatumFast(N);
1513         transdatums[1] = Float8GetDatumFast(sumX);
1514         transdatums[2] = Float8GetDatumFast(sumX2);
1515
1516         result = construct_array(transdatums, 3,
1517                                                          FLOAT8OID,
1518                                                  sizeof(float8), false /* float8 byval */ , 'd');
1519
1520         PG_RETURN_ARRAYTYPE_P(result);
1521 }
1522
1523 Datum
1524 float4_accum(PG_FUNCTION_ARGS)
1525 {
1526         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
1527         float4          newval4 = PG_GETARG_FLOAT4(1);
1528         float8     *transvalues;
1529         float8          N,
1530                                 sumX,
1531                                 sumX2,
1532                                 newval;
1533         Datum           transdatums[3];
1534         ArrayType  *result;
1535
1536         transvalues = check_float8_array(transarray, "float4_accum");
1537         N = transvalues[0];
1538         sumX = transvalues[1];
1539         sumX2 = transvalues[2];
1540
1541         /* Do arithmetic in float8 for best accuracy */
1542         newval = newval4;
1543
1544         N += 1.0;
1545         sumX += newval;
1546         sumX2 += newval * newval;
1547
1548         transdatums[0] = Float8GetDatumFast(N);
1549         transdatums[1] = Float8GetDatumFast(sumX);
1550         transdatums[2] = Float8GetDatumFast(sumX2);
1551
1552         result = construct_array(transdatums, 3,
1553                                                          FLOAT8OID,
1554                                                  sizeof(float8), false /* float8 byval */ , 'd');
1555
1556         PG_RETURN_ARRAYTYPE_P(result);
1557 }
1558
1559 Datum
1560 float8_avg(PG_FUNCTION_ARGS)
1561 {
1562         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
1563         float8     *transvalues;
1564         float8          N,
1565                                 sumX;
1566
1567         transvalues = check_float8_array(transarray, "float8_avg");
1568         N = transvalues[0];
1569         sumX = transvalues[1];
1570         /* ignore sumX2 */
1571
1572         /* SQL92 defines AVG of no values to be NULL */
1573         if (N == 0.0)
1574                 PG_RETURN_NULL();
1575
1576         PG_RETURN_FLOAT8(sumX / N);
1577 }
1578
1579 Datum
1580 float8_variance(PG_FUNCTION_ARGS)
1581 {
1582         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
1583         float8     *transvalues;
1584         float8          N,
1585                                 sumX,
1586                                 sumX2,
1587                                 numerator;
1588
1589         transvalues = check_float8_array(transarray, "float8_variance");
1590         N = transvalues[0];
1591         sumX = transvalues[1];
1592         sumX2 = transvalues[2];
1593
1594         /* We define VARIANCE of no values to be NULL, of 1 value to be 0 */
1595         if (N == 0.0)
1596                 PG_RETURN_NULL();
1597
1598         if (N <= 1.0)
1599                 PG_RETURN_FLOAT8(0.0);
1600
1601         numerator = N * sumX2 - sumX * sumX;
1602
1603         /* Watch out for roundoff error producing a negative numerator */
1604         if (numerator <= 0.0)
1605                 PG_RETURN_FLOAT8(0.0);
1606
1607         PG_RETURN_FLOAT8(numerator / (N * (N - 1.0)));
1608 }
1609
1610 Datum
1611 float8_stddev(PG_FUNCTION_ARGS)
1612 {
1613         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
1614         float8     *transvalues;
1615         float8          N,
1616                                 sumX,
1617                                 sumX2,
1618                                 numerator;
1619
1620         transvalues = check_float8_array(transarray, "float8_stddev");
1621         N = transvalues[0];
1622         sumX = transvalues[1];
1623         sumX2 = transvalues[2];
1624
1625         /* We define STDDEV of no values to be NULL, of 1 value to be 0 */
1626         if (N == 0.0)
1627                 PG_RETURN_NULL();
1628
1629         if (N <= 1.0)
1630                 PG_RETURN_FLOAT8(0.0);
1631
1632         numerator = N * sumX2 - sumX * sumX;
1633
1634         /* Watch out for roundoff error producing a negative numerator */
1635         if (numerator <= 0.0)
1636                 PG_RETURN_FLOAT8(0.0);
1637
1638         PG_RETURN_FLOAT8(sqrt(numerator / (N * (N - 1.0))));
1639 }
1640
1641
1642 /*
1643  *              ====================================
1644  *              MIXED-PRECISION ARITHMETIC OPERATORS
1645  *              ====================================
1646  */
1647
1648 /*
1649  *              float48pl               - returns arg1 + arg2
1650  *              float48mi               - returns arg1 - arg2
1651  *              float48mul              - returns arg1 * arg2
1652  *              float48div              - returns arg1 / arg2
1653  */
1654 Datum
1655 float48pl(PG_FUNCTION_ARGS)
1656 {
1657         float4          arg1 = PG_GETARG_FLOAT4(0);
1658         float8          arg2 = PG_GETARG_FLOAT8(1);
1659         float8          result;
1660
1661         result = arg1 + arg2;
1662         CheckFloat8Val(result);
1663         PG_RETURN_FLOAT8(result);
1664 }
1665
1666 Datum
1667 float48mi(PG_FUNCTION_ARGS)
1668 {
1669         float4          arg1 = PG_GETARG_FLOAT4(0);
1670         float8          arg2 = PG_GETARG_FLOAT8(1);
1671         float8          result;
1672
1673         result = arg1 - arg2;
1674         CheckFloat8Val(result);
1675         PG_RETURN_FLOAT8(result);
1676 }
1677
1678 Datum
1679 float48mul(PG_FUNCTION_ARGS)
1680 {
1681         float4          arg1 = PG_GETARG_FLOAT4(0);
1682         float8          arg2 = PG_GETARG_FLOAT8(1);
1683         float8          result;
1684
1685         result = arg1 * arg2;
1686         CheckFloat8Val(result);
1687         PG_RETURN_FLOAT8(result);
1688 }
1689
1690 Datum
1691 float48div(PG_FUNCTION_ARGS)
1692 {
1693         float4          arg1 = PG_GETARG_FLOAT4(0);
1694         float8          arg2 = PG_GETARG_FLOAT8(1);
1695         float8          result;
1696
1697         if (arg2 == 0.0)
1698                 elog(ERROR, "float48div: divide by zero");
1699
1700         result = arg1 / arg2;
1701         CheckFloat8Val(result);
1702         PG_RETURN_FLOAT8(result);
1703 }
1704
1705 /*
1706  *              float84pl               - returns arg1 + arg2
1707  *              float84mi               - returns arg1 - arg2
1708  *              float84mul              - returns arg1 * arg2
1709  *              float84div              - returns arg1 / arg2
1710  */
1711 Datum
1712 float84pl(PG_FUNCTION_ARGS)
1713 {
1714         float8          arg1 = PG_GETARG_FLOAT8(0);
1715         float4          arg2 = PG_GETARG_FLOAT4(1);
1716         float8          result;
1717
1718         result = arg1 + arg2;
1719
1720         CheckFloat8Val(result);
1721         PG_RETURN_FLOAT8(result);
1722 }
1723
1724 Datum
1725 float84mi(PG_FUNCTION_ARGS)
1726 {
1727         float8          arg1 = PG_GETARG_FLOAT8(0);
1728         float4          arg2 = PG_GETARG_FLOAT4(1);
1729         float8          result;
1730
1731         result = arg1 - arg2;
1732
1733         CheckFloat8Val(result);
1734         PG_RETURN_FLOAT8(result);
1735 }
1736
1737 Datum
1738 float84mul(PG_FUNCTION_ARGS)
1739 {
1740         float8          arg1 = PG_GETARG_FLOAT8(0);
1741         float4          arg2 = PG_GETARG_FLOAT4(1);
1742         float8          result;
1743
1744         result = arg1 * arg2;
1745
1746         CheckFloat8Val(result);
1747         PG_RETURN_FLOAT8(result);
1748 }
1749
1750 Datum
1751 float84div(PG_FUNCTION_ARGS)
1752 {
1753         float8          arg1 = PG_GETARG_FLOAT8(0);
1754         float4          arg2 = PG_GETARG_FLOAT4(1);
1755         float8          result;
1756
1757         if (arg2 == 0.0)
1758                 elog(ERROR, "float84div: divide by zero");
1759
1760         result = arg1 / arg2;
1761
1762         CheckFloat8Val(result);
1763         PG_RETURN_FLOAT8(result);
1764 }
1765
1766 /*
1767  *              ====================
1768  *              COMPARISON OPERATORS
1769  *              ====================
1770  */
1771
1772 /*
1773  *              float48{eq,ne,lt,le,gt,ge}              - float4/float8 comparison operations
1774  */
1775 Datum
1776 float48eq(PG_FUNCTION_ARGS)
1777 {
1778         float4          arg1 = PG_GETARG_FLOAT4(0);
1779         float8          arg2 = PG_GETARG_FLOAT8(1);
1780
1781         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
1782 }
1783
1784 Datum
1785 float48ne(PG_FUNCTION_ARGS)
1786 {
1787         float4          arg1 = PG_GETARG_FLOAT4(0);
1788         float8          arg2 = PG_GETARG_FLOAT8(1);
1789
1790         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
1791 }
1792
1793 Datum
1794 float48lt(PG_FUNCTION_ARGS)
1795 {
1796         float4          arg1 = PG_GETARG_FLOAT4(0);
1797         float8          arg2 = PG_GETARG_FLOAT8(1);
1798
1799         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
1800 }
1801
1802 Datum
1803 float48le(PG_FUNCTION_ARGS)
1804 {
1805         float4          arg1 = PG_GETARG_FLOAT4(0);
1806         float8          arg2 = PG_GETARG_FLOAT8(1);
1807
1808         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
1809 }
1810
1811 Datum
1812 float48gt(PG_FUNCTION_ARGS)
1813 {
1814         float4          arg1 = PG_GETARG_FLOAT4(0);
1815         float8          arg2 = PG_GETARG_FLOAT8(1);
1816
1817         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
1818 }
1819
1820 Datum
1821 float48ge(PG_FUNCTION_ARGS)
1822 {
1823         float4          arg1 = PG_GETARG_FLOAT4(0);
1824         float8          arg2 = PG_GETARG_FLOAT8(1);
1825
1826         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
1827 }
1828
1829 /*
1830  *              float84{eq,ne,lt,le,gt,ge}              - float8/float4 comparison operations
1831  */
1832 Datum
1833 float84eq(PG_FUNCTION_ARGS)
1834 {
1835         float8          arg1 = PG_GETARG_FLOAT8(0);
1836         float4          arg2 = PG_GETARG_FLOAT4(1);
1837
1838         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
1839 }
1840
1841 Datum
1842 float84ne(PG_FUNCTION_ARGS)
1843 {
1844         float8          arg1 = PG_GETARG_FLOAT8(0);
1845         float4          arg2 = PG_GETARG_FLOAT4(1);
1846
1847         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
1848 }
1849
1850 Datum
1851 float84lt(PG_FUNCTION_ARGS)
1852 {
1853         float8          arg1 = PG_GETARG_FLOAT8(0);
1854         float4          arg2 = PG_GETARG_FLOAT4(1);
1855
1856         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
1857 }
1858
1859 Datum
1860 float84le(PG_FUNCTION_ARGS)
1861 {
1862         float8          arg1 = PG_GETARG_FLOAT8(0);
1863         float4          arg2 = PG_GETARG_FLOAT4(1);
1864
1865         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
1866 }
1867
1868 Datum
1869 float84gt(PG_FUNCTION_ARGS)
1870 {
1871         float8          arg1 = PG_GETARG_FLOAT8(0);
1872         float4          arg2 = PG_GETARG_FLOAT4(1);
1873
1874         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
1875 }
1876
1877 Datum
1878 float84ge(PG_FUNCTION_ARGS)
1879 {
1880         float8          arg1 = PG_GETARG_FLOAT8(0);
1881         float4          arg2 = PG_GETARG_FLOAT4(1);
1882
1883         PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
1884 }
1885
1886 /* ========== PRIVATE ROUTINES ========== */
1887
1888 /* From "fdlibm" @ netlib.att.com */
1889
1890 #ifndef HAVE_RINT
1891
1892 /* @(#)s_rint.c 5.1 93/09/24 */
1893 /*
1894  * ====================================================
1895  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
1896  *
1897  * Developed at SunPro, a Sun Microsystems, Inc. business.
1898  * Permission to use, copy, modify, and distribute this
1899  * software is freely granted, provided that this notice
1900  * is preserved.
1901  * ====================================================
1902  */
1903
1904 /*
1905  * rint(x)
1906  * Return x rounded to integral value according to the prevailing
1907  * rounding mode.
1908  * Method:
1909  *              Using floating addition.
1910  * Exception:
1911  *              Inexact flag raised if x not equal to rint(x).
1912  */
1913
1914 static const double one = 1.0,
1915                         TWO52[2] = {
1916         4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
1917         -4.50359962737049600000e+15,    /* 0xC3300000, 0x00000000 */
1918 };
1919
1920 static double
1921 rint(double x)
1922 {
1923         int                     i0,
1924                                 n0,
1925                                 j0,
1926                                 sx;
1927         unsigned        i,
1928                                 i1;
1929         double          w,
1930                                 t;
1931
1932         n0 = (*((int *) &one) >> 29) ^ 1;
1933         i0 = *(n0 + (int *) &x);
1934         sx = (i0 >> 31) & 1;
1935         i1 = *(1 - n0 + (int *) &x);
1936         j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
1937         if (j0 < 20)
1938         {
1939                 if (j0 < 0)
1940                 {
1941                         if (((i0 & 0x7fffffff) | i1) == 0)
1942                                 return x;
1943                         i1 |= (i0 & 0x0fffff);
1944                         i0 &= 0xfffe0000;
1945                         i0 |= ((i1 | -i1) >> 12) & 0x80000;
1946                         *(n0 + (int *) &x) = i0;
1947                         w = TWO52[sx] + x;
1948                         t = w - TWO52[sx];
1949                         i0 = *(n0 + (int *) &t);
1950                         *(n0 + (int *) &t) = (i0 & 0x7fffffff) | (sx << 31);
1951                         return t;
1952                 }
1953                 else
1954                 {
1955                         i = (0x000fffff) >> j0;
1956                         if (((i0 & i) | i1) == 0)
1957                                 return x;               /* x is integral */
1958                         i >>= 1;
1959                         if (((i0 & i) | i1) != 0)
1960                         {
1961                                 if (j0 == 19)
1962                                         i1 = 0x40000000;
1963                                 else
1964                                         i0 = (i0 & (~i)) | ((0x20000) >> j0);
1965                         }
1966                 }
1967         }
1968         else if (j0 > 51)
1969         {
1970                 if (j0 == 0x400)
1971                         return x + x;           /* inf or NaN */
1972                 else
1973                         return x;                       /* x is integral */
1974         }
1975         else
1976         {
1977                 i = ((unsigned) (0xffffffff)) >> (j0 - 20);
1978                 if ((i1 & i) == 0)
1979                         return x;                       /* x is integral */
1980                 i >>= 1;
1981                 if ((i1 & i) != 0)
1982                         i1 = (i1 & (~i)) | ((0x40000000) >> (j0 - 20));
1983         }
1984         *(n0 + (int *) &x) = i0;
1985         *(1 - n0 + (int *) &x) = i1;
1986         w = TWO52[sx] + x;
1987         return w - TWO52[sx];
1988 }
1989 #endif   /* !HAVE_RINT */
1990
1991 #ifndef HAVE_CBRT
1992
1993 static double
1994 cbrt(double x)
1995 {
1996         int                     isneg = (x < 0.0);
1997         double          tmpres = pow(fabs(x), (double) 1.0 / (double) 3.0);
1998
1999         return isneg ? -tmpres : tmpres;
2000 }
2001
2002 #endif   /* !HAVE_CBRT */