1 /*------------------------------------------------------------------------- 2 * 3 * float.c 4 * Functions for the built-in floating-point types. 5 * 6 * Portions Copyright (c) 1996-2018, PostgreSQL Global Development Group 7 * Portions Copyright (c) 1994, Regents of the University of California 8 * 9 * 10 * IDENTIFICATION 11 * src/backend/utils/adt/float.c 12 * 13 *------------------------------------------------------------------------- 14 */ 15 #include "postgres.h" 16 17 #include <ctype.h> 18 #include <float.h> 19 #include <math.h> 20 #include <limits.h> 21 22 #include "catalog/pg_type.h" 23 #include "common/int.h" 24 #include "libpq/pqformat.h" 25 #include "utils/array.h" 26 #include "utils/builtins.h" 27 #include "utils/sortsupport.h" 28 29 30 #ifndef M_PI 31 /* from my RH5.2 gcc math.h file - thomas 2000-04-03 */ 32 #define M_PI 3.14159265358979323846 33 #endif 34 35 /* Radians per degree, a.k.a. PI / 180 */ 36 #define RADIANS_PER_DEGREE 0.0174532925199432957692 37 38 /* Visual C++ etc lacks NAN, and won't accept 0.0/0.0. NAN definition from 39 * http://msdn.microsoft.com/library/default.asp?url=/library/en-us/vclang/html/vclrfNotNumberNANItems.asp 40 */ 41 #if defined(WIN32) && !defined(NAN) 42 static const uint32 nan[2] = {0xffffffff, 0x7fffffff}; 43 44 #define NAN (*(const double *) nan) 45 #endif 46 47 /* 48 * check to see if a float4/8 val has underflowed or overflowed 49 */ 50 #define CHECKFLOATVAL(val, inf_is_valid, zero_is_valid) \ 51 do { \ 52 if (isinf(val) && !(inf_is_valid)) \ 53 ereport(ERROR, \ 54 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), \ 55 errmsg("value out of range: overflow"))); \ 56 \ 57 if ((val) == 0.0 && !(zero_is_valid)) \ 58 ereport(ERROR, \ 59 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), \ 60 errmsg("value out of range: underflow"))); \ 61 } while(0) 62 63 64 /* Configurable GUC parameter */ 65 int extra_float_digits = 0; /* Added to DBL_DIG or FLT_DIG */ 66 67 /* Cached constants for degree-based trig functions */ 68 static bool degree_consts_set = false; 69 static float8 sin_30 = 0; 70 static float8 one_minus_cos_60 = 0; 71 static float8 asin_0_5 = 0; 72 static float8 acos_0_5 = 0; 73 static float8 atan_1_0 = 0; 74 static float8 tan_45 = 0; 75 static float8 cot_45 = 0; 76 77 /* 78 * These are intentionally not static; don't "fix" them. They will never 79 * be referenced by other files, much less changed; but we don't want the 80 * compiler to know that, else it might try to precompute expressions 81 * involving them. See comments for init_degree_constants(). 82 */ 83 float8 degree_c_thirty = 30.0; 84 float8 degree_c_forty_five = 45.0; 85 float8 degree_c_sixty = 60.0; 86 float8 degree_c_one_half = 0.5; 87 float8 degree_c_one = 1.0; 88 89 /* Local function prototypes */ 90 static double sind_q1(double x); 91 static double cosd_q1(double x); 92 static void init_degree_constants(void); 93 94 #ifndef HAVE_CBRT 95 /* 96 * Some machines (in particular, some versions of AIX) have an extern 97 * declaration for cbrt() in <math.h> but fail to provide the actual 98 * function, which causes configure to not set HAVE_CBRT. Furthermore, 99 * their compilers spit up at the mismatch between extern declaration 100 * and static definition. We work around that here by the expedient 101 * of a #define to make the actual name of the static function different. 102 */ 103 #define cbrt my_cbrt 104 static double cbrt(double x); 105 #endif /* HAVE_CBRT */ 106 107 108 /* 109 * Routines to provide reasonably platform-independent handling of 110 * infinity and NaN. We assume that isinf() and isnan() are available 111 * and work per spec. (On some platforms, we have to supply our own; 112 * see src/port.) However, generating an Infinity or NaN in the first 113 * place is less well standardized; pre-C99 systems tend not to have C99's 114 * INFINITY and NAN macros. We centralize our workarounds for this here. 115 */ 116 117 double 118 get_float8_infinity(void) 119 { 120 #ifdef INFINITY 121 /* C99 standard way */ 122 return (double) INFINITY; 123 #else 124 125 /* 126 * On some platforms, HUGE_VAL is an infinity, elsewhere it's just the 127 * largest normal double. We assume forcing an overflow will get us a 128 * true infinity. 129 */ 130 return (double) (HUGE_VAL * HUGE_VAL); 131 #endif 132 } 133 134 /* 135 * The funny placements of the two #pragmas is necessary because of a 136 * long lived bug in the Microsoft compilers. 137 * See http://support.microsoft.com/kb/120968/en-us for details 138 */ 139 #if (_MSC_VER >= 1800) 140 #pragma warning(disable:4756) 141 #endif 142 float 143 get_float4_infinity(void) 144 { 145 #ifdef INFINITY 146 /* C99 standard way */ 147 return (float) INFINITY; 148 #else 149 #if (_MSC_VER >= 1800) 150 #pragma warning(default:4756) 151 #endif 152 153 /* 154 * On some platforms, HUGE_VAL is an infinity, elsewhere it's just the 155 * largest normal double. We assume forcing an overflow will get us a 156 * true infinity. 157 */ 158 return (float) (HUGE_VAL * HUGE_VAL); 159 #endif 160 } 161 162 double 163 get_float8_nan(void) 164 { 165 /* (double) NAN doesn't work on some NetBSD/MIPS releases */ 166 #if defined(NAN) && !(defined(__NetBSD__) && defined(__mips__)) 167 /* C99 standard way */ 168 return (double) NAN; 169 #else 170 /* Assume we can get a NAN via zero divide */ 171 return (double) (0.0 / 0.0); 172 #endif 173 } 174 175 float 176 get_float4_nan(void) 177 { 178 #ifdef NAN 179 /* C99 standard way */ 180 return (float) NAN; 181 #else 182 /* Assume we can get a NAN via zero divide */ 183 return (float) (0.0 / 0.0); 184 #endif 185 } 186 187 188 /* 189 * Returns -1 if 'val' represents negative infinity, 1 if 'val' 190 * represents (positive) infinity, and 0 otherwise. On some platforms, 191 * this is equivalent to the isinf() macro, but not everywhere: C99 192 * does not specify that isinf() needs to distinguish between positive 193 * and negative infinity. 194 */ 195 int 196 is_infinite(double val) 197 { 198 int inf = isinf(val); 199 200 if (inf == 0) 201 return 0; 202 else if (val > 0) 203 return 1; 204 else 205 return -1; 206 } 207 208 209 /* ========== USER I/O ROUTINES ========== */ 210 211 212 /* 213 * float4in - converts "num" to float4 214 */ 215 Datum 216 float4in(PG_FUNCTION_ARGS) 217 { 218 char *num = PG_GETARG_CSTRING(0); 219 char *orig_num; 220 double val; 221 char *endptr; 222 223 /* 224 * endptr points to the first character _after_ the sequence we recognized 225 * as a valid floating point number. orig_num points to the original input 226 * string. 227 */ 228 orig_num = num; 229 230 /* skip leading whitespace */ 231 while (*num != '\0' && isspace((unsigned char) *num)) 232 num++; 233 234 /* 235 * Check for an empty-string input to begin with, to avoid the vagaries of 236 * strtod() on different platforms. 237 */ 238 if (*num == '\0') 239 ereport(ERROR, 240 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 241 errmsg("invalid input syntax for type %s: \"%s\"", 242 "real", orig_num))); 243 244 errno = 0; 245 val = strtod(num, &endptr); 246 247 /* did we not see anything that looks like a double? */ 248 if (endptr == num || errno != 0) 249 { 250 int save_errno = errno; 251 252 /* 253 * C99 requires that strtod() accept NaN, [+-]Infinity, and [+-]Inf, 254 * but not all platforms support all of these (and some accept them 255 * but set ERANGE anyway...) Therefore, we check for these inputs 256 * ourselves if strtod() fails. 257 * 258 * Note: C99 also requires hexadecimal input as well as some extended 259 * forms of NaN, but we consider these forms unportable and don't try 260 * to support them. You can use 'em if your strtod() takes 'em. 261 */ 262 if (pg_strncasecmp(num, "NaN", 3) == 0) 263 { 264 val = get_float4_nan(); 265 endptr = num + 3; 266 } 267 else if (pg_strncasecmp(num, "Infinity", 8) == 0) 268 { 269 val = get_float4_infinity(); 270 endptr = num + 8; 271 } 272 else if (pg_strncasecmp(num, "+Infinity", 9) == 0) 273 { 274 val = get_float4_infinity(); 275 endptr = num + 9; 276 } 277 else if (pg_strncasecmp(num, "-Infinity", 9) == 0) 278 { 279 val = -get_float4_infinity(); 280 endptr = num + 9; 281 } 282 else if (pg_strncasecmp(num, "inf", 3) == 0) 283 { 284 val = get_float4_infinity(); 285 endptr = num + 3; 286 } 287 else if (pg_strncasecmp(num, "+inf", 4) == 0) 288 { 289 val = get_float4_infinity(); 290 endptr = num + 4; 291 } 292 else if (pg_strncasecmp(num, "-inf", 4) == 0) 293 { 294 val = -get_float4_infinity(); 295 endptr = num + 4; 296 } 297 else if (save_errno == ERANGE) 298 { 299 /* 300 * Some platforms return ERANGE for denormalized numbers (those 301 * that are not zero, but are too close to zero to have full 302 * precision). We'd prefer not to throw error for that, so try to 303 * detect whether it's a "real" out-of-range condition by checking 304 * to see if the result is zero or huge. 305 */ 306 if (val == 0.0 || val >= HUGE_VAL || val <= -HUGE_VAL) 307 ereport(ERROR, 308 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 309 errmsg("\"%s\" is out of range for type real", 310 orig_num))); 311 } 312 else 313 ereport(ERROR, 314 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 315 errmsg("invalid input syntax for type %s: \"%s\"", 316 "real", orig_num))); 317 } 318 #ifdef HAVE_BUGGY_SOLARIS_STRTOD 319 else 320 { 321 /* 322 * Many versions of Solaris have a bug wherein strtod sets endptr to 323 * point one byte beyond the end of the string when given "inf" or 324 * "infinity". 325 */ 326 if (endptr != num && endptr[-1] == '\0') 327 endptr--; 328 } 329 #endif /* HAVE_BUGGY_SOLARIS_STRTOD */ 330 331 /* skip trailing whitespace */ 332 while (*endptr != '\0' && isspace((unsigned char) *endptr)) 333 endptr++; 334 335 /* if there is any junk left at the end of the string, bail out */ 336 if (*endptr != '\0') 337 ereport(ERROR, 338 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 339 errmsg("invalid input syntax for type %s: \"%s\"", 340 "real", orig_num))); 341 342 /* 343 * if we get here, we have a legal double, still need to check to see if 344 * it's a legal float4 345 */ 346 CHECKFLOATVAL((float4) val, isinf(val), val == 0); 347 348 PG_RETURN_FLOAT4((float4) val); 349 } 350 351 /* 352 * float4out - converts a float4 number to a string 353 * using a standard output format 354 */ 355 Datum 356 float4out(PG_FUNCTION_ARGS) 357 { 358 float4 num = PG_GETARG_FLOAT4(0); 359 char *ascii; 360 361 if (isnan(num)) 362 PG_RETURN_CSTRING(pstrdup("NaN")); 363 364 switch (is_infinite(num)) 365 { 366 case 1: 367 ascii = pstrdup("Infinity"); 368 break; 369 case -1: 370 ascii = pstrdup("-Infinity"); 371 break; 372 default: 373 { 374 int ndig = FLT_DIG + extra_float_digits; 375 376 if (ndig < 1) 377 ndig = 1; 378 379 ascii = psprintf("%.*g", ndig, num); 380 } 381 } 382 383 PG_RETURN_CSTRING(ascii); 384 } 385 386 /* 387 * float4recv - converts external binary format to float4 388 */ 389 Datum 390 float4recv(PG_FUNCTION_ARGS) 391 { 392 StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); 393 394 PG_RETURN_FLOAT4(pq_getmsgfloat4(buf)); 395 } 396 397 /* 398 * float4send - converts float4 to binary format 399 */ 400 Datum 401 float4send(PG_FUNCTION_ARGS) 402 { 403 float4 num = PG_GETARG_FLOAT4(0); 404 StringInfoData buf; 405 406 pq_begintypsend(&buf); 407 pq_sendfloat4(&buf, num); 408 PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); 409 } 410 411 /* 412 * float8in - converts "num" to float8 413 */ 414 Datum 415 float8in(PG_FUNCTION_ARGS) 416 { 417 char *num = PG_GETARG_CSTRING(0); 418 419 PG_RETURN_FLOAT8(float8in_internal(num, NULL, "double precision", num)); 420 } 421 422 /* 423 * float8in_internal - guts of float8in() 424 * 425 * This is exposed for use by functions that want a reasonably 426 * platform-independent way of inputting doubles. The behavior is 427 * essentially like strtod + ereport on error, but note the following 428 * differences: 429 * 1. Both leading and trailing whitespace are skipped. 430 * 2. If endptr_p is NULL, we throw error if there's trailing junk. 431 * Otherwise, it's up to the caller to complain about trailing junk. 432 * 3. In event of a syntax error, the report mentions the given type_name 433 * and prints orig_string as the input; this is meant to support use of 434 * this function with types such as "box" and "point", where what we are 435 * parsing here is just a substring of orig_string. 436 * 437 * "num" could validly be declared "const char *", but that results in an 438 * unreasonable amount of extra casting both here and in callers, so we don't. 439 */ 440 double 441 float8in_internal(char *num, char **endptr_p, 442 const char *type_name, const char *orig_string) 443 { 444 double val; 445 char *endptr; 446 447 /* skip leading whitespace */ 448 while (*num != '\0' && isspace((unsigned char) *num)) 449 num++; 450 451 /* 452 * Check for an empty-string input to begin with, to avoid the vagaries of 453 * strtod() on different platforms. 454 */ 455 if (*num == '\0') 456 ereport(ERROR, 457 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 458 errmsg("invalid input syntax for type %s: \"%s\"", 459 type_name, orig_string))); 460 461 errno = 0; 462 val = strtod(num, &endptr); 463 464 /* did we not see anything that looks like a double? */ 465 if (endptr == num || errno != 0) 466 { 467 int save_errno = errno; 468 469 /* 470 * C99 requires that strtod() accept NaN, [+-]Infinity, and [+-]Inf, 471 * but not all platforms support all of these (and some accept them 472 * but set ERANGE anyway...) Therefore, we check for these inputs 473 * ourselves if strtod() fails. 474 * 475 * Note: C99 also requires hexadecimal input as well as some extended 476 * forms of NaN, but we consider these forms unportable and don't try 477 * to support them. You can use 'em if your strtod() takes 'em. 478 */ 479 if (pg_strncasecmp(num, "NaN", 3) == 0) 480 { 481 val = get_float8_nan(); 482 endptr = num + 3; 483 } 484 else if (pg_strncasecmp(num, "Infinity", 8) == 0) 485 { 486 val = get_float8_infinity(); 487 endptr = num + 8; 488 } 489 else if (pg_strncasecmp(num, "+Infinity", 9) == 0) 490 { 491 val = get_float8_infinity(); 492 endptr = num + 9; 493 } 494 else if (pg_strncasecmp(num, "-Infinity", 9) == 0) 495 { 496 val = -get_float8_infinity(); 497 endptr = num + 9; 498 } 499 else if (pg_strncasecmp(num, "inf", 3) == 0) 500 { 501 val = get_float8_infinity(); 502 endptr = num + 3; 503 } 504 else if (pg_strncasecmp(num, "+inf", 4) == 0) 505 { 506 val = get_float8_infinity(); 507 endptr = num + 4; 508 } 509 else if (pg_strncasecmp(num, "-inf", 4) == 0) 510 { 511 val = -get_float8_infinity(); 512 endptr = num + 4; 513 } 514 else if (save_errno == ERANGE) 515 { 516 /* 517 * Some platforms return ERANGE for denormalized numbers (those 518 * that are not zero, but are too close to zero to have full 519 * precision). We'd prefer not to throw error for that, so try to 520 * detect whether it's a "real" out-of-range condition by checking 521 * to see if the result is zero or huge. 522 * 523 * On error, we intentionally complain about double precision not 524 * the given type name, and we print only the part of the string 525 * that is the current number. 526 */ 527 if (val == 0.0 || val >= HUGE_VAL || val <= -HUGE_VAL) 528 { 529 char *errnumber = pstrdup(num); 530 531 errnumber[endptr - num] = '\0'; 532 ereport(ERROR, 533 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 534 errmsg("\"%s\" is out of range for type double precision", 535 errnumber))); 536 } 537 } 538 else 539 ereport(ERROR, 540 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 541 errmsg("invalid input syntax for type %s: \"%s\"", 542 type_name, orig_string))); 543 } 544 #ifdef HAVE_BUGGY_SOLARIS_STRTOD 545 else 546 { 547 /* 548 * Many versions of Solaris have a bug wherein strtod sets endptr to 549 * point one byte beyond the end of the string when given "inf" or 550 * "infinity". 551 */ 552 if (endptr != num && endptr[-1] == '\0') 553 endptr--; 554 } 555 #endif /* HAVE_BUGGY_SOLARIS_STRTOD */ 556 557 /* skip trailing whitespace */ 558 while (*endptr != '\0' && isspace((unsigned char) *endptr)) 559 endptr++; 560 561 /* report stopping point if wanted, else complain if not end of string */ 562 if (endptr_p) 563 *endptr_p = endptr; 564 else if (*endptr != '\0') 565 ereport(ERROR, 566 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 567 errmsg("invalid input syntax for type %s: \"%s\"", 568 type_name, orig_string))); 569 570 return val; 571 } 572 573 /* 574 * float8out - converts float8 number to a string 575 * using a standard output format 576 */ 577 Datum 578 float8out(PG_FUNCTION_ARGS) 579 { 580 float8 num = PG_GETARG_FLOAT8(0); 581 582 PG_RETURN_CSTRING(float8out_internal(num)); 583 } 584 585 /* 586 * float8out_internal - guts of float8out() 587 * 588 * This is exposed for use by functions that want a reasonably 589 * platform-independent way of outputting doubles. 590 * The result is always palloc'd. 591 */ 592 char * 593 float8out_internal(double num) 594 { 595 char *ascii; 596 597 if (isnan(num)) 598 return pstrdup("NaN"); 599 600 switch (is_infinite(num)) 601 { 602 case 1: 603 ascii = pstrdup("Infinity"); 604 break; 605 case -1: 606 ascii = pstrdup("-Infinity"); 607 break; 608 default: 609 { 610 int ndig = DBL_DIG + extra_float_digits; 611 612 if (ndig < 1) 613 ndig = 1; 614 615 ascii = psprintf("%.*g", ndig, num); 616 } 617 } 618 619 return ascii; 620 } 621 622 /* 623 * float8recv - converts external binary format to float8 624 */ 625 Datum 626 float8recv(PG_FUNCTION_ARGS) 627 { 628 StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); 629 630 PG_RETURN_FLOAT8(pq_getmsgfloat8(buf)); 631 } 632 633 /* 634 * float8send - converts float8 to binary format 635 */ 636 Datum 637 float8send(PG_FUNCTION_ARGS) 638 { 639 float8 num = PG_GETARG_FLOAT8(0); 640 StringInfoData buf; 641 642 pq_begintypsend(&buf); 643 pq_sendfloat8(&buf, num); 644 PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); 645 } 646 647 648 /* ========== PUBLIC ROUTINES ========== */ 649 650 651 /* 652 * ====================== 653 * FLOAT4 BASE OPERATIONS 654 * ====================== 655 */ 656 657 /* 658 * float4abs - returns |arg1| (absolute value) 659 */ 660 Datum 661 float4abs(PG_FUNCTION_ARGS) 662 { 663 float4 arg1 = PG_GETARG_FLOAT4(0); 664 665 PG_RETURN_FLOAT4((float4) fabs(arg1)); 666 } 667 668 /* 669 * float4um - returns -arg1 (unary minus) 670 */ 671 Datum 672 float4um(PG_FUNCTION_ARGS) 673 { 674 float4 arg1 = PG_GETARG_FLOAT4(0); 675 float4 result; 676 677 result = -arg1; 678 PG_RETURN_FLOAT4(result); 679 } 680 681 Datum 682 float4up(PG_FUNCTION_ARGS) 683 { 684 float4 arg = PG_GETARG_FLOAT4(0); 685 686 PG_RETURN_FLOAT4(arg); 687 } 688 689 Datum 690 float4larger(PG_FUNCTION_ARGS) 691 { 692 float4 arg1 = PG_GETARG_FLOAT4(0); 693 float4 arg2 = PG_GETARG_FLOAT4(1); 694 float4 result; 695 696 if (float4_cmp_internal(arg1, arg2) > 0) 697 result = arg1; 698 else 699 result = arg2; 700 PG_RETURN_FLOAT4(result); 701 } 702 703 Datum 704 float4smaller(PG_FUNCTION_ARGS) 705 { 706 float4 arg1 = PG_GETARG_FLOAT4(0); 707 float4 arg2 = PG_GETARG_FLOAT4(1); 708 float4 result; 709 710 if (float4_cmp_internal(arg1, arg2) < 0) 711 result = arg1; 712 else 713 result = arg2; 714 PG_RETURN_FLOAT4(result); 715 } 716 717 /* 718 * ====================== 719 * FLOAT8 BASE OPERATIONS 720 * ====================== 721 */ 722 723 /* 724 * float8abs - returns |arg1| (absolute value) 725 */ 726 Datum 727 float8abs(PG_FUNCTION_ARGS) 728 { 729 float8 arg1 = PG_GETARG_FLOAT8(0); 730 731 PG_RETURN_FLOAT8(fabs(arg1)); 732 } 733 734 735 /* 736 * float8um - returns -arg1 (unary minus) 737 */ 738 Datum 739 float8um(PG_FUNCTION_ARGS) 740 { 741 float8 arg1 = PG_GETARG_FLOAT8(0); 742 float8 result; 743 744 result = -arg1; 745 PG_RETURN_FLOAT8(result); 746 } 747 748 Datum 749 float8up(PG_FUNCTION_ARGS) 750 { 751 float8 arg = PG_GETARG_FLOAT8(0); 752 753 PG_RETURN_FLOAT8(arg); 754 } 755 756 Datum 757 float8larger(PG_FUNCTION_ARGS) 758 { 759 float8 arg1 = PG_GETARG_FLOAT8(0); 760 float8 arg2 = PG_GETARG_FLOAT8(1); 761 float8 result; 762 763 if (float8_cmp_internal(arg1, arg2) > 0) 764 result = arg1; 765 else 766 result = arg2; 767 PG_RETURN_FLOAT8(result); 768 } 769 770 Datum 771 float8smaller(PG_FUNCTION_ARGS) 772 { 773 float8 arg1 = PG_GETARG_FLOAT8(0); 774 float8 arg2 = PG_GETARG_FLOAT8(1); 775 float8 result; 776 777 if (float8_cmp_internal(arg1, arg2) < 0) 778 result = arg1; 779 else 780 result = arg2; 781 PG_RETURN_FLOAT8(result); 782 } 783 784 785 /* 786 * ==================== 787 * ARITHMETIC OPERATORS 788 * ==================== 789 */ 790 791 /* 792 * float4pl - returns arg1 + arg2 793 * float4mi - returns arg1 - arg2 794 * float4mul - returns arg1 * arg2 795 * float4div - returns arg1 / arg2 796 */ 797 Datum 798 float4pl(PG_FUNCTION_ARGS) 799 { 800 float4 arg1 = PG_GETARG_FLOAT4(0); 801 float4 arg2 = PG_GETARG_FLOAT4(1); 802 float4 result; 803 804 result = arg1 + arg2; 805 806 /* 807 * There isn't any way to check for underflow of addition/subtraction 808 * because numbers near the underflow value have already been rounded to 809 * the point where we can't detect that the two values were originally 810 * different, e.g. on x86, '1e-45'::float4 == '2e-45'::float4 == 811 * 1.4013e-45. 812 */ 813 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true); 814 PG_RETURN_FLOAT4(result); 815 } 816 817 Datum 818 float4mi(PG_FUNCTION_ARGS) 819 { 820 float4 arg1 = PG_GETARG_FLOAT4(0); 821 float4 arg2 = PG_GETARG_FLOAT4(1); 822 float4 result; 823 824 result = arg1 - arg2; 825 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true); 826 PG_RETURN_FLOAT4(result); 827 } 828 829 Datum 830 float4mul(PG_FUNCTION_ARGS) 831 { 832 float4 arg1 = PG_GETARG_FLOAT4(0); 833 float4 arg2 = PG_GETARG_FLOAT4(1); 834 float4 result; 835 836 result = arg1 * arg2; 837 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), 838 arg1 == 0 || arg2 == 0); 839 PG_RETURN_FLOAT4(result); 840 } 841 842 Datum 843 float4div(PG_FUNCTION_ARGS) 844 { 845 float4 arg1 = PG_GETARG_FLOAT4(0); 846 float4 arg2 = PG_GETARG_FLOAT4(1); 847 float4 result; 848 849 if (arg2 == 0.0) 850 ereport(ERROR, 851 (errcode(ERRCODE_DIVISION_BY_ZERO), 852 errmsg("division by zero"))); 853 854 result = arg1 / arg2; 855 856 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0); 857 PG_RETURN_FLOAT4(result); 858 } 859 860 /* 861 * float8pl - returns arg1 + arg2 862 * float8mi - returns arg1 - arg2 863 * float8mul - returns arg1 * arg2 864 * float8div - returns arg1 / arg2 865 */ 866 Datum 867 float8pl(PG_FUNCTION_ARGS) 868 { 869 float8 arg1 = PG_GETARG_FLOAT8(0); 870 float8 arg2 = PG_GETARG_FLOAT8(1); 871 float8 result; 872 873 result = arg1 + arg2; 874 875 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true); 876 PG_RETURN_FLOAT8(result); 877 } 878 879 Datum 880 float8mi(PG_FUNCTION_ARGS) 881 { 882 float8 arg1 = PG_GETARG_FLOAT8(0); 883 float8 arg2 = PG_GETARG_FLOAT8(1); 884 float8 result; 885 886 result = arg1 - arg2; 887 888 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true); 889 PG_RETURN_FLOAT8(result); 890 } 891 892 Datum 893 float8mul(PG_FUNCTION_ARGS) 894 { 895 float8 arg1 = PG_GETARG_FLOAT8(0); 896 float8 arg2 = PG_GETARG_FLOAT8(1); 897 float8 result; 898 899 result = arg1 * arg2; 900 901 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), 902 arg1 == 0 || arg2 == 0); 903 PG_RETURN_FLOAT8(result); 904 } 905 906 Datum 907 float8div(PG_FUNCTION_ARGS) 908 { 909 float8 arg1 = PG_GETARG_FLOAT8(0); 910 float8 arg2 = PG_GETARG_FLOAT8(1); 911 float8 result; 912 913 if (arg2 == 0.0) 914 ereport(ERROR, 915 (errcode(ERRCODE_DIVISION_BY_ZERO), 916 errmsg("division by zero"))); 917 918 result = arg1 / arg2; 919 920 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0); 921 PG_RETURN_FLOAT8(result); 922 } 923 924 925 /* 926 * ==================== 927 * COMPARISON OPERATORS 928 * ==================== 929 */ 930 931 /* 932 * float4{eq,ne,lt,le,gt,ge} - float4/float4 comparison operations 933 */ 934 int 935 float4_cmp_internal(float4 a, float4 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 float4eq(PG_FUNCTION_ARGS) 966 { 967 float4 arg1 = PG_GETARG_FLOAT4(0); 968 float4 arg2 = PG_GETARG_FLOAT4(1); 969 970 PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) == 0); 971 } 972 973 Datum 974 float4ne(PG_FUNCTION_ARGS) 975 { 976 float4 arg1 = PG_GETARG_FLOAT4(0); 977 float4 arg2 = PG_GETARG_FLOAT4(1); 978 979 PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) != 0); 980 } 981 982 Datum 983 float4lt(PG_FUNCTION_ARGS) 984 { 985 float4 arg1 = PG_GETARG_FLOAT4(0); 986 float4 arg2 = PG_GETARG_FLOAT4(1); 987 988 PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) < 0); 989 } 990 991 Datum 992 float4le(PG_FUNCTION_ARGS) 993 { 994 float4 arg1 = PG_GETARG_FLOAT4(0); 995 float4 arg2 = PG_GETARG_FLOAT4(1); 996 997 PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) <= 0); 998 } 999 1000 Datum 1001 float4gt(PG_FUNCTION_ARGS) 1002 { 1003 float4 arg1 = PG_GETARG_FLOAT4(0); 1004 float4 arg2 = PG_GETARG_FLOAT4(1); 1005 1006 PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) > 0); 1007 } 1008 1009 Datum 1010 float4ge(PG_FUNCTION_ARGS) 1011 { 1012 float4 arg1 = PG_GETARG_FLOAT4(0); 1013 float4 arg2 = PG_GETARG_FLOAT4(1); 1014 1015 PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) >= 0); 1016 } 1017 1018 Datum 1019 btfloat4cmp(PG_FUNCTION_ARGS) 1020 { 1021 float4 arg1 = PG_GETARG_FLOAT4(0); 1022 float4 arg2 = PG_GETARG_FLOAT4(1); 1023 1024 PG_RETURN_INT32(float4_cmp_internal(arg1, arg2)); 1025 } 1026 1027 static int 1028 btfloat4fastcmp(Datum x, Datum y, SortSupport ssup) 1029 { 1030 float4 arg1 = DatumGetFloat4(x); 1031 float4 arg2 = DatumGetFloat4(y); 1032 1033 return float4_cmp_internal(arg1, arg2); 1034 } 1035 1036 Datum 1037 btfloat4sortsupport(PG_FUNCTION_ARGS) 1038 { 1039 SortSupport ssup = (SortSupport) PG_GETARG_POINTER(0); 1040 1041 ssup->comparator = btfloat4fastcmp; 1042 PG_RETURN_VOID(); 1043 } 1044 1045 /* 1046 * float8{eq,ne,lt,le,gt,ge} - float8/float8 comparison operations 1047 */ 1048 int 1049 float8_cmp_internal(float8 a, float8 b) 1050 { 1051 /* 1052 * We consider all NANs to be equal and larger than any non-NAN. This is 1053 * somewhat arbitrary; the important thing is to have a consistent sort 1054 * order. 1055 */ 1056 if (isnan(a)) 1057 { 1058 if (isnan(b)) 1059 return 0; /* NAN = NAN */ 1060 else 1061 return 1; /* NAN > non-NAN */ 1062 } 1063 else if (isnan(b)) 1064 { 1065 return -1; /* non-NAN < NAN */ 1066 } 1067 else 1068 { 1069 if (a > b) 1070 return 1; 1071 else if (a < b) 1072 return -1; 1073 else 1074 return 0; 1075 } 1076 } 1077 1078 Datum 1079 float8eq(PG_FUNCTION_ARGS) 1080 { 1081 float8 arg1 = PG_GETARG_FLOAT8(0); 1082 float8 arg2 = PG_GETARG_FLOAT8(1); 1083 1084 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0); 1085 } 1086 1087 Datum 1088 float8ne(PG_FUNCTION_ARGS) 1089 { 1090 float8 arg1 = PG_GETARG_FLOAT8(0); 1091 float8 arg2 = PG_GETARG_FLOAT8(1); 1092 1093 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0); 1094 } 1095 1096 Datum 1097 float8lt(PG_FUNCTION_ARGS) 1098 { 1099 float8 arg1 = PG_GETARG_FLOAT8(0); 1100 float8 arg2 = PG_GETARG_FLOAT8(1); 1101 1102 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0); 1103 } 1104 1105 Datum 1106 float8le(PG_FUNCTION_ARGS) 1107 { 1108 float8 arg1 = PG_GETARG_FLOAT8(0); 1109 float8 arg2 = PG_GETARG_FLOAT8(1); 1110 1111 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0); 1112 } 1113 1114 Datum 1115 float8gt(PG_FUNCTION_ARGS) 1116 { 1117 float8 arg1 = PG_GETARG_FLOAT8(0); 1118 float8 arg2 = PG_GETARG_FLOAT8(1); 1119 1120 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0); 1121 } 1122 1123 Datum 1124 float8ge(PG_FUNCTION_ARGS) 1125 { 1126 float8 arg1 = PG_GETARG_FLOAT8(0); 1127 float8 arg2 = PG_GETARG_FLOAT8(1); 1128 1129 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0); 1130 } 1131 1132 Datum 1133 btfloat8cmp(PG_FUNCTION_ARGS) 1134 { 1135 float8 arg1 = PG_GETARG_FLOAT8(0); 1136 float8 arg2 = PG_GETARG_FLOAT8(1); 1137 1138 PG_RETURN_INT32(float8_cmp_internal(arg1, arg2)); 1139 } 1140 1141 static int 1142 btfloat8fastcmp(Datum x, Datum y, SortSupport ssup) 1143 { 1144 float8 arg1 = DatumGetFloat8(x); 1145 float8 arg2 = DatumGetFloat8(y); 1146 1147 return float8_cmp_internal(arg1, arg2); 1148 } 1149 1150 Datum 1151 btfloat8sortsupport(PG_FUNCTION_ARGS) 1152 { 1153 SortSupport ssup = (SortSupport) PG_GETARG_POINTER(0); 1154 1155 ssup->comparator = btfloat8fastcmp; 1156 PG_RETURN_VOID(); 1157 } 1158 1159 Datum 1160 btfloat48cmp(PG_FUNCTION_ARGS) 1161 { 1162 float4 arg1 = PG_GETARG_FLOAT4(0); 1163 float8 arg2 = PG_GETARG_FLOAT8(1); 1164 1165 /* widen float4 to float8 and then compare */ 1166 PG_RETURN_INT32(float8_cmp_internal(arg1, arg2)); 1167 } 1168 1169 Datum 1170 btfloat84cmp(PG_FUNCTION_ARGS) 1171 { 1172 float8 arg1 = PG_GETARG_FLOAT8(0); 1173 float4 arg2 = PG_GETARG_FLOAT4(1); 1174 1175 /* widen float4 to float8 and then compare */ 1176 PG_RETURN_INT32(float8_cmp_internal(arg1, arg2)); 1177 } 1178 1179 /* 1180 * in_range support function for float8. 1181 * 1182 * Note: we needn't supply a float8_float4 variant, as implicit coercion 1183 * of the offset value takes care of that scenario just as well. 1184 */ 1185 Datum 1186 in_range_float8_float8(PG_FUNCTION_ARGS) 1187 { 1188 float8 val = PG_GETARG_FLOAT8(0); 1189 float8 base = PG_GETARG_FLOAT8(1); 1190 float8 offset = PG_GETARG_FLOAT8(2); 1191 bool sub = PG_GETARG_BOOL(3); 1192 bool less = PG_GETARG_BOOL(4); 1193 float8 sum; 1194 1195 /* 1196 * Reject negative or NaN offset. Negative is per spec, and NaN is 1197 * because appropriate semantics for that seem non-obvious. 1198 */ 1199 if (isnan(offset) || offset < 0) 1200 ereport(ERROR, 1201 (errcode(ERRCODE_INVALID_PRECEDING_OR_FOLLOWING_SIZE), 1202 errmsg("invalid preceding or following size in window function"))); 1203 1204 /* 1205 * Deal with cases where val and/or base is NaN, following the rule that 1206 * NaN sorts after non-NaN (cf float8_cmp_internal). The offset cannot 1207 * affect the conclusion. 1208 */ 1209 if (isnan(val)) 1210 { 1211 if (isnan(base)) 1212 PG_RETURN_BOOL(true); /* NAN = NAN */ 1213 else 1214 PG_RETURN_BOOL(!less); /* NAN > non-NAN */ 1215 } 1216 else if (isnan(base)) 1217 { 1218 PG_RETURN_BOOL(less); /* non-NAN < NAN */ 1219 } 1220 1221 /* 1222 * Deal with infinite offset (necessarily +inf, at this point). We must 1223 * special-case this because if base happens to be -inf, their sum would 1224 * be NaN, which is an overflow-ish condition we should avoid. 1225 */ 1226 if (isinf(offset)) 1227 { 1228 PG_RETURN_BOOL(sub ? !less : less); 1229 } 1230 1231 /* 1232 * Otherwise it should be safe to compute base +/- offset. We trust the 1233 * FPU to cope if base is +/-inf or the true sum would overflow, and 1234 * produce a suitably signed infinity, which will compare properly against 1235 * val whether or not that's infinity. 1236 */ 1237 if (sub) 1238 sum = base - offset; 1239 else 1240 sum = base + offset; 1241 1242 if (less) 1243 PG_RETURN_BOOL(val <= sum); 1244 else 1245 PG_RETURN_BOOL(val >= sum); 1246 } 1247 1248 /* 1249 * in_range support function for float4. 1250 * 1251 * We would need a float4_float8 variant in any case, so we supply that and 1252 * let implicit coercion take care of the float4_float4 case. 1253 */ 1254 Datum 1255 in_range_float4_float8(PG_FUNCTION_ARGS) 1256 { 1257 float4 val = PG_GETARG_FLOAT4(0); 1258 float4 base = PG_GETARG_FLOAT4(1); 1259 float8 offset = PG_GETARG_FLOAT8(2); 1260 bool sub = PG_GETARG_BOOL(3); 1261 bool less = PG_GETARG_BOOL(4); 1262 float8 sum; 1263 1264 /* 1265 * Reject negative or NaN offset. Negative is per spec, and NaN is 1266 * because appropriate semantics for that seem non-obvious. 1267 */ 1268 if (isnan(offset) || offset < 0) 1269 ereport(ERROR, 1270 (errcode(ERRCODE_INVALID_PRECEDING_OR_FOLLOWING_SIZE), 1271 errmsg("invalid preceding or following size in window function"))); 1272 1273 /* 1274 * Deal with cases where val and/or base is NaN, following the rule that 1275 * NaN sorts after non-NaN (cf float8_cmp_internal). The offset cannot 1276 * affect the conclusion. 1277 */ 1278 if (isnan(val)) 1279 { 1280 if (isnan(base)) 1281 PG_RETURN_BOOL(true); /* NAN = NAN */ 1282 else 1283 PG_RETURN_BOOL(!less); /* NAN > non-NAN */ 1284 } 1285 else if (isnan(base)) 1286 { 1287 PG_RETURN_BOOL(less); /* non-NAN < NAN */ 1288 } 1289 1290 /* 1291 * Deal with infinite offset (necessarily +inf, at this point). We must 1292 * special-case this because if base happens to be -inf, their sum would 1293 * be NaN, which is an overflow-ish condition we should avoid. 1294 */ 1295 if (isinf(offset)) 1296 { 1297 PG_RETURN_BOOL(sub ? !less : less); 1298 } 1299 1300 /* 1301 * Otherwise it should be safe to compute base +/- offset. We trust the 1302 * FPU to cope if base is +/-inf or the true sum would overflow, and 1303 * produce a suitably signed infinity, which will compare properly against 1304 * val whether or not that's infinity. 1305 */ 1306 if (sub) 1307 sum = base - offset; 1308 else 1309 sum = base + offset; 1310 1311 if (less) 1312 PG_RETURN_BOOL(val <= sum); 1313 else 1314 PG_RETURN_BOOL(val >= sum); 1315 } 1316 1317 1318 /* 1319 * =================== 1320 * CONVERSION ROUTINES 1321 * =================== 1322 */ 1323 1324 /* 1325 * ftod - converts a float4 number to a float8 number 1326 */ 1327 Datum 1328 ftod(PG_FUNCTION_ARGS) 1329 { 1330 float4 num = PG_GETARG_FLOAT4(0); 1331 1332 PG_RETURN_FLOAT8((float8) num); 1333 } 1334 1335 1336 /* 1337 * dtof - converts a float8 number to a float4 number 1338 */ 1339 Datum 1340 dtof(PG_FUNCTION_ARGS) 1341 { 1342 float8 num = PG_GETARG_FLOAT8(0); 1343 1344 CHECKFLOATVAL((float4) num, isinf(num), num == 0); 1345 1346 PG_RETURN_FLOAT4((float4) num); 1347 } 1348 1349 1350 /* 1351 * dtoi4 - converts a float8 number to an int4 number 1352 */ 1353 Datum 1354 dtoi4(PG_FUNCTION_ARGS) 1355 { 1356 float8 num = PG_GETARG_FLOAT8(0); 1357 1358 /* 1359 * Get rid of any fractional part in the input. This is so we don't fail 1360 * on just-out-of-range values that would round into range. Note 1361 * assumption that rint() will pass through a NaN or Inf unchanged. 1362 */ 1363 num = rint(num); 1364 1365 /* Range check */ 1366 if (unlikely(isnan(num) || !FLOAT8_FITS_IN_INT32(num))) 1367 ereport(ERROR, 1368 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 1369 errmsg("integer out of range"))); 1370 1371 PG_RETURN_INT32((int32) num); 1372 } 1373 1374 1375 /* 1376 * dtoi2 - converts a float8 number to an int2 number 1377 */ 1378 Datum 1379 dtoi2(PG_FUNCTION_ARGS) 1380 { 1381 float8 num = PG_GETARG_FLOAT8(0); 1382 1383 /* 1384 * Get rid of any fractional part in the input. This is so we don't fail 1385 * on just-out-of-range values that would round into range. Note 1386 * assumption that rint() will pass through a NaN or Inf unchanged. 1387 */ 1388 num = rint(num); 1389 1390 /* Range check */ 1391 if (unlikely(isnan(num) || !FLOAT8_FITS_IN_INT16(num))) 1392 ereport(ERROR, 1393 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 1394 errmsg("smallint out of range"))); 1395 1396 PG_RETURN_INT16((int16) num); 1397 } 1398 1399 1400 /* 1401 * i4tod - converts an int4 number to a float8 number 1402 */ 1403 Datum 1404 i4tod(PG_FUNCTION_ARGS) 1405 { 1406 int32 num = PG_GETARG_INT32(0); 1407 1408 PG_RETURN_FLOAT8((float8) num); 1409 } 1410 1411 1412 /* 1413 * i2tod - converts an int2 number to a float8 number 1414 */ 1415 Datum 1416 i2tod(PG_FUNCTION_ARGS) 1417 { 1418 int16 num = PG_GETARG_INT16(0); 1419 1420 PG_RETURN_FLOAT8((float8) num); 1421 } 1422 1423 1424 /* 1425 * ftoi4 - converts a float4 number to an int4 number 1426 */ 1427 Datum 1428 ftoi4(PG_FUNCTION_ARGS) 1429 { 1430 float4 num = PG_GETARG_FLOAT4(0); 1431 1432 /* 1433 * Get rid of any fractional part in the input. This is so we don't fail 1434 * on just-out-of-range values that would round into range. Note 1435 * assumption that rint() will pass through a NaN or Inf unchanged. 1436 */ 1437 num = rint(num); 1438 1439 /* Range check */ 1440 if (unlikely(isnan(num) || !FLOAT4_FITS_IN_INT32(num))) 1441 ereport(ERROR, 1442 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 1443 errmsg("integer out of range"))); 1444 1445 PG_RETURN_INT32((int32) num); 1446 } 1447 1448 1449 /* 1450 * ftoi2 - converts a float4 number to an int2 number 1451 */ 1452 Datum 1453 ftoi2(PG_FUNCTION_ARGS) 1454 { 1455 float4 num = PG_GETARG_FLOAT4(0); 1456 1457 /* 1458 * Get rid of any fractional part in the input. This is so we don't fail 1459 * on just-out-of-range values that would round into range. Note 1460 * assumption that rint() will pass through a NaN or Inf unchanged. 1461 */ 1462 num = rint(num); 1463 1464 /* Range check */ 1465 if (unlikely(isnan(num) || !FLOAT4_FITS_IN_INT16(num))) 1466 ereport(ERROR, 1467 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 1468 errmsg("smallint out of range"))); 1469 1470 PG_RETURN_INT16((int16) num); 1471 } 1472 1473 1474 /* 1475 * i4tof - converts an int4 number to a float4 number 1476 */ 1477 Datum 1478 i4tof(PG_FUNCTION_ARGS) 1479 { 1480 int32 num = PG_GETARG_INT32(0); 1481 1482 PG_RETURN_FLOAT4((float4) num); 1483 } 1484 1485 1486 /* 1487 * i2tof - converts an int2 number to a float4 number 1488 */ 1489 Datum 1490 i2tof(PG_FUNCTION_ARGS) 1491 { 1492 int16 num = PG_GETARG_INT16(0); 1493 1494 PG_RETURN_FLOAT4((float4) num); 1495 } 1496 1497 1498 /* 1499 * ======================= 1500 * RANDOM FLOAT8 OPERATORS 1501 * ======================= 1502 */ 1503 1504 /* 1505 * dround - returns ROUND(arg1) 1506 */ 1507 Datum 1508 dround(PG_FUNCTION_ARGS) 1509 { 1510 float8 arg1 = PG_GETARG_FLOAT8(0); 1511 1512 PG_RETURN_FLOAT8(rint(arg1)); 1513 } 1514 1515 /* 1516 * dceil - returns the smallest integer greater than or 1517 * equal to the specified float 1518 */ 1519 Datum 1520 dceil(PG_FUNCTION_ARGS) 1521 { 1522 float8 arg1 = PG_GETARG_FLOAT8(0); 1523 1524 PG_RETURN_FLOAT8(ceil(arg1)); 1525 } 1526 1527 /* 1528 * dfloor - returns the largest integer lesser than or 1529 * equal to the specified float 1530 */ 1531 Datum 1532 dfloor(PG_FUNCTION_ARGS) 1533 { 1534 float8 arg1 = PG_GETARG_FLOAT8(0); 1535 1536 PG_RETURN_FLOAT8(floor(arg1)); 1537 } 1538 1539 /* 1540 * dsign - returns -1 if the argument is less than 0, 0 1541 * if the argument is equal to 0, and 1 if the 1542 * argument is greater than zero. 1543 */ 1544 Datum 1545 dsign(PG_FUNCTION_ARGS) 1546 { 1547 float8 arg1 = PG_GETARG_FLOAT8(0); 1548 float8 result; 1549 1550 if (arg1 > 0) 1551 result = 1.0; 1552 else if (arg1 < 0) 1553 result = -1.0; 1554 else 1555 result = 0.0; 1556 1557 PG_RETURN_FLOAT8(result); 1558 } 1559 1560 /* 1561 * dtrunc - returns truncation-towards-zero of arg1, 1562 * arg1 >= 0 ... the greatest integer less 1563 * than or equal to arg1 1564 * arg1 < 0 ... the least integer greater 1565 * than or equal to arg1 1566 */ 1567 Datum 1568 dtrunc(PG_FUNCTION_ARGS) 1569 { 1570 float8 arg1 = PG_GETARG_FLOAT8(0); 1571 float8 result; 1572 1573 if (arg1 >= 0) 1574 result = floor(arg1); 1575 else 1576 result = -floor(-arg1); 1577 1578 PG_RETURN_FLOAT8(result); 1579 } 1580 1581 1582 /* 1583 * dsqrt - returns square root of arg1 1584 */ 1585 Datum 1586 dsqrt(PG_FUNCTION_ARGS) 1587 { 1588 float8 arg1 = PG_GETARG_FLOAT8(0); 1589 float8 result; 1590 1591 if (arg1 < 0) 1592 ereport(ERROR, 1593 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION), 1594 errmsg("cannot take square root of a negative number"))); 1595 1596 result = sqrt(arg1); 1597 1598 CHECKFLOATVAL(result, isinf(arg1), arg1 == 0); 1599 PG_RETURN_FLOAT8(result); 1600 } 1601 1602 1603 /* 1604 * dcbrt - returns cube root of arg1 1605 */ 1606 Datum 1607 dcbrt(PG_FUNCTION_ARGS) 1608 { 1609 float8 arg1 = PG_GETARG_FLOAT8(0); 1610 float8 result; 1611 1612 result = cbrt(arg1); 1613 CHECKFLOATVAL(result, isinf(arg1), arg1 == 0); 1614 PG_RETURN_FLOAT8(result); 1615 } 1616 1617 1618 /* 1619 * dpow - returns pow(arg1,arg2) 1620 */ 1621 Datum 1622 dpow(PG_FUNCTION_ARGS) 1623 { 1624 float8 arg1 = PG_GETARG_FLOAT8(0); 1625 float8 arg2 = PG_GETARG_FLOAT8(1); 1626 float8 result; 1627 1628 /* 1629 * The POSIX spec says that NaN ^ 0 = 1, and 1 ^ NaN = 1, while all other 1630 * cases with NaN inputs yield NaN (with no error). Many older platforms 1631 * get one or more of these cases wrong, so deal with them via explicit 1632 * logic rather than trusting pow(3). 1633 */ 1634 if (isnan(arg1)) 1635 { 1636 if (isnan(arg2) || arg2 != 0.0) 1637 PG_RETURN_FLOAT8(get_float8_nan()); 1638 PG_RETURN_FLOAT8(1.0); 1639 } 1640 if (isnan(arg2)) 1641 { 1642 if (arg1 != 1.0) 1643 PG_RETURN_FLOAT8(get_float8_nan()); 1644 PG_RETURN_FLOAT8(1.0); 1645 } 1646 1647 /* 1648 * The SQL spec requires that we emit a particular SQLSTATE error code for 1649 * certain error conditions. Specifically, we don't return a 1650 * divide-by-zero error code for 0 ^ -1. 1651 */ 1652 if (arg1 == 0 && arg2 < 0) 1653 ereport(ERROR, 1654 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION), 1655 errmsg("zero raised to a negative power is undefined"))); 1656 if (arg1 < 0 && floor(arg2) != arg2) 1657 ereport(ERROR, 1658 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION), 1659 errmsg("a negative number raised to a non-integer power yields a complex result"))); 1660 1661 /* 1662 * pow() sets errno only on some platforms, depending on whether it 1663 * follows _IEEE_, _POSIX_, _XOPEN_, or _SVID_, so we try to avoid using 1664 * errno. However, some platform/CPU combinations return errno == EDOM 1665 * and result == NaN for negative arg1 and very large arg2 (they must be 1666 * using something different from our floor() test to decide it's 1667 * invalid). Other platforms (HPPA) return errno == ERANGE and a large 1668 * (HUGE_VAL) but finite result to signal overflow. 1669 */ 1670 errno = 0; 1671 result = pow(arg1, arg2); 1672 if (errno == EDOM && isnan(result)) 1673 { 1674 if ((fabs(arg1) > 1 && arg2 >= 0) || (fabs(arg1) < 1 && arg2 < 0)) 1675 /* The sign of Inf is not significant in this case. */ 1676 result = get_float8_infinity(); 1677 else if (fabs(arg1) != 1) 1678 result = 0; 1679 else 1680 result = 1; 1681 } 1682 else if (errno == ERANGE && result != 0 && !isinf(result)) 1683 result = get_float8_infinity(); 1684 1685 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0); 1686 PG_RETURN_FLOAT8(result); 1687 } 1688 1689 1690 /* 1691 * dexp - returns the exponential function of arg1 1692 */ 1693 Datum 1694 dexp(PG_FUNCTION_ARGS) 1695 { 1696 float8 arg1 = PG_GETARG_FLOAT8(0); 1697 float8 result; 1698 1699 errno = 0; 1700 result = exp(arg1); 1701 if (errno == ERANGE && result != 0 && !isinf(result)) 1702 result = get_float8_infinity(); 1703 1704 CHECKFLOATVAL(result, isinf(arg1), false); 1705 PG_RETURN_FLOAT8(result); 1706 } 1707 1708 1709 /* 1710 * dlog1 - returns the natural logarithm of arg1 1711 */ 1712 Datum 1713 dlog1(PG_FUNCTION_ARGS) 1714 { 1715 float8 arg1 = PG_GETARG_FLOAT8(0); 1716 float8 result; 1717 1718 /* 1719 * Emit particular SQLSTATE error codes for ln(). This is required by the 1720 * SQL standard. 1721 */ 1722 if (arg1 == 0.0) 1723 ereport(ERROR, 1724 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG), 1725 errmsg("cannot take logarithm of zero"))); 1726 if (arg1 < 0) 1727 ereport(ERROR, 1728 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG), 1729 errmsg("cannot take logarithm of a negative number"))); 1730 1731 result = log(arg1); 1732 1733 CHECKFLOATVAL(result, isinf(arg1), arg1 == 1); 1734 PG_RETURN_FLOAT8(result); 1735 } 1736 1737 1738 /* 1739 * dlog10 - returns the base 10 logarithm of arg1 1740 */ 1741 Datum 1742 dlog10(PG_FUNCTION_ARGS) 1743 { 1744 float8 arg1 = PG_GETARG_FLOAT8(0); 1745 float8 result; 1746 1747 /* 1748 * Emit particular SQLSTATE error codes for log(). The SQL spec doesn't 1749 * define log(), but it does define ln(), so it makes sense to emit the 1750 * same error code for an analogous error condition. 1751 */ 1752 if (arg1 == 0.0) 1753 ereport(ERROR, 1754 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG), 1755 errmsg("cannot take logarithm of zero"))); 1756 if (arg1 < 0) 1757 ereport(ERROR, 1758 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG), 1759 errmsg("cannot take logarithm of a negative number"))); 1760 1761 result = log10(arg1); 1762 1763 CHECKFLOATVAL(result, isinf(arg1), arg1 == 1); 1764 PG_RETURN_FLOAT8(result); 1765 } 1766 1767 1768 /* 1769 * dacos - returns the arccos of arg1 (radians) 1770 */ 1771 Datum 1772 dacos(PG_FUNCTION_ARGS) 1773 { 1774 float8 arg1 = PG_GETARG_FLOAT8(0); 1775 float8 result; 1776 1777 /* Per the POSIX spec, return NaN if the input is NaN */ 1778 if (isnan(arg1)) 1779 PG_RETURN_FLOAT8(get_float8_nan()); 1780 1781 /* 1782 * The principal branch of the inverse cosine function maps values in the 1783 * range [-1, 1] to values in the range [0, Pi], so we should reject any 1784 * inputs outside that range and the result will always be finite. 1785 */ 1786 if (arg1 < -1.0 || arg1 > 1.0) 1787 ereport(ERROR, 1788 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 1789 errmsg("input is out of range"))); 1790 1791 result = acos(arg1); 1792 1793 CHECKFLOATVAL(result, false, true); 1794 PG_RETURN_FLOAT8(result); 1795 } 1796 1797 1798 /* 1799 * dasin - returns the arcsin of arg1 (radians) 1800 */ 1801 Datum 1802 dasin(PG_FUNCTION_ARGS) 1803 { 1804 float8 arg1 = PG_GETARG_FLOAT8(0); 1805 float8 result; 1806 1807 /* Per the POSIX spec, return NaN if the input is NaN */ 1808 if (isnan(arg1)) 1809 PG_RETURN_FLOAT8(get_float8_nan()); 1810 1811 /* 1812 * The principal branch of the inverse sine function maps values in the 1813 * range [-1, 1] to values in the range [-Pi/2, Pi/2], so we should reject 1814 * any inputs outside that range and the result will always be finite. 1815 */ 1816 if (arg1 < -1.0 || arg1 > 1.0) 1817 ereport(ERROR, 1818 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 1819 errmsg("input is out of range"))); 1820 1821 result = asin(arg1); 1822 1823 CHECKFLOATVAL(result, false, true); 1824 PG_RETURN_FLOAT8(result); 1825 } 1826 1827 1828 /* 1829 * datan - returns the arctan of arg1 (radians) 1830 */ 1831 Datum 1832 datan(PG_FUNCTION_ARGS) 1833 { 1834 float8 arg1 = PG_GETARG_FLOAT8(0); 1835 float8 result; 1836 1837 /* Per the POSIX spec, return NaN if the input is NaN */ 1838 if (isnan(arg1)) 1839 PG_RETURN_FLOAT8(get_float8_nan()); 1840 1841 /* 1842 * The principal branch of the inverse tangent function maps all inputs to 1843 * values in the range [-Pi/2, Pi/2], so the result should always be 1844 * finite, even if the input is infinite. 1845 */ 1846 result = atan(arg1); 1847 1848 CHECKFLOATVAL(result, false, true); 1849 PG_RETURN_FLOAT8(result); 1850 } 1851 1852 1853 /* 1854 * atan2 - returns the arctan of arg1/arg2 (radians) 1855 */ 1856 Datum 1857 datan2(PG_FUNCTION_ARGS) 1858 { 1859 float8 arg1 = PG_GETARG_FLOAT8(0); 1860 float8 arg2 = PG_GETARG_FLOAT8(1); 1861 float8 result; 1862 1863 /* Per the POSIX spec, return NaN if either input is NaN */ 1864 if (isnan(arg1) || isnan(arg2)) 1865 PG_RETURN_FLOAT8(get_float8_nan()); 1866 1867 /* 1868 * atan2 maps all inputs to values in the range [-Pi, Pi], so the result 1869 * should always be finite, even if the inputs are infinite. 1870 */ 1871 result = atan2(arg1, arg2); 1872 1873 CHECKFLOATVAL(result, false, true); 1874 PG_RETURN_FLOAT8(result); 1875 } 1876 1877 1878 /* 1879 * dcos - returns the cosine of arg1 (radians) 1880 */ 1881 Datum 1882 dcos(PG_FUNCTION_ARGS) 1883 { 1884 float8 arg1 = PG_GETARG_FLOAT8(0); 1885 float8 result; 1886 1887 /* Per the POSIX spec, return NaN if the input is NaN */ 1888 if (isnan(arg1)) 1889 PG_RETURN_FLOAT8(get_float8_nan()); 1890 1891 /* 1892 * cos() is periodic and so theoretically can work for all finite inputs, 1893 * but some implementations may choose to throw error if the input is so 1894 * large that there are no significant digits in the result. So we should 1895 * check for errors. POSIX allows an error to be reported either via 1896 * errno or via fetestexcept(), but currently we only support checking 1897 * errno. (fetestexcept() is rumored to report underflow unreasonably 1898 * early on some platforms, so it's not clear that believing it would be a 1899 * net improvement anyway.) 1900 * 1901 * For infinite inputs, POSIX specifies that the trigonometric functions 1902 * should return a domain error; but we won't notice that unless the 1903 * platform reports via errno, so also explicitly test for infinite 1904 * inputs. 1905 */ 1906 errno = 0; 1907 result = cos(arg1); 1908 if (errno != 0 || isinf(arg1)) 1909 ereport(ERROR, 1910 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 1911 errmsg("input is out of range"))); 1912 1913 CHECKFLOATVAL(result, false, true); 1914 PG_RETURN_FLOAT8(result); 1915 } 1916 1917 1918 /* 1919 * dcot - returns the cotangent of arg1 (radians) 1920 */ 1921 Datum 1922 dcot(PG_FUNCTION_ARGS) 1923 { 1924 float8 arg1 = PG_GETARG_FLOAT8(0); 1925 float8 result; 1926 1927 /* Per the POSIX spec, return NaN if the input is NaN */ 1928 if (isnan(arg1)) 1929 PG_RETURN_FLOAT8(get_float8_nan()); 1930 1931 /* Be sure to throw an error if the input is infinite --- see dcos() */ 1932 errno = 0; 1933 result = tan(arg1); 1934 if (errno != 0 || isinf(arg1)) 1935 ereport(ERROR, 1936 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 1937 errmsg("input is out of range"))); 1938 1939 result = 1.0 / result; 1940 CHECKFLOATVAL(result, true /* cot(0) == Inf */ , true); 1941 PG_RETURN_FLOAT8(result); 1942 } 1943 1944 1945 /* 1946 * dsin - returns the sine of arg1 (radians) 1947 */ 1948 Datum 1949 dsin(PG_FUNCTION_ARGS) 1950 { 1951 float8 arg1 = PG_GETARG_FLOAT8(0); 1952 float8 result; 1953 1954 /* Per the POSIX spec, return NaN if the input is NaN */ 1955 if (isnan(arg1)) 1956 PG_RETURN_FLOAT8(get_float8_nan()); 1957 1958 /* Be sure to throw an error if the input is infinite --- see dcos() */ 1959 errno = 0; 1960 result = sin(arg1); 1961 if (errno != 0 || isinf(arg1)) 1962 ereport(ERROR, 1963 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 1964 errmsg("input is out of range"))); 1965 1966 CHECKFLOATVAL(result, false, true); 1967 PG_RETURN_FLOAT8(result); 1968 } 1969 1970 1971 /* 1972 * dtan - returns the tangent of arg1 (radians) 1973 */ 1974 Datum 1975 dtan(PG_FUNCTION_ARGS) 1976 { 1977 float8 arg1 = PG_GETARG_FLOAT8(0); 1978 float8 result; 1979 1980 /* Per the POSIX spec, return NaN if the input is NaN */ 1981 if (isnan(arg1)) 1982 PG_RETURN_FLOAT8(get_float8_nan()); 1983 1984 /* Be sure to throw an error if the input is infinite --- see dcos() */ 1985 errno = 0; 1986 result = tan(arg1); 1987 if (errno != 0 || isinf(arg1)) 1988 ereport(ERROR, 1989 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 1990 errmsg("input is out of range"))); 1991 1992 CHECKFLOATVAL(result, true /* tan(pi/2) == Inf */ , true); 1993 PG_RETURN_FLOAT8(result); 1994 } 1995 1996 1997 /* ========== DEGREE-BASED TRIGONOMETRIC FUNCTIONS ========== */ 1998 1999 2000 /* 2001 * Initialize the cached constants declared at the head of this file 2002 * (sin_30 etc). The fact that we need those at all, let alone need this 2003 * Rube-Goldberg-worthy method of initializing them, is because there are 2004 * compilers out there that will precompute expressions such as sin(constant) 2005 * using a sin() function different from what will be used at runtime. If we 2006 * want exact results, we must ensure that none of the scaling constants used 2007 * in the degree-based trig functions are computed that way. To do so, we 2008 * compute them from the variables degree_c_thirty etc, which are also really 2009 * constants, but the compiler cannot assume that. 2010 * 2011 * Other hazards we are trying to forestall with this kluge include the 2012 * possibility that compilers will rearrange the expressions, or compute 2013 * some intermediate results in registers wider than a standard double. 2014 * 2015 * In the places where we use these constants, the typical pattern is like 2016 * volatile float8 sin_x = sin(x * RADIANS_PER_DEGREE); 2017 * return (sin_x / sin_30); 2018 * where we hope to get a value of exactly 1.0 from the division when x = 30. 2019 * The volatile temporary variable is needed on machines with wide float 2020 * registers, to ensure that the result of sin(x) is rounded to double width 2021 * the same as the value of sin_30 has been. Experimentation with gcc shows 2022 * that marking the temp variable volatile is necessary to make the store and 2023 * reload actually happen; hopefully the same trick works for other compilers. 2024 * (gcc's documentation suggests using the -ffloat-store compiler switch to 2025 * ensure this, but that is compiler-specific and it also pessimizes code in 2026 * many places where we don't care about this.) 2027 */ 2028 static void 2029 init_degree_constants(void) 2030 { 2031 sin_30 = sin(degree_c_thirty * RADIANS_PER_DEGREE); 2032 one_minus_cos_60 = 1.0 - cos(degree_c_sixty * RADIANS_PER_DEGREE); 2033 asin_0_5 = asin(degree_c_one_half); 2034 acos_0_5 = acos(degree_c_one_half); 2035 atan_1_0 = atan(degree_c_one); 2036 tan_45 = sind_q1(degree_c_forty_five) / cosd_q1(degree_c_forty_five); 2037 cot_45 = cosd_q1(degree_c_forty_five) / sind_q1(degree_c_forty_five); 2038 degree_consts_set = true; 2039 } 2040 2041 #define INIT_DEGREE_CONSTANTS() \ 2042 do { \ 2043 if (!degree_consts_set) \ 2044 init_degree_constants(); \ 2045 } while(0) 2046 2047 2048 /* 2049 * asind_q1 - returns the inverse sine of x in degrees, for x in 2050 * the range [0, 1]. The result is an angle in the 2051 * first quadrant --- [0, 90] degrees. 2052 * 2053 * For the 3 special case inputs (0, 0.5 and 1), this 2054 * function will return exact values (0, 30 and 90 2055 * degrees respectively). 2056 */ 2057 static double 2058 asind_q1(double x) 2059 { 2060 /* 2061 * Stitch together inverse sine and cosine functions for the ranges [0, 2062 * 0.5] and (0.5, 1]. Each expression below is guaranteed to return 2063 * exactly 30 for x=0.5, so the result is a continuous monotonic function 2064 * over the full range. 2065 */ 2066 if (x <= 0.5) 2067 { 2068 volatile float8 asin_x = asin(x); 2069 2070 return (asin_x / asin_0_5) * 30.0; 2071 } 2072 else 2073 { 2074 volatile float8 acos_x = acos(x); 2075 2076 return 90.0 - (acos_x / acos_0_5) * 60.0; 2077 } 2078 } 2079 2080 2081 /* 2082 * acosd_q1 - returns the inverse cosine of x in degrees, for x in 2083 * the range [0, 1]. The result is an angle in the 2084 * first quadrant --- [0, 90] degrees. 2085 * 2086 * For the 3 special case inputs (0, 0.5 and 1), this 2087 * function will return exact values (0, 60 and 90 2088 * degrees respectively). 2089 */ 2090 static double 2091 acosd_q1(double x) 2092 { 2093 /* 2094 * Stitch together inverse sine and cosine functions for the ranges [0, 2095 * 0.5] and (0.5, 1]. Each expression below is guaranteed to return 2096 * exactly 60 for x=0.5, so the result is a continuous monotonic function 2097 * over the full range. 2098 */ 2099 if (x <= 0.5) 2100 { 2101 volatile float8 asin_x = asin(x); 2102 2103 return 90.0 - (asin_x / asin_0_5) * 30.0; 2104 } 2105 else 2106 { 2107 volatile float8 acos_x = acos(x); 2108 2109 return (acos_x / acos_0_5) * 60.0; 2110 } 2111 } 2112 2113 2114 /* 2115 * dacosd - returns the arccos of arg1 (degrees) 2116 */ 2117 Datum 2118 dacosd(PG_FUNCTION_ARGS) 2119 { 2120 float8 arg1 = PG_GETARG_FLOAT8(0); 2121 float8 result; 2122 2123 /* Per the POSIX spec, return NaN if the input is NaN */ 2124 if (isnan(arg1)) 2125 PG_RETURN_FLOAT8(get_float8_nan()); 2126 2127 INIT_DEGREE_CONSTANTS(); 2128 2129 /* 2130 * The principal branch of the inverse cosine function maps values in the 2131 * range [-1, 1] to values in the range [0, 180], so we should reject any 2132 * inputs outside that range and the result will always be finite. 2133 */ 2134 if (arg1 < -1.0 || arg1 > 1.0) 2135 ereport(ERROR, 2136 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 2137 errmsg("input is out of range"))); 2138 2139 if (arg1 >= 0.0) 2140 result = acosd_q1(arg1); 2141 else 2142 result = 90.0 + asind_q1(-arg1); 2143 2144 CHECKFLOATVAL(result, false, true); 2145 PG_RETURN_FLOAT8(result); 2146 } 2147 2148 2149 /* 2150 * dasind - returns the arcsin of arg1 (degrees) 2151 */ 2152 Datum 2153 dasind(PG_FUNCTION_ARGS) 2154 { 2155 float8 arg1 = PG_GETARG_FLOAT8(0); 2156 float8 result; 2157 2158 /* Per the POSIX spec, return NaN if the input is NaN */ 2159 if (isnan(arg1)) 2160 PG_RETURN_FLOAT8(get_float8_nan()); 2161 2162 INIT_DEGREE_CONSTANTS(); 2163 2164 /* 2165 * The principal branch of the inverse sine function maps values in the 2166 * range [-1, 1] to values in the range [-90, 90], so we should reject any 2167 * inputs outside that range and the result will always be finite. 2168 */ 2169 if (arg1 < -1.0 || arg1 > 1.0) 2170 ereport(ERROR, 2171 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 2172 errmsg("input is out of range"))); 2173 2174 if (arg1 >= 0.0) 2175 result = asind_q1(arg1); 2176 else 2177 result = -asind_q1(-arg1); 2178 2179 CHECKFLOATVAL(result, false, true); 2180 PG_RETURN_FLOAT8(result); 2181 } 2182 2183 2184 /* 2185 * datand - returns the arctan of arg1 (degrees) 2186 */ 2187 Datum 2188 datand(PG_FUNCTION_ARGS) 2189 { 2190 float8 arg1 = PG_GETARG_FLOAT8(0); 2191 float8 result; 2192 volatile float8 atan_arg1; 2193 2194 /* Per the POSIX spec, return NaN if the input is NaN */ 2195 if (isnan(arg1)) 2196 PG_RETURN_FLOAT8(get_float8_nan()); 2197 2198 INIT_DEGREE_CONSTANTS(); 2199 2200 /* 2201 * The principal branch of the inverse tangent function maps all inputs to 2202 * values in the range [-90, 90], so the result should always be finite, 2203 * even if the input is infinite. Additionally, we take care to ensure 2204 * than when arg1 is 1, the result is exactly 45. 2205 */ 2206 atan_arg1 = atan(arg1); 2207 result = (atan_arg1 / atan_1_0) * 45.0; 2208 2209 CHECKFLOATVAL(result, false, true); 2210 PG_RETURN_FLOAT8(result); 2211 } 2212 2213 2214 /* 2215 * atan2d - returns the arctan of arg1/arg2 (degrees) 2216 */ 2217 Datum 2218 datan2d(PG_FUNCTION_ARGS) 2219 { 2220 float8 arg1 = PG_GETARG_FLOAT8(0); 2221 float8 arg2 = PG_GETARG_FLOAT8(1); 2222 float8 result; 2223 volatile float8 atan2_arg1_arg2; 2224 2225 /* Per the POSIX spec, return NaN if either input is NaN */ 2226 if (isnan(arg1) || isnan(arg2)) 2227 PG_RETURN_FLOAT8(get_float8_nan()); 2228 2229 INIT_DEGREE_CONSTANTS(); 2230 2231 /* 2232 * atan2d maps all inputs to values in the range [-180, 180], so the 2233 * result should always be finite, even if the inputs are infinite. 2234 * 2235 * Note: this coding assumes that atan(1.0) is a suitable scaling constant 2236 * to get an exact result from atan2(). This might well fail on us at 2237 * some point, requiring us to decide exactly what inputs we think we're 2238 * going to guarantee an exact result for. 2239 */ 2240 atan2_arg1_arg2 = atan2(arg1, arg2); 2241 result = (atan2_arg1_arg2 / atan_1_0) * 45.0; 2242 2243 CHECKFLOATVAL(result, false, true); 2244 PG_RETURN_FLOAT8(result); 2245 } 2246 2247 2248 /* 2249 * sind_0_to_30 - returns the sine of an angle that lies between 0 and 2250 * 30 degrees. This will return exactly 0 when x is 0, 2251 * and exactly 0.5 when x is 30 degrees. 2252 */ 2253 static double 2254 sind_0_to_30(double x) 2255 { 2256 volatile float8 sin_x = sin(x * RADIANS_PER_DEGREE); 2257 2258 return (sin_x / sin_30) / 2.0; 2259 } 2260 2261 2262 /* 2263 * cosd_0_to_60 - returns the cosine of an angle that lies between 0 2264 * and 60 degrees. This will return exactly 1 when x 2265 * is 0, and exactly 0.5 when x is 60 degrees. 2266 */ 2267 static double 2268 cosd_0_to_60(double x) 2269 { 2270 volatile float8 one_minus_cos_x = 1.0 - cos(x * RADIANS_PER_DEGREE); 2271 2272 return 1.0 - (one_minus_cos_x / one_minus_cos_60) / 2.0; 2273 } 2274 2275 2276 /* 2277 * sind_q1 - returns the sine of an angle in the first quadrant 2278 * (0 to 90 degrees). 2279 */ 2280 static double 2281 sind_q1(double x) 2282 { 2283 /* 2284 * Stitch together the sine and cosine functions for the ranges [0, 30] 2285 * and (30, 90]. These guarantee to return exact answers at their 2286 * endpoints, so the overall result is a continuous monotonic function 2287 * that gives exact results when x = 0, 30 and 90 degrees. 2288 */ 2289 if (x <= 30.0) 2290 return sind_0_to_30(x); 2291 else 2292 return cosd_0_to_60(90.0 - x); 2293 } 2294 2295 2296 /* 2297 * cosd_q1 - returns the cosine of an angle in the first quadrant 2298 * (0 to 90 degrees). 2299 */ 2300 static double 2301 cosd_q1(double x) 2302 { 2303 /* 2304 * Stitch together the sine and cosine functions for the ranges [0, 60] 2305 * and (60, 90]. These guarantee to return exact answers at their 2306 * endpoints, so the overall result is a continuous monotonic function 2307 * that gives exact results when x = 0, 60 and 90 degrees. 2308 */ 2309 if (x <= 60.0) 2310 return cosd_0_to_60(x); 2311 else 2312 return sind_0_to_30(90.0 - x); 2313 } 2314 2315 2316 /* 2317 * dcosd - returns the cosine of arg1 (degrees) 2318 */ 2319 Datum 2320 dcosd(PG_FUNCTION_ARGS) 2321 { 2322 float8 arg1 = PG_GETARG_FLOAT8(0); 2323 float8 result; 2324 int sign = 1; 2325 2326 /* 2327 * Per the POSIX spec, return NaN if the input is NaN and throw an error 2328 * if the input is infinite. 2329 */ 2330 if (isnan(arg1)) 2331 PG_RETURN_FLOAT8(get_float8_nan()); 2332 2333 if (isinf(arg1)) 2334 ereport(ERROR, 2335 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 2336 errmsg("input is out of range"))); 2337 2338 INIT_DEGREE_CONSTANTS(); 2339 2340 /* Reduce the range of the input to [0,90] degrees */ 2341 arg1 = fmod(arg1, 360.0); 2342 2343 if (arg1 < 0.0) 2344 { 2345 /* cosd(-x) = cosd(x) */ 2346 arg1 = -arg1; 2347 } 2348 2349 if (arg1 > 180.0) 2350 { 2351 /* cosd(360-x) = cosd(x) */ 2352 arg1 = 360.0 - arg1; 2353 } 2354 2355 if (arg1 > 90.0) 2356 { 2357 /* cosd(180-x) = -cosd(x) */ 2358 arg1 = 180.0 - arg1; 2359 sign = -sign; 2360 } 2361 2362 result = sign * cosd_q1(arg1); 2363 2364 CHECKFLOATVAL(result, false, true); 2365 PG_RETURN_FLOAT8(result); 2366 } 2367 2368 2369 /* 2370 * dcotd - returns the cotangent of arg1 (degrees) 2371 */ 2372 Datum 2373 dcotd(PG_FUNCTION_ARGS) 2374 { 2375 float8 arg1 = PG_GETARG_FLOAT8(0); 2376 float8 result; 2377 volatile float8 cot_arg1; 2378 int sign = 1; 2379 2380 /* 2381 * Per the POSIX spec, return NaN if the input is NaN and throw an error 2382 * if the input is infinite. 2383 */ 2384 if (isnan(arg1)) 2385 PG_RETURN_FLOAT8(get_float8_nan()); 2386 2387 if (isinf(arg1)) 2388 ereport(ERROR, 2389 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 2390 errmsg("input is out of range"))); 2391 2392 INIT_DEGREE_CONSTANTS(); 2393 2394 /* Reduce the range of the input to [0,90] degrees */ 2395 arg1 = fmod(arg1, 360.0); 2396 2397 if (arg1 < 0.0) 2398 { 2399 /* cotd(-x) = -cotd(x) */ 2400 arg1 = -arg1; 2401 sign = -sign; 2402 } 2403 2404 if (arg1 > 180.0) 2405 { 2406 /* cotd(360-x) = -cotd(x) */ 2407 arg1 = 360.0 - arg1; 2408 sign = -sign; 2409 } 2410 2411 if (arg1 > 90.0) 2412 { 2413 /* cotd(180-x) = -cotd(x) */ 2414 arg1 = 180.0 - arg1; 2415 sign = -sign; 2416 } 2417 2418 cot_arg1 = cosd_q1(arg1) / sind_q1(arg1); 2419 result = sign * (cot_arg1 / cot_45); 2420 2421 /* 2422 * On some machines we get cotd(270) = minus zero, but this isn't always 2423 * true. For portability, and because the user constituency for this 2424 * function probably doesn't want minus zero, force it to plain zero. 2425 */ 2426 if (result == 0.0) 2427 result = 0.0; 2428 2429 CHECKFLOATVAL(result, true /* cotd(0) == Inf */ , true); 2430 PG_RETURN_FLOAT8(result); 2431 } 2432 2433 2434 /* 2435 * dsind - returns the sine of arg1 (degrees) 2436 */ 2437 Datum 2438 dsind(PG_FUNCTION_ARGS) 2439 { 2440 float8 arg1 = PG_GETARG_FLOAT8(0); 2441 float8 result; 2442 int sign = 1; 2443 2444 /* 2445 * Per the POSIX spec, return NaN if the input is NaN and throw an error 2446 * if the input is infinite. 2447 */ 2448 if (isnan(arg1)) 2449 PG_RETURN_FLOAT8(get_float8_nan()); 2450 2451 if (isinf(arg1)) 2452 ereport(ERROR, 2453 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 2454 errmsg("input is out of range"))); 2455 2456 INIT_DEGREE_CONSTANTS(); 2457 2458 /* Reduce the range of the input to [0,90] degrees */ 2459 arg1 = fmod(arg1, 360.0); 2460 2461 if (arg1 < 0.0) 2462 { 2463 /* sind(-x) = -sind(x) */ 2464 arg1 = -arg1; 2465 sign = -sign; 2466 } 2467 2468 if (arg1 > 180.0) 2469 { 2470 /* sind(360-x) = -sind(x) */ 2471 arg1 = 360.0 - arg1; 2472 sign = -sign; 2473 } 2474 2475 if (arg1 > 90.0) 2476 { 2477 /* sind(180-x) = sind(x) */ 2478 arg1 = 180.0 - arg1; 2479 } 2480 2481 result = sign * sind_q1(arg1); 2482 2483 CHECKFLOATVAL(result, false, true); 2484 PG_RETURN_FLOAT8(result); 2485 } 2486 2487 2488 /* 2489 * dtand - returns the tangent of arg1 (degrees) 2490 */ 2491 Datum 2492 dtand(PG_FUNCTION_ARGS) 2493 { 2494 float8 arg1 = PG_GETARG_FLOAT8(0); 2495 float8 result; 2496 volatile float8 tan_arg1; 2497 int sign = 1; 2498 2499 /* 2500 * Per the POSIX spec, return NaN if the input is NaN and throw an error 2501 * if the input is infinite. 2502 */ 2503 if (isnan(arg1)) 2504 PG_RETURN_FLOAT8(get_float8_nan()); 2505 2506 if (isinf(arg1)) 2507 ereport(ERROR, 2508 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 2509 errmsg("input is out of range"))); 2510 2511 INIT_DEGREE_CONSTANTS(); 2512 2513 /* Reduce the range of the input to [0,90] degrees */ 2514 arg1 = fmod(arg1, 360.0); 2515 2516 if (arg1 < 0.0) 2517 { 2518 /* tand(-x) = -tand(x) */ 2519 arg1 = -arg1; 2520 sign = -sign; 2521 } 2522 2523 if (arg1 > 180.0) 2524 { 2525 /* tand(360-x) = -tand(x) */ 2526 arg1 = 360.0 - arg1; 2527 sign = -sign; 2528 } 2529 2530 if (arg1 > 90.0) 2531 { 2532 /* tand(180-x) = -tand(x) */ 2533 arg1 = 180.0 - arg1; 2534 sign = -sign; 2535 } 2536 2537 tan_arg1 = sind_q1(arg1) / cosd_q1(arg1); 2538 result = sign * (tan_arg1 / tan_45); 2539 2540 /* 2541 * On some machines we get tand(180) = minus zero, but this isn't always 2542 * true. For portability, and because the user constituency for this 2543 * function probably doesn't want minus zero, force it to plain zero. 2544 */ 2545 if (result == 0.0) 2546 result = 0.0; 2547 2548 CHECKFLOATVAL(result, true /* tand(90) == Inf */ , true); 2549 PG_RETURN_FLOAT8(result); 2550 } 2551 2552 2553 /* 2554 * degrees - returns degrees converted from radians 2555 */ 2556 Datum 2557 degrees(PG_FUNCTION_ARGS) 2558 { 2559 float8 arg1 = PG_GETARG_FLOAT8(0); 2560 float8 result; 2561 2562 result = arg1 / RADIANS_PER_DEGREE; 2563 2564 CHECKFLOATVAL(result, isinf(arg1), arg1 == 0); 2565 PG_RETURN_FLOAT8(result); 2566 } 2567 2568 2569 /* 2570 * dpi - returns the constant PI 2571 */ 2572 Datum 2573 dpi(PG_FUNCTION_ARGS) 2574 { 2575 PG_RETURN_FLOAT8(M_PI); 2576 } 2577 2578 2579 /* 2580 * radians - returns radians converted from degrees 2581 */ 2582 Datum 2583 radians(PG_FUNCTION_ARGS) 2584 { 2585 float8 arg1 = PG_GETARG_FLOAT8(0); 2586 float8 result; 2587 2588 result = arg1 * RADIANS_PER_DEGREE; 2589 2590 CHECKFLOATVAL(result, isinf(arg1), arg1 == 0); 2591 PG_RETURN_FLOAT8(result); 2592 } 2593 2594 2595 /* 2596 * drandom - returns a random number 2597 */ 2598 Datum 2599 drandom(PG_FUNCTION_ARGS) 2600 { 2601 float8 result; 2602 2603 /* result [0.0 - 1.0) */ 2604 result = (double) random() / ((double) MAX_RANDOM_VALUE + 1); 2605 2606 PG_RETURN_FLOAT8(result); 2607 } 2608 2609 2610 /* 2611 * setseed - set seed for the random number generator 2612 */ 2613 Datum 2614 setseed(PG_FUNCTION_ARGS) 2615 { 2616 float8 seed = PG_GETARG_FLOAT8(0); 2617 int iseed; 2618 2619 if (seed < -1 || seed > 1) 2620 elog(ERROR, "setseed parameter %f out of range [-1,1]", seed); 2621 2622 iseed = (int) (seed * MAX_RANDOM_VALUE); 2623 srandom((unsigned int) iseed); 2624 2625 PG_RETURN_VOID(); 2626 } 2627 2628 2629 2630 /* 2631 * ========================= 2632 * FLOAT AGGREGATE OPERATORS 2633 * ========================= 2634 * 2635 * float8_accum - accumulate for AVG(), variance aggregates, etc. 2636 * float4_accum - same, but input data is float4 2637 * float8_avg - produce final result for float AVG() 2638 * float8_var_samp - produce final result for float VAR_SAMP() 2639 * float8_var_pop - produce final result for float VAR_POP() 2640 * float8_stddev_samp - produce final result for float STDDEV_SAMP() 2641 * float8_stddev_pop - produce final result for float STDDEV_POP() 2642 * 2643 * The transition datatype for all these aggregates is a 3-element array 2644 * of float8, holding the values N, sum(X), sum(X*X) in that order. 2645 * 2646 * Note that we represent N as a float to avoid having to build a special 2647 * datatype. Given a reasonable floating-point implementation, there should 2648 * be no accuracy loss unless N exceeds 2 ^ 52 or so (by which time the 2649 * user will have doubtless lost interest anyway...) 2650 */ 2651 2652 static float8 * 2653 check_float8_array(ArrayType *transarray, const char *caller, int n) 2654 { 2655 /* 2656 * We expect the input to be an N-element float array; verify that. We 2657 * don't need to use deconstruct_array() since the array data is just 2658 * going to look like a C array of N float8 values. 2659 */ 2660 if (ARR_NDIM(transarray) != 1 || 2661 ARR_DIMS(transarray)[0] != n || 2662 ARR_HASNULL(transarray) || 2663 ARR_ELEMTYPE(transarray) != FLOAT8OID) 2664 elog(ERROR, "%s: expected %d-element float8 array", caller, n); 2665 return (float8 *) ARR_DATA_PTR(transarray); 2666 } 2667 2668 /* 2669 * float8_combine 2670 * 2671 * An aggregate combine function used to combine two 3 fields 2672 * aggregate transition data into a single transition data. 2673 * This function is used only in two stage aggregation and 2674 * shouldn't be called outside aggregate context. 2675 */ 2676 Datum 2677 float8_combine(PG_FUNCTION_ARGS) 2678 { 2679 ArrayType *transarray1 = PG_GETARG_ARRAYTYPE_P(0); 2680 ArrayType *transarray2 = PG_GETARG_ARRAYTYPE_P(1); 2681 float8 *transvalues1; 2682 float8 *transvalues2; 2683 float8 N, 2684 sumX, 2685 sumX2; 2686 2687 if (!AggCheckCallContext(fcinfo, NULL)) 2688 elog(ERROR, "aggregate function called in non-aggregate context"); 2689 2690 transvalues1 = check_float8_array(transarray1, "float8_combine", 3); 2691 N = transvalues1[0]; 2692 sumX = transvalues1[1]; 2693 sumX2 = transvalues1[2]; 2694 2695 transvalues2 = check_float8_array(transarray2, "float8_combine", 3); 2696 2697 N += transvalues2[0]; 2698 sumX += transvalues2[1]; 2699 CHECKFLOATVAL(sumX, isinf(transvalues1[1]) || isinf(transvalues2[1]), 2700 true); 2701 sumX2 += transvalues2[2]; 2702 CHECKFLOATVAL(sumX2, isinf(transvalues1[2]) || isinf(transvalues2[2]), 2703 true); 2704 2705 transvalues1[0] = N; 2706 transvalues1[1] = sumX; 2707 transvalues1[2] = sumX2; 2708 2709 PG_RETURN_ARRAYTYPE_P(transarray1); 2710 } 2711 2712 Datum 2713 float8_accum(PG_FUNCTION_ARGS) 2714 { 2715 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 2716 float8 newval = PG_GETARG_FLOAT8(1); 2717 float8 *transvalues; 2718 float8 N, 2719 sumX, 2720 sumX2; 2721 2722 transvalues = check_float8_array(transarray, "float8_accum", 3); 2723 N = transvalues[0]; 2724 sumX = transvalues[1]; 2725 sumX2 = transvalues[2]; 2726 2727 N += 1.0; 2728 sumX += newval; 2729 CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newval), true); 2730 sumX2 += newval * newval; 2731 CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newval), true); 2732 2733 /* 2734 * If we're invoked as an aggregate, we can cheat and modify our first 2735 * parameter in-place to reduce palloc overhead. Otherwise we construct a 2736 * new array with the updated transition data and return it. 2737 */ 2738 if (AggCheckCallContext(fcinfo, NULL)) 2739 { 2740 transvalues[0] = N; 2741 transvalues[1] = sumX; 2742 transvalues[2] = sumX2; 2743 2744 PG_RETURN_ARRAYTYPE_P(transarray); 2745 } 2746 else 2747 { 2748 Datum transdatums[3]; 2749 ArrayType *result; 2750 2751 transdatums[0] = Float8GetDatumFast(N); 2752 transdatums[1] = Float8GetDatumFast(sumX); 2753 transdatums[2] = Float8GetDatumFast(sumX2); 2754 2755 result = construct_array(transdatums, 3, 2756 FLOAT8OID, 2757 sizeof(float8), FLOAT8PASSBYVAL, 'd'); 2758 2759 PG_RETURN_ARRAYTYPE_P(result); 2760 } 2761 } 2762 2763 Datum 2764 float4_accum(PG_FUNCTION_ARGS) 2765 { 2766 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 2767 2768 /* do computations as float8 */ 2769 float8 newval = PG_GETARG_FLOAT4(1); 2770 float8 *transvalues; 2771 float8 N, 2772 sumX, 2773 sumX2; 2774 2775 transvalues = check_float8_array(transarray, "float4_accum", 3); 2776 N = transvalues[0]; 2777 sumX = transvalues[1]; 2778 sumX2 = transvalues[2]; 2779 2780 N += 1.0; 2781 sumX += newval; 2782 CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newval), true); 2783 sumX2 += newval * newval; 2784 CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newval), true); 2785 2786 /* 2787 * If we're invoked as an aggregate, we can cheat and modify our first 2788 * parameter in-place to reduce palloc overhead. Otherwise we construct a 2789 * new array with the updated transition data and return it. 2790 */ 2791 if (AggCheckCallContext(fcinfo, NULL)) 2792 { 2793 transvalues[0] = N; 2794 transvalues[1] = sumX; 2795 transvalues[2] = sumX2; 2796 2797 PG_RETURN_ARRAYTYPE_P(transarray); 2798 } 2799 else 2800 { 2801 Datum transdatums[3]; 2802 ArrayType *result; 2803 2804 transdatums[0] = Float8GetDatumFast(N); 2805 transdatums[1] = Float8GetDatumFast(sumX); 2806 transdatums[2] = Float8GetDatumFast(sumX2); 2807 2808 result = construct_array(transdatums, 3, 2809 FLOAT8OID, 2810 sizeof(float8), FLOAT8PASSBYVAL, 'd'); 2811 2812 PG_RETURN_ARRAYTYPE_P(result); 2813 } 2814 } 2815 2816 Datum 2817 float8_avg(PG_FUNCTION_ARGS) 2818 { 2819 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 2820 float8 *transvalues; 2821 float8 N, 2822 sumX; 2823 2824 transvalues = check_float8_array(transarray, "float8_avg", 3); 2825 N = transvalues[0]; 2826 sumX = transvalues[1]; 2827 /* ignore sumX2 */ 2828 2829 /* SQL defines AVG of no values to be NULL */ 2830 if (N == 0.0) 2831 PG_RETURN_NULL(); 2832 2833 PG_RETURN_FLOAT8(sumX / N); 2834 } 2835 2836 Datum 2837 float8_var_pop(PG_FUNCTION_ARGS) 2838 { 2839 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 2840 float8 *transvalues; 2841 float8 N, 2842 sumX, 2843 sumX2, 2844 numerator; 2845 2846 transvalues = check_float8_array(transarray, "float8_var_pop", 3); 2847 N = transvalues[0]; 2848 sumX = transvalues[1]; 2849 sumX2 = transvalues[2]; 2850 2851 /* Population variance is undefined when N is 0, so return NULL */ 2852 if (N == 0.0) 2853 PG_RETURN_NULL(); 2854 2855 numerator = N * sumX2 - sumX * sumX; 2856 CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true); 2857 2858 /* Watch out for roundoff error producing a negative numerator */ 2859 if (numerator <= 0.0) 2860 PG_RETURN_FLOAT8(0.0); 2861 2862 PG_RETURN_FLOAT8(numerator / (N * N)); 2863 } 2864 2865 Datum 2866 float8_var_samp(PG_FUNCTION_ARGS) 2867 { 2868 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 2869 float8 *transvalues; 2870 float8 N, 2871 sumX, 2872 sumX2, 2873 numerator; 2874 2875 transvalues = check_float8_array(transarray, "float8_var_samp", 3); 2876 N = transvalues[0]; 2877 sumX = transvalues[1]; 2878 sumX2 = transvalues[2]; 2879 2880 /* Sample variance is undefined when N is 0 or 1, so return NULL */ 2881 if (N <= 1.0) 2882 PG_RETURN_NULL(); 2883 2884 numerator = N * sumX2 - sumX * sumX; 2885 CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true); 2886 2887 /* Watch out for roundoff error producing a negative numerator */ 2888 if (numerator <= 0.0) 2889 PG_RETURN_FLOAT8(0.0); 2890 2891 PG_RETURN_FLOAT8(numerator / (N * (N - 1.0))); 2892 } 2893 2894 Datum 2895 float8_stddev_pop(PG_FUNCTION_ARGS) 2896 { 2897 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 2898 float8 *transvalues; 2899 float8 N, 2900 sumX, 2901 sumX2, 2902 numerator; 2903 2904 transvalues = check_float8_array(transarray, "float8_stddev_pop", 3); 2905 N = transvalues[0]; 2906 sumX = transvalues[1]; 2907 sumX2 = transvalues[2]; 2908 2909 /* Population stddev is undefined when N is 0, so return NULL */ 2910 if (N == 0.0) 2911 PG_RETURN_NULL(); 2912 2913 numerator = N * sumX2 - sumX * sumX; 2914 CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true); 2915 2916 /* Watch out for roundoff error producing a negative numerator */ 2917 if (numerator <= 0.0) 2918 PG_RETURN_FLOAT8(0.0); 2919 2920 PG_RETURN_FLOAT8(sqrt(numerator / (N * N))); 2921 } 2922 2923 Datum 2924 float8_stddev_samp(PG_FUNCTION_ARGS) 2925 { 2926 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 2927 float8 *transvalues; 2928 float8 N, 2929 sumX, 2930 sumX2, 2931 numerator; 2932 2933 transvalues = check_float8_array(transarray, "float8_stddev_samp", 3); 2934 N = transvalues[0]; 2935 sumX = transvalues[1]; 2936 sumX2 = transvalues[2]; 2937 2938 /* Sample stddev is undefined when N is 0 or 1, so return NULL */ 2939 if (N <= 1.0) 2940 PG_RETURN_NULL(); 2941 2942 numerator = N * sumX2 - sumX * sumX; 2943 CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true); 2944 2945 /* Watch out for roundoff error producing a negative numerator */ 2946 if (numerator <= 0.0) 2947 PG_RETURN_FLOAT8(0.0); 2948 2949 PG_RETURN_FLOAT8(sqrt(numerator / (N * (N - 1.0)))); 2950 } 2951 2952 /* 2953 * ========================= 2954 * SQL2003 BINARY AGGREGATES 2955 * ========================= 2956 * 2957 * The transition datatype for all these aggregates is a 6-element array of 2958 * float8, holding the values N, sum(X), sum(X*X), sum(Y), sum(Y*Y), sum(X*Y) 2959 * in that order. Note that Y is the first argument to the aggregates! 2960 * 2961 * It might seem attractive to optimize this by having multiple accumulator 2962 * functions that only calculate the sums actually needed. But on most 2963 * modern machines, a couple of extra floating-point multiplies will be 2964 * insignificant compared to the other per-tuple overhead, so I've chosen 2965 * to minimize code space instead. 2966 */ 2967 2968 Datum 2969 float8_regr_accum(PG_FUNCTION_ARGS) 2970 { 2971 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 2972 float8 newvalY = PG_GETARG_FLOAT8(1); 2973 float8 newvalX = PG_GETARG_FLOAT8(2); 2974 float8 *transvalues; 2975 float8 N, 2976 sumX, 2977 sumX2, 2978 sumY, 2979 sumY2, 2980 sumXY; 2981 2982 transvalues = check_float8_array(transarray, "float8_regr_accum", 6); 2983 N = transvalues[0]; 2984 sumX = transvalues[1]; 2985 sumX2 = transvalues[2]; 2986 sumY = transvalues[3]; 2987 sumY2 = transvalues[4]; 2988 sumXY = transvalues[5]; 2989 2990 N += 1.0; 2991 sumX += newvalX; 2992 CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newvalX), true); 2993 sumX2 += newvalX * newvalX; 2994 CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newvalX), true); 2995 sumY += newvalY; 2996 CHECKFLOATVAL(sumY, isinf(transvalues[3]) || isinf(newvalY), true); 2997 sumY2 += newvalY * newvalY; 2998 CHECKFLOATVAL(sumY2, isinf(transvalues[4]) || isinf(newvalY), true); 2999 sumXY += newvalX * newvalY; 3000 CHECKFLOATVAL(sumXY, isinf(transvalues[5]) || isinf(newvalX) || 3001 isinf(newvalY), true); 3002 3003 /* 3004 * If we're invoked as an aggregate, we can cheat and modify our first 3005 * parameter in-place to reduce palloc overhead. Otherwise we construct a 3006 * new array with the updated transition data and return it. 3007 */ 3008 if (AggCheckCallContext(fcinfo, NULL)) 3009 { 3010 transvalues[0] = N; 3011 transvalues[1] = sumX; 3012 transvalues[2] = sumX2; 3013 transvalues[3] = sumY; 3014 transvalues[4] = sumY2; 3015 transvalues[5] = sumXY; 3016 3017 PG_RETURN_ARRAYTYPE_P(transarray); 3018 } 3019 else 3020 { 3021 Datum transdatums[6]; 3022 ArrayType *result; 3023 3024 transdatums[0] = Float8GetDatumFast(N); 3025 transdatums[1] = Float8GetDatumFast(sumX); 3026 transdatums[2] = Float8GetDatumFast(sumX2); 3027 transdatums[3] = Float8GetDatumFast(sumY); 3028 transdatums[4] = Float8GetDatumFast(sumY2); 3029 transdatums[5] = Float8GetDatumFast(sumXY); 3030 3031 result = construct_array(transdatums, 6, 3032 FLOAT8OID, 3033 sizeof(float8), FLOAT8PASSBYVAL, 'd'); 3034 3035 PG_RETURN_ARRAYTYPE_P(result); 3036 } 3037 } 3038 3039 /* 3040 * float8_regr_combine 3041 * 3042 * An aggregate combine function used to combine two 6 fields 3043 * aggregate transition data into a single transition data. 3044 * This function is used only in two stage aggregation and 3045 * shouldn't be called outside aggregate context. 3046 */ 3047 Datum 3048 float8_regr_combine(PG_FUNCTION_ARGS) 3049 { 3050 ArrayType *transarray1 = PG_GETARG_ARRAYTYPE_P(0); 3051 ArrayType *transarray2 = PG_GETARG_ARRAYTYPE_P(1); 3052 float8 *transvalues1; 3053 float8 *transvalues2; 3054 float8 N, 3055 sumX, 3056 sumX2, 3057 sumY, 3058 sumY2, 3059 sumXY; 3060 3061 if (!AggCheckCallContext(fcinfo, NULL)) 3062 elog(ERROR, "aggregate function called in non-aggregate context"); 3063 3064 transvalues1 = check_float8_array(transarray1, "float8_regr_combine", 6); 3065 N = transvalues1[0]; 3066 sumX = transvalues1[1]; 3067 sumX2 = transvalues1[2]; 3068 sumY = transvalues1[3]; 3069 sumY2 = transvalues1[4]; 3070 sumXY = transvalues1[5]; 3071 3072 transvalues2 = check_float8_array(transarray2, "float8_regr_combine", 6); 3073 3074 N += transvalues2[0]; 3075 sumX += transvalues2[1]; 3076 CHECKFLOATVAL(sumX, isinf(transvalues1[1]) || isinf(transvalues2[1]), 3077 true); 3078 sumX2 += transvalues2[2]; 3079 CHECKFLOATVAL(sumX2, isinf(transvalues1[2]) || isinf(transvalues2[2]), 3080 true); 3081 sumY += transvalues2[3]; 3082 CHECKFLOATVAL(sumY, isinf(transvalues1[3]) || isinf(transvalues2[3]), 3083 true); 3084 sumY2 += transvalues2[4]; 3085 CHECKFLOATVAL(sumY2, isinf(transvalues1[4]) || isinf(transvalues2[4]), 3086 true); 3087 sumXY += transvalues2[5]; 3088 CHECKFLOATVAL(sumXY, isinf(transvalues1[5]) || isinf(transvalues2[5]), 3089 true); 3090 3091 transvalues1[0] = N; 3092 transvalues1[1] = sumX; 3093 transvalues1[2] = sumX2; 3094 transvalues1[3] = sumY; 3095 transvalues1[4] = sumY2; 3096 transvalues1[5] = sumXY; 3097 3098 PG_RETURN_ARRAYTYPE_P(transarray1); 3099 } 3100 3101 3102 Datum 3103 float8_regr_sxx(PG_FUNCTION_ARGS) 3104 { 3105 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 3106 float8 *transvalues; 3107 float8 N, 3108 sumX, 3109 sumX2, 3110 numerator; 3111 3112 transvalues = check_float8_array(transarray, "float8_regr_sxx", 6); 3113 N = transvalues[0]; 3114 sumX = transvalues[1]; 3115 sumX2 = transvalues[2]; 3116 3117 /* if N is 0 we should return NULL */ 3118 if (N < 1.0) 3119 PG_RETURN_NULL(); 3120 3121 numerator = N * sumX2 - sumX * sumX; 3122 CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true); 3123 3124 /* Watch out for roundoff error producing a negative numerator */ 3125 if (numerator <= 0.0) 3126 PG_RETURN_FLOAT8(0.0); 3127 3128 PG_RETURN_FLOAT8(numerator / N); 3129 } 3130 3131 Datum 3132 float8_regr_syy(PG_FUNCTION_ARGS) 3133 { 3134 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 3135 float8 *transvalues; 3136 float8 N, 3137 sumY, 3138 sumY2, 3139 numerator; 3140 3141 transvalues = check_float8_array(transarray, "float8_regr_syy", 6); 3142 N = transvalues[0]; 3143 sumY = transvalues[3]; 3144 sumY2 = transvalues[4]; 3145 3146 /* if N is 0 we should return NULL */ 3147 if (N < 1.0) 3148 PG_RETURN_NULL(); 3149 3150 numerator = N * sumY2 - sumY * sumY; 3151 CHECKFLOATVAL(numerator, isinf(sumY2) || isinf(sumY), true); 3152 3153 /* Watch out for roundoff error producing a negative numerator */ 3154 if (numerator <= 0.0) 3155 PG_RETURN_FLOAT8(0.0); 3156 3157 PG_RETURN_FLOAT8(numerator / N); 3158 } 3159 3160 Datum 3161 float8_regr_sxy(PG_FUNCTION_ARGS) 3162 { 3163 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 3164 float8 *transvalues; 3165 float8 N, 3166 sumX, 3167 sumY, 3168 sumXY, 3169 numerator; 3170 3171 transvalues = check_float8_array(transarray, "float8_regr_sxy", 6); 3172 N = transvalues[0]; 3173 sumX = transvalues[1]; 3174 sumY = transvalues[3]; 3175 sumXY = transvalues[5]; 3176 3177 /* if N is 0 we should return NULL */ 3178 if (N < 1.0) 3179 PG_RETURN_NULL(); 3180 3181 numerator = N * sumXY - sumX * sumY; 3182 CHECKFLOATVAL(numerator, isinf(sumXY) || isinf(sumX) || 3183 isinf(sumY), true); 3184 3185 /* A negative result is valid here */ 3186 3187 PG_RETURN_FLOAT8(numerator / N); 3188 } 3189 3190 Datum 3191 float8_regr_avgx(PG_FUNCTION_ARGS) 3192 { 3193 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 3194 float8 *transvalues; 3195 float8 N, 3196 sumX; 3197 3198 transvalues = check_float8_array(transarray, "float8_regr_avgx", 6); 3199 N = transvalues[0]; 3200 sumX = transvalues[1]; 3201 3202 /* if N is 0 we should return NULL */ 3203 if (N < 1.0) 3204 PG_RETURN_NULL(); 3205 3206 PG_RETURN_FLOAT8(sumX / N); 3207 } 3208 3209 Datum 3210 float8_regr_avgy(PG_FUNCTION_ARGS) 3211 { 3212 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 3213 float8 *transvalues; 3214 float8 N, 3215 sumY; 3216 3217 transvalues = check_float8_array(transarray, "float8_regr_avgy", 6); 3218 N = transvalues[0]; 3219 sumY = transvalues[3]; 3220 3221 /* if N is 0 we should return NULL */ 3222 if (N < 1.0) 3223 PG_RETURN_NULL(); 3224 3225 PG_RETURN_FLOAT8(sumY / N); 3226 } 3227 3228 Datum 3229 float8_covar_pop(PG_FUNCTION_ARGS) 3230 { 3231 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 3232 float8 *transvalues; 3233 float8 N, 3234 sumX, 3235 sumY, 3236 sumXY, 3237 numerator; 3238 3239 transvalues = check_float8_array(transarray, "float8_covar_pop", 6); 3240 N = transvalues[0]; 3241 sumX = transvalues[1]; 3242 sumY = transvalues[3]; 3243 sumXY = transvalues[5]; 3244 3245 /* if N is 0 we should return NULL */ 3246 if (N < 1.0) 3247 PG_RETURN_NULL(); 3248 3249 numerator = N * sumXY - sumX * sumY; 3250 CHECKFLOATVAL(numerator, isinf(sumXY) || isinf(sumX) || 3251 isinf(sumY), true); 3252 3253 PG_RETURN_FLOAT8(numerator / (N * N)); 3254 } 3255 3256 Datum 3257 float8_covar_samp(PG_FUNCTION_ARGS) 3258 { 3259 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 3260 float8 *transvalues; 3261 float8 N, 3262 sumX, 3263 sumY, 3264 sumXY, 3265 numerator; 3266 3267 transvalues = check_float8_array(transarray, "float8_covar_samp", 6); 3268 N = transvalues[0]; 3269 sumX = transvalues[1]; 3270 sumY = transvalues[3]; 3271 sumXY = transvalues[5]; 3272 3273 /* if N is <= 1 we should return NULL */ 3274 if (N < 2.0) 3275 PG_RETURN_NULL(); 3276 3277 numerator = N * sumXY - sumX * sumY; 3278 CHECKFLOATVAL(numerator, isinf(sumXY) || isinf(sumX) || 3279 isinf(sumY), true); 3280 3281 PG_RETURN_FLOAT8(numerator / (N * (N - 1.0))); 3282 } 3283 3284 Datum 3285 float8_corr(PG_FUNCTION_ARGS) 3286 { 3287 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 3288 float8 *transvalues; 3289 float8 N, 3290 sumX, 3291 sumX2, 3292 sumY, 3293 sumY2, 3294 sumXY, 3295 numeratorX, 3296 numeratorY, 3297 numeratorXY; 3298 3299 transvalues = check_float8_array(transarray, "float8_corr", 6); 3300 N = transvalues[0]; 3301 sumX = transvalues[1]; 3302 sumX2 = transvalues[2]; 3303 sumY = transvalues[3]; 3304 sumY2 = transvalues[4]; 3305 sumXY = transvalues[5]; 3306 3307 /* if N is 0 we should return NULL */ 3308 if (N < 1.0) 3309 PG_RETURN_NULL(); 3310 3311 numeratorX = N * sumX2 - sumX * sumX; 3312 CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true); 3313 numeratorY = N * sumY2 - sumY * sumY; 3314 CHECKFLOATVAL(numeratorY, isinf(sumY2) || isinf(sumY), true); 3315 numeratorXY = N * sumXY - sumX * sumY; 3316 CHECKFLOATVAL(numeratorXY, isinf(sumXY) || isinf(sumX) || 3317 isinf(sumY), true); 3318 if (numeratorX <= 0 || numeratorY <= 0) 3319 PG_RETURN_NULL(); 3320 3321 PG_RETURN_FLOAT8(numeratorXY / sqrt(numeratorX * numeratorY)); 3322 } 3323 3324 Datum 3325 float8_regr_r2(PG_FUNCTION_ARGS) 3326 { 3327 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 3328 float8 *transvalues; 3329 float8 N, 3330 sumX, 3331 sumX2, 3332 sumY, 3333 sumY2, 3334 sumXY, 3335 numeratorX, 3336 numeratorY, 3337 numeratorXY; 3338 3339 transvalues = check_float8_array(transarray, "float8_regr_r2", 6); 3340 N = transvalues[0]; 3341 sumX = transvalues[1]; 3342 sumX2 = transvalues[2]; 3343 sumY = transvalues[3]; 3344 sumY2 = transvalues[4]; 3345 sumXY = transvalues[5]; 3346 3347 /* if N is 0 we should return NULL */ 3348 if (N < 1.0) 3349 PG_RETURN_NULL(); 3350 3351 numeratorX = N * sumX2 - sumX * sumX; 3352 CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true); 3353 numeratorY = N * sumY2 - sumY * sumY; 3354 CHECKFLOATVAL(numeratorY, isinf(sumY2) || isinf(sumY), true); 3355 numeratorXY = N * sumXY - sumX * sumY; 3356 CHECKFLOATVAL(numeratorXY, isinf(sumXY) || isinf(sumX) || 3357 isinf(sumY), true); 3358 if (numeratorX <= 0) 3359 PG_RETURN_NULL(); 3360 /* per spec, horizontal line produces 1.0 */ 3361 if (numeratorY <= 0) 3362 PG_RETURN_FLOAT8(1.0); 3363 3364 PG_RETURN_FLOAT8((numeratorXY * numeratorXY) / 3365 (numeratorX * numeratorY)); 3366 } 3367 3368 Datum 3369 float8_regr_slope(PG_FUNCTION_ARGS) 3370 { 3371 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 3372 float8 *transvalues; 3373 float8 N, 3374 sumX, 3375 sumX2, 3376 sumY, 3377 sumXY, 3378 numeratorX, 3379 numeratorXY; 3380 3381 transvalues = check_float8_array(transarray, "float8_regr_slope", 6); 3382 N = transvalues[0]; 3383 sumX = transvalues[1]; 3384 sumX2 = transvalues[2]; 3385 sumY = transvalues[3]; 3386 sumXY = transvalues[5]; 3387 3388 /* if N is 0 we should return NULL */ 3389 if (N < 1.0) 3390 PG_RETURN_NULL(); 3391 3392 numeratorX = N * sumX2 - sumX * sumX; 3393 CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true); 3394 numeratorXY = N * sumXY - sumX * sumY; 3395 CHECKFLOATVAL(numeratorXY, isinf(sumXY) || isinf(sumX) || 3396 isinf(sumY), true); 3397 if (numeratorX <= 0) 3398 PG_RETURN_NULL(); 3399 3400 PG_RETURN_FLOAT8(numeratorXY / numeratorX); 3401 } 3402 3403 Datum 3404 float8_regr_intercept(PG_FUNCTION_ARGS) 3405 { 3406 ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); 3407 float8 *transvalues; 3408 float8 N, 3409 sumX, 3410 sumX2, 3411 sumY, 3412 sumXY, 3413 numeratorX, 3414 numeratorXXY; 3415 3416 transvalues = check_float8_array(transarray, "float8_regr_intercept", 6); 3417 N = transvalues[0]; 3418 sumX = transvalues[1]; 3419 sumX2 = transvalues[2]; 3420 sumY = transvalues[3]; 3421 sumXY = transvalues[5]; 3422 3423 /* if N is 0 we should return NULL */ 3424 if (N < 1.0) 3425 PG_RETURN_NULL(); 3426 3427 numeratorX = N * sumX2 - sumX * sumX; 3428 CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true); 3429 numeratorXXY = sumY * sumX2 - sumX * sumXY; 3430 CHECKFLOATVAL(numeratorXXY, isinf(sumY) || isinf(sumX2) || 3431 isinf(sumX) || isinf(sumXY), true); 3432 if (numeratorX <= 0) 3433 PG_RETURN_NULL(); 3434 3435 PG_RETURN_FLOAT8(numeratorXXY / numeratorX); 3436 } 3437 3438 3439 /* 3440 * ==================================== 3441 * MIXED-PRECISION ARITHMETIC OPERATORS 3442 * ==================================== 3443 */ 3444 3445 /* 3446 * float48pl - returns arg1 + arg2 3447 * float48mi - returns arg1 - arg2 3448 * float48mul - returns arg1 * arg2 3449 * float48div - returns arg1 / arg2 3450 */ 3451 Datum 3452 float48pl(PG_FUNCTION_ARGS) 3453 { 3454 float4 arg1 = PG_GETARG_FLOAT4(0); 3455 float8 arg2 = PG_GETARG_FLOAT8(1); 3456 float8 result; 3457 3458 result = arg1 + arg2; 3459 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true); 3460 PG_RETURN_FLOAT8(result); 3461 } 3462 3463 Datum 3464 float48mi(PG_FUNCTION_ARGS) 3465 { 3466 float4 arg1 = PG_GETARG_FLOAT4(0); 3467 float8 arg2 = PG_GETARG_FLOAT8(1); 3468 float8 result; 3469 3470 result = arg1 - arg2; 3471 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true); 3472 PG_RETURN_FLOAT8(result); 3473 } 3474 3475 Datum 3476 float48mul(PG_FUNCTION_ARGS) 3477 { 3478 float4 arg1 = PG_GETARG_FLOAT4(0); 3479 float8 arg2 = PG_GETARG_FLOAT8(1); 3480 float8 result; 3481 3482 result = arg1 * arg2; 3483 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), 3484 arg1 == 0 || arg2 == 0); 3485 PG_RETURN_FLOAT8(result); 3486 } 3487 3488 Datum 3489 float48div(PG_FUNCTION_ARGS) 3490 { 3491 float4 arg1 = PG_GETARG_FLOAT4(0); 3492 float8 arg2 = PG_GETARG_FLOAT8(1); 3493 float8 result; 3494 3495 if (arg2 == 0.0) 3496 ereport(ERROR, 3497 (errcode(ERRCODE_DIVISION_BY_ZERO), 3498 errmsg("division by zero"))); 3499 3500 result = arg1 / arg2; 3501 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0); 3502 PG_RETURN_FLOAT8(result); 3503 } 3504 3505 /* 3506 * float84pl - returns arg1 + arg2 3507 * float84mi - returns arg1 - arg2 3508 * float84mul - returns arg1 * arg2 3509 * float84div - returns arg1 / arg2 3510 */ 3511 Datum 3512 float84pl(PG_FUNCTION_ARGS) 3513 { 3514 float8 arg1 = PG_GETARG_FLOAT8(0); 3515 float4 arg2 = PG_GETARG_FLOAT4(1); 3516 float8 result; 3517 3518 result = arg1 + arg2; 3519 3520 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true); 3521 PG_RETURN_FLOAT8(result); 3522 } 3523 3524 Datum 3525 float84mi(PG_FUNCTION_ARGS) 3526 { 3527 float8 arg1 = PG_GETARG_FLOAT8(0); 3528 float4 arg2 = PG_GETARG_FLOAT4(1); 3529 float8 result; 3530 3531 result = arg1 - arg2; 3532 3533 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true); 3534 PG_RETURN_FLOAT8(result); 3535 } 3536 3537 Datum 3538 float84mul(PG_FUNCTION_ARGS) 3539 { 3540 float8 arg1 = PG_GETARG_FLOAT8(0); 3541 float4 arg2 = PG_GETARG_FLOAT4(1); 3542 float8 result; 3543 3544 result = arg1 * arg2; 3545 3546 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), 3547 arg1 == 0 || arg2 == 0); 3548 PG_RETURN_FLOAT8(result); 3549 } 3550 3551 Datum 3552 float84div(PG_FUNCTION_ARGS) 3553 { 3554 float8 arg1 = PG_GETARG_FLOAT8(0); 3555 float4 arg2 = PG_GETARG_FLOAT4(1); 3556 float8 result; 3557 3558 if (arg2 == 0.0) 3559 ereport(ERROR, 3560 (errcode(ERRCODE_DIVISION_BY_ZERO), 3561 errmsg("division by zero"))); 3562 3563 result = arg1 / arg2; 3564 3565 CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0); 3566 PG_RETURN_FLOAT8(result); 3567 } 3568 3569 /* 3570 * ==================== 3571 * COMPARISON OPERATORS 3572 * ==================== 3573 */ 3574 3575 /* 3576 * float48{eq,ne,lt,le,gt,ge} - float4/float8 comparison operations 3577 */ 3578 Datum 3579 float48eq(PG_FUNCTION_ARGS) 3580 { 3581 float4 arg1 = PG_GETARG_FLOAT4(0); 3582 float8 arg2 = PG_GETARG_FLOAT8(1); 3583 3584 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0); 3585 } 3586 3587 Datum 3588 float48ne(PG_FUNCTION_ARGS) 3589 { 3590 float4 arg1 = PG_GETARG_FLOAT4(0); 3591 float8 arg2 = PG_GETARG_FLOAT8(1); 3592 3593 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0); 3594 } 3595 3596 Datum 3597 float48lt(PG_FUNCTION_ARGS) 3598 { 3599 float4 arg1 = PG_GETARG_FLOAT4(0); 3600 float8 arg2 = PG_GETARG_FLOAT8(1); 3601 3602 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0); 3603 } 3604 3605 Datum 3606 float48le(PG_FUNCTION_ARGS) 3607 { 3608 float4 arg1 = PG_GETARG_FLOAT4(0); 3609 float8 arg2 = PG_GETARG_FLOAT8(1); 3610 3611 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0); 3612 } 3613 3614 Datum 3615 float48gt(PG_FUNCTION_ARGS) 3616 { 3617 float4 arg1 = PG_GETARG_FLOAT4(0); 3618 float8 arg2 = PG_GETARG_FLOAT8(1); 3619 3620 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0); 3621 } 3622 3623 Datum 3624 float48ge(PG_FUNCTION_ARGS) 3625 { 3626 float4 arg1 = PG_GETARG_FLOAT4(0); 3627 float8 arg2 = PG_GETARG_FLOAT8(1); 3628 3629 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0); 3630 } 3631 3632 /* 3633 * float84{eq,ne,lt,le,gt,ge} - float8/float4 comparison operations 3634 */ 3635 Datum 3636 float84eq(PG_FUNCTION_ARGS) 3637 { 3638 float8 arg1 = PG_GETARG_FLOAT8(0); 3639 float4 arg2 = PG_GETARG_FLOAT4(1); 3640 3641 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0); 3642 } 3643 3644 Datum 3645 float84ne(PG_FUNCTION_ARGS) 3646 { 3647 float8 arg1 = PG_GETARG_FLOAT8(0); 3648 float4 arg2 = PG_GETARG_FLOAT4(1); 3649 3650 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0); 3651 } 3652 3653 Datum 3654 float84lt(PG_FUNCTION_ARGS) 3655 { 3656 float8 arg1 = PG_GETARG_FLOAT8(0); 3657 float4 arg2 = PG_GETARG_FLOAT4(1); 3658 3659 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0); 3660 } 3661 3662 Datum 3663 float84le(PG_FUNCTION_ARGS) 3664 { 3665 float8 arg1 = PG_GETARG_FLOAT8(0); 3666 float4 arg2 = PG_GETARG_FLOAT4(1); 3667 3668 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0); 3669 } 3670 3671 Datum 3672 float84gt(PG_FUNCTION_ARGS) 3673 { 3674 float8 arg1 = PG_GETARG_FLOAT8(0); 3675 float4 arg2 = PG_GETARG_FLOAT4(1); 3676 3677 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0); 3678 } 3679 3680 Datum 3681 float84ge(PG_FUNCTION_ARGS) 3682 { 3683 float8 arg1 = PG_GETARG_FLOAT8(0); 3684 float4 arg2 = PG_GETARG_FLOAT4(1); 3685 3686 PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0); 3687 } 3688 3689 /* 3690 * Implements the float8 version of the width_bucket() function 3691 * defined by SQL2003. See also width_bucket_numeric(). 3692 * 3693 * 'bound1' and 'bound2' are the lower and upper bounds of the 3694 * histogram's range, respectively. 'count' is the number of buckets 3695 * in the histogram. width_bucket() returns an integer indicating the 3696 * bucket number that 'operand' belongs to in an equiwidth histogram 3697 * with the specified characteristics. An operand smaller than the 3698 * lower bound is assigned to bucket 0. An operand greater than the 3699 * upper bound is assigned to an additional bucket (with number 3700 * count+1). We don't allow "NaN" for any of the float8 inputs, and we 3701 * don't allow either of the histogram bounds to be +/- infinity. 3702 */ 3703 Datum 3704 width_bucket_float8(PG_FUNCTION_ARGS) 3705 { 3706 float8 operand = PG_GETARG_FLOAT8(0); 3707 float8 bound1 = PG_GETARG_FLOAT8(1); 3708 float8 bound2 = PG_GETARG_FLOAT8(2); 3709 int32 count = PG_GETARG_INT32(3); 3710 int32 result; 3711 3712 if (count <= 0.0) 3713 ereport(ERROR, 3714 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION), 3715 errmsg("count must be greater than zero"))); 3716 3717 if (isnan(operand) || isnan(bound1) || isnan(bound2)) 3718 ereport(ERROR, 3719 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION), 3720 errmsg("operand, lower bound, and upper bound cannot be NaN"))); 3721 3722 /* Note that we allow "operand" to be infinite */ 3723 if (isinf(bound1) || isinf(bound2)) 3724 ereport(ERROR, 3725 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION), 3726 errmsg("lower and upper bounds must be finite"))); 3727 3728 if (bound1 < bound2) 3729 { 3730 if (operand < bound1) 3731 result = 0; 3732 else if (operand >= bound2) 3733 { 3734 if (pg_add_s32_overflow(count, 1, &result)) 3735 ereport(ERROR, 3736 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 3737 errmsg("integer out of range"))); 3738 } 3739 else 3740 result = ((float8) count * (operand - bound1) / (bound2 - bound1)) + 1; 3741 } 3742 else if (bound1 > bound2) 3743 { 3744 if (operand > bound1) 3745 result = 0; 3746 else if (operand <= bound2) 3747 { 3748 if (pg_add_s32_overflow(count, 1, &result)) 3749 ereport(ERROR, 3750 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), 3751 errmsg("integer out of range"))); 3752 } 3753 else 3754 result = ((float8) count * (bound1 - operand) / (bound1 - bound2)) + 1; 3755 } 3756 else 3757 { 3758 ereport(ERROR, 3759 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION), 3760 errmsg("lower bound cannot equal upper bound"))); 3761 result = 0; /* keep the compiler quiet */ 3762 } 3763 3764 PG_RETURN_INT32(result); 3765 } 3766 3767 /* ========== PRIVATE ROUTINES ========== */ 3768 3769 #ifndef HAVE_CBRT 3770 3771 static double 3772 cbrt(double x) 3773 { 3774 int isneg = (x < 0.0); 3775 double absx = fabs(x); 3776 double tmpres = pow(absx, (double) 1.0 / (double) 3.0); 3777 3778 /* 3779 * The result is somewhat inaccurate --- not really pow()'s fault, as the 3780 * exponent it's handed contains roundoff error. We can improve the 3781 * accuracy by doing one iteration of Newton's formula. Beware of zero 3782 * input however. 3783 */ 3784 if (tmpres > 0.0) 3785 tmpres -= (tmpres - absx / (tmpres * tmpres)) / (double) 3.0; 3786 3787 return isneg ? -tmpres : tmpres; 3788 } 3789 3790 #endif /* !HAVE_CBRT */ 3791