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