float.c source code [PostgreSQL/src/backend/utils/adt/float.c]

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

Browse the source code of PostgreSQL/src/backend/utils/adt/float.c