1
/*-------------------------------------------------------------------------
4
* An exact numeric data type for the Postgres database system
6
* Original coding 1998, Jan Wieck. Heavily revised 2003, Tom Lane.
8
* Many of the algorithmic ideas are borrowed from David M. Smith's "FM"
9
* multiple-precision math library, most recently published as Algorithm
10
* 786: Multiple-Precision Complex Arithmetic and Functions, ACM
11
* Transactions on Mathematical Software, Vol. 24, No. 4, December 1998,
14
* Copyright (c) 1998-2005, PostgreSQL Global Development Group
17
* $PostgreSQL: pgsql/src/backend/utils/adt/numeric.c,v 1.81 2005-01-01 05:43:07 momjian Exp $
19
*-------------------------------------------------------------------------
30
#include "catalog/pg_type.h"
31
#include "libpq/pqformat.h"
32
#include "utils/array.h"
33
#include "utils/builtins.h"
34
#include "utils/int8.h"
35
#include "utils/numeric.h"
38
* Uncomment the following to enable compilation of dump_numeric()
39
* and dump_var() and to get a dump of any result produced by make_result().
48
* Numeric values are represented in a base-NBASE floating point format.
49
* Each "digit" ranges from 0 to NBASE-1. The type NumericDigit is signed
50
* and wide enough to store a digit. We assume that NBASE*NBASE can fit in
51
* an int. Although the purely calculational routines could handle any even
52
* NBASE that's less than sqrt(INT_MAX), in practice we are only interested
53
* in NBASE a power of ten, so that I/O conversions and decimal rounding
54
* are easy. Also, it's actually more efficient if NBASE is rather less than
55
* sqrt(INT_MAX), so that there is "headroom" for mul_var and div_var to
56
* postpone processing carries.
63
#define DEC_DIGITS 1 /* decimal digits per NBASE digit */
64
#define MUL_GUARD_DIGITS 4 /* these are measured in NBASE digits */
65
#define DIV_GUARD_DIGITS 8
67
typedef signed char NumericDigit;
73
#define DEC_DIGITS 2 /* decimal digits per NBASE digit */
74
#define MUL_GUARD_DIGITS 3 /* these are measured in NBASE digits */
75
#define DIV_GUARD_DIGITS 6
77
typedef signed char NumericDigit;
82
#define HALF_NBASE 5000
83
#define DEC_DIGITS 4 /* decimal digits per NBASE digit */
84
#define MUL_GUARD_DIGITS 2 /* these are measured in NBASE digits */
85
#define DIV_GUARD_DIGITS 4
87
typedef int16 NumericDigit;
92
* The value represented by a NumericVar is determined by the sign, weight,
93
* ndigits, and digits[] array.
94
* Note: the first digit of a NumericVar's value is assumed to be multiplied
95
* by NBASE ** weight. Another way to say it is that there are weight+1
96
* digits before the decimal point. It is possible to have weight < 0.
98
* buf points at the physical start of the palloc'd digit buffer for the
99
* NumericVar. digits points at the first digit in actual use (the one
100
* with the specified weight). We normally leave an unused digit or two
101
* (preset to zeroes) between buf and digits, so that there is room to store
102
* a carry out of the top digit without special pushups. We just need to
103
* decrement digits (and increment weight) to make room for the carry digit.
104
* (There is no such extra space in a numeric value stored in the database,
105
* only in a NumericVar in memory.)
107
* If buf is NULL then the digit buffer isn't actually palloc'd and should
108
* not be freed --- see the constants below for an example.
110
* dscale, or display scale, is the nominal precision expressed as number
111
* of digits after the decimal point (it must always be >= 0 at present).
112
* dscale may be more than the number of physically stored fractional digits,
113
* implying that we have suppressed storage of significant trailing zeroes.
114
* It should never be less than the number of stored digits, since that would
115
* imply hiding digits that are present. NOTE that dscale is always expressed
116
* in *decimal* digits, and so it may correspond to a fractional number of
117
* base-NBASE digits --- divide by DEC_DIGITS to convert to NBASE digits.
119
* rscale, or result scale, is the target precision for a computation.
120
* Like dscale it is expressed as number of *decimal* digits after the decimal
121
* point, and is always >= 0 at present.
122
* Note that rscale is not stored in variables --- it's figured on-the-fly
123
* from the dscales of the inputs.
125
* NB: All the variable-level functions are written in a style that makes it
126
* possible to give one and the same variable as argument and destination.
127
* This is feasible because the digit buffer is separate from the variable.
130
typedef struct NumericVar
132
int ndigits; /* # of digits in digits[] - can be 0! */
133
int weight; /* weight of first digit */
134
int sign; /* NUMERIC_POS, NUMERIC_NEG, or
136
int dscale; /* display scale */
137
NumericDigit *buf; /* start of palloc'd space for digits[] */
138
NumericDigit *digits; /* base-NBASE digits */
143
* Some preinitialized constants
146
static NumericDigit const_zero_data[1] = {0};
147
static NumericVar const_zero =
148
{0, 0, NUMERIC_POS, 0, NULL, const_zero_data};
150
static NumericDigit const_one_data[1] = {1};
151
static NumericVar const_one =
152
{1, 0, NUMERIC_POS, 0, NULL, const_one_data};
154
static NumericDigit const_two_data[1] = {2};
155
static NumericVar const_two =
156
{1, 0, NUMERIC_POS, 0, NULL, const_two_data};
159
static NumericDigit const_zero_point_five_data[1] = {5000};
161
#elif DEC_DIGITS == 2
162
static NumericDigit const_zero_point_five_data[1] = {50};
164
#elif DEC_DIGITS == 1
165
static NumericDigit const_zero_point_five_data[1] = {5};
167
static NumericVar const_zero_point_five =
168
{1, -1, NUMERIC_POS, 1, NULL, const_zero_point_five_data};
171
static NumericDigit const_zero_point_nine_data[1] = {9000};
173
#elif DEC_DIGITS == 2
174
static NumericDigit const_zero_point_nine_data[1] = {90};
176
#elif DEC_DIGITS == 1
177
static NumericDigit const_zero_point_nine_data[1] = {9};
179
static NumericVar const_zero_point_nine =
180
{1, -1, NUMERIC_POS, 1, NULL, const_zero_point_nine_data};
183
static NumericDigit const_zero_point_01_data[1] = {100};
184
static NumericVar const_zero_point_01 =
185
{1, -1, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
187
#elif DEC_DIGITS == 2
188
static NumericDigit const_zero_point_01_data[1] = {1};
189
static NumericVar const_zero_point_01 =
190
{1, -1, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
192
#elif DEC_DIGITS == 1
193
static NumericDigit const_zero_point_01_data[1] = {1};
194
static NumericVar const_zero_point_01 =
195
{1, -2, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
199
static NumericDigit const_one_point_one_data[2] = {1, 1000};
201
#elif DEC_DIGITS == 2
202
static NumericDigit const_one_point_one_data[2] = {1, 10};
204
#elif DEC_DIGITS == 1
205
static NumericDigit const_one_point_one_data[2] = {1, 1};
207
static NumericVar const_one_point_one =
208
{2, 0, NUMERIC_POS, 1, NULL, const_one_point_one_data};
210
static NumericVar const_nan =
211
{0, 0, NUMERIC_NAN, 0, NULL, NULL};
214
static const int round_powers[4] = {0, 1000, 100, 10};
224
static void dump_numeric(const char *str, Numeric num);
225
static void dump_var(const char *str, NumericVar *var);
228
#define dump_numeric(s,n)
229
#define dump_var(s,v)
232
#define digitbuf_alloc(ndigits) \
233
((NumericDigit *) palloc((ndigits) * sizeof(NumericDigit)))
234
#define digitbuf_free(buf) \
240
#define init_var(v) MemSetAligned(v, 0, sizeof(NumericVar))
242
static void alloc_var(NumericVar *var, int ndigits);
243
static void free_var(NumericVar *var);
244
static void zero_var(NumericVar *var);
246
static void set_var_from_str(const char *str, NumericVar *dest);
247
static void set_var_from_num(Numeric value, NumericVar *dest);
248
static void set_var_from_var(NumericVar *value, NumericVar *dest);
249
static char *get_str_from_var(NumericVar *var, int dscale);
251
static Numeric make_result(NumericVar *var);
253
static void apply_typmod(NumericVar *var, int32 typmod);
255
static int32 numericvar_to_int4(NumericVar *var);
256
static bool numericvar_to_int8(NumericVar *var, int64 *result);
257
static void int8_to_numericvar(int64 val, NumericVar *var);
258
static double numeric_to_double_no_overflow(Numeric num);
259
static double numericvar_to_double_no_overflow(NumericVar *var);
261
static int cmp_numerics(Numeric num1, Numeric num2);
262
static int cmp_var(NumericVar *var1, NumericVar *var2);
263
static void add_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
264
static void sub_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
265
static void mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
267
static void div_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
269
static int select_div_scale(NumericVar *var1, NumericVar *var2);
270
static void mod_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
271
static void ceil_var(NumericVar *var, NumericVar *result);
272
static void floor_var(NumericVar *var, NumericVar *result);
274
static void sqrt_var(NumericVar *arg, NumericVar *result, int rscale);
275
static void exp_var(NumericVar *arg, NumericVar *result, int rscale);
276
static void exp_var_internal(NumericVar *arg, NumericVar *result, int rscale);
277
static void ln_var(NumericVar *arg, NumericVar *result, int rscale);
278
static void log_var(NumericVar *base, NumericVar *num, NumericVar *result);
279
static void power_var(NumericVar *base, NumericVar *exp, NumericVar *result);
280
static void power_var_int(NumericVar *base, int exp, NumericVar *result,
283
static int cmp_abs(NumericVar *var1, NumericVar *var2);
284
static void add_abs(NumericVar *var1, NumericVar *var2, NumericVar *result);
285
static void sub_abs(NumericVar *var1, NumericVar *var2, NumericVar *result);
286
static void round_var(NumericVar *var, int rscale);
287
static void trunc_var(NumericVar *var, int rscale);
288
static void strip_var(NumericVar *var);
289
static void compute_bucket(Numeric operand, Numeric bound1, Numeric bound2,
290
NumericVar *count_var, NumericVar *result_var);
293
/* ----------------------------------------------------------------------
295
* Input-, output- and rounding-functions
297
* ----------------------------------------------------------------------
304
* Input function for numeric data type
307
numeric_in(PG_FUNCTION_ARGS)
309
char *str = PG_GETARG_CSTRING(0);
312
Oid typelem = PG_GETARG_OID(1);
314
int32 typmod = PG_GETARG_INT32(2);
321
if (pg_strcasecmp(str, "NaN") == 0)
322
PG_RETURN_NUMERIC(make_result(&const_nan));
325
* Use set_var_from_str() to parse the input string and return it in
326
* the packed DB storage format
329
set_var_from_str(str, &value);
331
apply_typmod(&value, typmod);
333
res = make_result(&value);
336
PG_RETURN_NUMERIC(res);
343
* Output function for numeric data type
346
numeric_out(PG_FUNCTION_ARGS)
348
Numeric num = PG_GETARG_NUMERIC(0);
355
if (NUMERIC_IS_NAN(num))
356
PG_RETURN_CSTRING(pstrdup("NaN"));
359
* Get the number in the variable format.
361
* Even if we didn't need to change format, we'd still need to copy the
362
* value to have a modifiable copy for rounding. set_var_from_num()
363
* also guarantees there is extra digit space in case we produce a
364
* carry out from rounding.
367
set_var_from_num(num, &x);
369
str = get_str_from_var(&x, x.dscale);
373
PG_RETURN_CSTRING(str);
377
* numeric_recv - converts external binary format to numeric
379
* External format is a sequence of int16's:
380
* ndigits, weight, sign, dscale, NumericDigits.
383
numeric_recv(PG_FUNCTION_ARGS)
385
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
393
len = (uint16) pq_getmsgint(buf, sizeof(uint16));
394
if (len < 0 || len > NUMERIC_MAX_PRECISION + NUMERIC_MAX_RESULT_SCALE)
396
(errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
397
errmsg("invalid length in external \"numeric\" value")));
399
alloc_var(&value, len);
401
value.weight = (int16) pq_getmsgint(buf, sizeof(int16));
402
value.sign = (uint16) pq_getmsgint(buf, sizeof(uint16));
403
if (!(value.sign == NUMERIC_POS ||
404
value.sign == NUMERIC_NEG ||
405
value.sign == NUMERIC_NAN))
407
(errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
408
errmsg("invalid sign in external \"numeric\" value")));
410
value.dscale = (uint16) pq_getmsgint(buf, sizeof(uint16));
411
for (i = 0; i < len; i++)
413
NumericDigit d = pq_getmsgint(buf, sizeof(NumericDigit));
415
if (d < 0 || d >= NBASE)
417
(errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
418
errmsg("invalid digit in external \"numeric\" value")));
422
res = make_result(&value);
425
PG_RETURN_NUMERIC(res);
429
* numeric_send - converts numeric to binary format
432
numeric_send(PG_FUNCTION_ARGS)
434
Numeric num = PG_GETARG_NUMERIC(0);
440
set_var_from_num(num, &x);
442
pq_begintypsend(&buf);
444
pq_sendint(&buf, x.ndigits, sizeof(int16));
445
pq_sendint(&buf, x.weight, sizeof(int16));
446
pq_sendint(&buf, x.sign, sizeof(int16));
447
pq_sendint(&buf, x.dscale, sizeof(int16));
448
for (i = 0; i < x.ndigits; i++)
449
pq_sendint(&buf, x.digits[i], sizeof(NumericDigit));
453
PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
460
* This is a special function called by the Postgres database system
461
* before a value is stored in a tuple's attribute. The precision and
462
* scale of the attribute have to be applied on the value.
465
numeric (PG_FUNCTION_ARGS)
467
Numeric num = PG_GETARG_NUMERIC(0);
468
int32 typmod = PG_GETARG_INT32(1);
480
if (NUMERIC_IS_NAN(num))
481
PG_RETURN_NUMERIC(make_result(&const_nan));
484
* If the value isn't a valid type modifier, simply return a copy of
487
if (typmod < (int32) (VARHDRSZ))
489
new = (Numeric) palloc(num->varlen);
490
memcpy(new, num, num->varlen);
491
PG_RETURN_NUMERIC(new);
495
* Get the precision and scale out of the typmod value
497
tmp_typmod = typmod - VARHDRSZ;
498
precision = (tmp_typmod >> 16) & 0xffff;
499
scale = tmp_typmod & 0xffff;
500
maxdigits = precision - scale;
503
* If the number is certainly in bounds and due to the target scale no
504
* rounding could be necessary, just make a copy of the input and
505
* modify its scale fields. (Note we assume the existing dscale is
508
ddigits = (num->n_weight + 1) * DEC_DIGITS;
509
if (ddigits <= maxdigits && scale >= NUMERIC_DSCALE(num))
511
new = (Numeric) palloc(num->varlen);
512
memcpy(new, num, num->varlen);
513
new->n_sign_dscale = NUMERIC_SIGN(new) |
514
((uint16) scale & NUMERIC_DSCALE_MASK);
515
PG_RETURN_NUMERIC(new);
519
* We really need to fiddle with things - unpack the number into a
520
* variable and let apply_typmod() do it.
524
set_var_from_num(num, &var);
525
apply_typmod(&var, typmod);
526
new = make_result(&var);
530
PG_RETURN_NUMERIC(new);
534
/* ----------------------------------------------------------------------
536
* Sign manipulation, rounding and the like
538
* ----------------------------------------------------------------------
542
numeric_abs(PG_FUNCTION_ARGS)
544
Numeric num = PG_GETARG_NUMERIC(0);
550
if (NUMERIC_IS_NAN(num))
551
PG_RETURN_NUMERIC(make_result(&const_nan));
554
* Do it the easy way directly on the packed format
556
res = (Numeric) palloc(num->varlen);
557
memcpy(res, num, num->varlen);
559
res->n_sign_dscale = NUMERIC_POS | NUMERIC_DSCALE(num);
561
PG_RETURN_NUMERIC(res);
566
numeric_uminus(PG_FUNCTION_ARGS)
568
Numeric num = PG_GETARG_NUMERIC(0);
574
if (NUMERIC_IS_NAN(num))
575
PG_RETURN_NUMERIC(make_result(&const_nan));
578
* Do it the easy way directly on the packed format
580
res = (Numeric) palloc(num->varlen);
581
memcpy(res, num, num->varlen);
584
* The packed format is known to be totally zero digit trimmed always.
585
* So we can identify a ZERO by the fact that there are no digits at
586
* all. Do nothing to a zero.
588
if (num->varlen != NUMERIC_HDRSZ)
590
/* Else, flip the sign */
591
if (NUMERIC_SIGN(num) == NUMERIC_POS)
592
res->n_sign_dscale = NUMERIC_NEG | NUMERIC_DSCALE(num);
594
res->n_sign_dscale = NUMERIC_POS | NUMERIC_DSCALE(num);
597
PG_RETURN_NUMERIC(res);
602
numeric_uplus(PG_FUNCTION_ARGS)
604
Numeric num = PG_GETARG_NUMERIC(0);
607
res = (Numeric) palloc(num->varlen);
608
memcpy(res, num, num->varlen);
610
PG_RETURN_NUMERIC(res);
616
* returns -1 if the argument is less than 0, 0 if the argument is equal
617
* to 0, and 1 if the argument is greater than zero.
620
numeric_sign(PG_FUNCTION_ARGS)
622
Numeric num = PG_GETARG_NUMERIC(0);
629
if (NUMERIC_IS_NAN(num))
630
PG_RETURN_NUMERIC(make_result(&const_nan));
635
* The packed format is known to be totally zero digit trimmed always.
636
* So we can identify a ZERO by the fact that there are no digits at
639
if (num->varlen == NUMERIC_HDRSZ)
640
set_var_from_var(&const_zero, &result);
644
* And if there are some, we return a copy of ONE with the sign of
647
set_var_from_var(&const_one, &result);
648
result.sign = NUMERIC_SIGN(num);
651
res = make_result(&result);
654
PG_RETURN_NUMERIC(res);
661
* Round a value to have 'scale' digits after the decimal point.
662
* We allow negative 'scale', implying rounding before the decimal
663
* point --- Oracle interprets rounding that way.
666
numeric_round(PG_FUNCTION_ARGS)
668
Numeric num = PG_GETARG_NUMERIC(0);
669
int32 scale = PG_GETARG_INT32(1);
676
if (NUMERIC_IS_NAN(num))
677
PG_RETURN_NUMERIC(make_result(&const_nan));
680
* Limit the scale value to avoid possible overflow in calculations
682
scale = Max(scale, -NUMERIC_MAX_RESULT_SCALE);
683
scale = Min(scale, NUMERIC_MAX_RESULT_SCALE);
686
* Unpack the argument and round it at the proper digit position
689
set_var_from_num(num, &arg);
691
round_var(&arg, scale);
693
/* We don't allow negative output dscale */
698
* Return the rounded result
700
res = make_result(&arg);
703
PG_RETURN_NUMERIC(res);
710
* Truncate a value to have 'scale' digits after the decimal point.
711
* We allow negative 'scale', implying a truncation before the decimal
712
* point --- Oracle interprets truncation that way.
715
numeric_trunc(PG_FUNCTION_ARGS)
717
Numeric num = PG_GETARG_NUMERIC(0);
718
int32 scale = PG_GETARG_INT32(1);
725
if (NUMERIC_IS_NAN(num))
726
PG_RETURN_NUMERIC(make_result(&const_nan));
729
* Limit the scale value to avoid possible overflow in calculations
731
scale = Max(scale, -NUMERIC_MAX_RESULT_SCALE);
732
scale = Min(scale, NUMERIC_MAX_RESULT_SCALE);
735
* Unpack the argument and truncate it at the proper digit position
738
set_var_from_num(num, &arg);
740
trunc_var(&arg, scale);
742
/* We don't allow negative output dscale */
747
* Return the truncated result
749
res = make_result(&arg);
752
PG_RETURN_NUMERIC(res);
759
* Return the smallest integer greater than or equal to the argument
762
numeric_ceil(PG_FUNCTION_ARGS)
764
Numeric num = PG_GETARG_NUMERIC(0);
768
if (NUMERIC_IS_NAN(num))
769
PG_RETURN_NUMERIC(make_result(&const_nan));
773
set_var_from_num(num, &result);
774
ceil_var(&result, &result);
776
res = make_result(&result);
779
PG_RETURN_NUMERIC(res);
786
* Return the largest integer equal to or less than the argument
789
numeric_floor(PG_FUNCTION_ARGS)
791
Numeric num = PG_GETARG_NUMERIC(0);
795
if (NUMERIC_IS_NAN(num))
796
PG_RETURN_NUMERIC(make_result(&const_nan));
800
set_var_from_num(num, &result);
801
floor_var(&result, &result);
803
res = make_result(&result);
806
PG_RETURN_NUMERIC(res);
810
* width_bucket_numeric() -
812
* 'bound1' and 'bound2' are the lower and upper bounds of the
813
* histogram's range, respectively. 'count' is the number of buckets
814
* in the histogram. width_bucket() returns an integer indicating the
815
* bucket number that 'operand' belongs in for an equiwidth histogram
816
* with the specified characteristics. An operand smaller than the
817
* lower bound is assigned to bucket 0. An operand greater than the
818
* upper bound is assigned to an additional bucket (with number
822
width_bucket_numeric(PG_FUNCTION_ARGS)
824
Numeric operand = PG_GETARG_NUMERIC(0);
825
Numeric bound1 = PG_GETARG_NUMERIC(1);
826
Numeric bound2 = PG_GETARG_NUMERIC(2);
827
int32 count = PG_GETARG_INT32(3);
828
NumericVar count_var;
829
NumericVar result_var;
834
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
835
errmsg("count must be greater than zero")));
837
init_var(&result_var);
838
init_var(&count_var);
840
/* Convert 'count' to a numeric, for ease of use later */
841
int8_to_numericvar((int64) count, &count_var);
843
switch (cmp_numerics(bound1, bound2))
847
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
848
errmsg("lower bound cannot equal upper bound")));
850
/* bound1 < bound2 */
852
if (cmp_numerics(operand, bound1) < 0)
853
set_var_from_var(&const_zero, &result_var);
854
else if (cmp_numerics(operand, bound2) >= 0)
855
add_var(&count_var, &const_one, &result_var);
857
compute_bucket(operand, bound1, bound2,
858
&count_var, &result_var);
861
/* bound1 > bound2 */
863
if (cmp_numerics(operand, bound1) > 0)
864
set_var_from_var(&const_zero, &result_var);
865
else if (cmp_numerics(operand, bound2) <= 0)
866
add_var(&count_var, &const_one, &result_var);
868
compute_bucket(operand, bound1, bound2,
869
&count_var, &result_var);
873
result = numericvar_to_int4(&result_var);
875
free_var(&count_var);
876
free_var(&result_var);
878
PG_RETURN_INT32(result);
884
* If 'operand' is not outside the bucket range, determine the correct
885
* bucket for it to go. The calculations performed by this function
886
* are derived directly from the SQL2003 spec.
889
compute_bucket(Numeric operand, Numeric bound1, Numeric bound2,
890
NumericVar *count_var, NumericVar *result_var)
892
NumericVar bound1_var;
893
NumericVar bound2_var;
894
NumericVar operand_var;
896
init_var(&bound1_var);
897
init_var(&bound2_var);
898
init_var(&operand_var);
900
set_var_from_num(bound1, &bound1_var);
901
set_var_from_num(bound2, &bound2_var);
902
set_var_from_num(operand, &operand_var);
904
if (cmp_var(&bound1_var, &bound2_var) < 0)
906
sub_var(&operand_var, &bound1_var, &operand_var);
907
sub_var(&bound2_var, &bound1_var, &bound2_var);
908
div_var(&operand_var, &bound2_var, result_var,
909
select_div_scale(&operand_var, &bound2_var));
913
sub_var(&bound1_var, &operand_var, &operand_var);
914
sub_var(&bound1_var, &bound2_var, &bound1_var);
915
div_var(&operand_var, &bound1_var, result_var,
916
select_div_scale(&operand_var, &bound1_var));
919
mul_var(result_var, count_var, result_var,
920
result_var->dscale + count_var->dscale);
921
add_var(result_var, &const_one, result_var);
922
floor_var(result_var, result_var);
924
free_var(&bound1_var);
925
free_var(&bound2_var);
926
free_var(&operand_var);
929
/* ----------------------------------------------------------------------
931
* Comparison functions
933
* Note: btree indexes need these routines not to leak memory; therefore,
934
* be careful to free working copies of toasted datums. Most places don't
935
* need to be so careful.
936
* ----------------------------------------------------------------------
941
numeric_cmp(PG_FUNCTION_ARGS)
943
Numeric num1 = PG_GETARG_NUMERIC(0);
944
Numeric num2 = PG_GETARG_NUMERIC(1);
947
result = cmp_numerics(num1, num2);
949
PG_FREE_IF_COPY(num1, 0);
950
PG_FREE_IF_COPY(num2, 1);
952
PG_RETURN_INT32(result);
957
numeric_eq(PG_FUNCTION_ARGS)
959
Numeric num1 = PG_GETARG_NUMERIC(0);
960
Numeric num2 = PG_GETARG_NUMERIC(1);
963
result = cmp_numerics(num1, num2) == 0;
965
PG_FREE_IF_COPY(num1, 0);
966
PG_FREE_IF_COPY(num2, 1);
968
PG_RETURN_BOOL(result);
972
numeric_ne(PG_FUNCTION_ARGS)
974
Numeric num1 = PG_GETARG_NUMERIC(0);
975
Numeric num2 = PG_GETARG_NUMERIC(1);
978
result = cmp_numerics(num1, num2) != 0;
980
PG_FREE_IF_COPY(num1, 0);
981
PG_FREE_IF_COPY(num2, 1);
983
PG_RETURN_BOOL(result);
987
numeric_gt(PG_FUNCTION_ARGS)
989
Numeric num1 = PG_GETARG_NUMERIC(0);
990
Numeric num2 = PG_GETARG_NUMERIC(1);
993
result = cmp_numerics(num1, num2) > 0;
995
PG_FREE_IF_COPY(num1, 0);
996
PG_FREE_IF_COPY(num2, 1);
998
PG_RETURN_BOOL(result);
1002
numeric_ge(PG_FUNCTION_ARGS)
1004
Numeric num1 = PG_GETARG_NUMERIC(0);
1005
Numeric num2 = PG_GETARG_NUMERIC(1);
1008
result = cmp_numerics(num1, num2) >= 0;
1010
PG_FREE_IF_COPY(num1, 0);
1011
PG_FREE_IF_COPY(num2, 1);
1013
PG_RETURN_BOOL(result);
1017
numeric_lt(PG_FUNCTION_ARGS)
1019
Numeric num1 = PG_GETARG_NUMERIC(0);
1020
Numeric num2 = PG_GETARG_NUMERIC(1);
1023
result = cmp_numerics(num1, num2) < 0;
1025
PG_FREE_IF_COPY(num1, 0);
1026
PG_FREE_IF_COPY(num2, 1);
1028
PG_RETURN_BOOL(result);
1032
numeric_le(PG_FUNCTION_ARGS)
1034
Numeric num1 = PG_GETARG_NUMERIC(0);
1035
Numeric num2 = PG_GETARG_NUMERIC(1);
1038
result = cmp_numerics(num1, num2) <= 0;
1040
PG_FREE_IF_COPY(num1, 0);
1041
PG_FREE_IF_COPY(num2, 1);
1043
PG_RETURN_BOOL(result);
1047
cmp_numerics(Numeric num1, Numeric num2)
1052
* We consider all NANs to be equal and larger than any non-NAN. This
1053
* is somewhat arbitrary; the important thing is to have a consistent
1056
if (NUMERIC_IS_NAN(num1))
1058
if (NUMERIC_IS_NAN(num2))
1059
result = 0; /* NAN = NAN */
1061
result = 1; /* NAN > non-NAN */
1063
else if (NUMERIC_IS_NAN(num2))
1065
result = -1; /* non-NAN < NAN */
1075
set_var_from_num(num1, &arg1);
1076
set_var_from_num(num2, &arg2);
1078
result = cmp_var(&arg1, &arg2);
1088
/* ----------------------------------------------------------------------
1090
* Basic arithmetic functions
1092
* ----------------------------------------------------------------------
1102
numeric_add(PG_FUNCTION_ARGS)
1104
Numeric num1 = PG_GETARG_NUMERIC(0);
1105
Numeric num2 = PG_GETARG_NUMERIC(1);
1114
if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1115
PG_RETURN_NUMERIC(make_result(&const_nan));
1118
* Unpack the values, let add_var() compute the result and return it.
1124
set_var_from_num(num1, &arg1);
1125
set_var_from_num(num2, &arg2);
1127
add_var(&arg1, &arg2, &result);
1129
res = make_result(&result);
1135
PG_RETURN_NUMERIC(res);
1142
* Subtract one numeric from another
1145
numeric_sub(PG_FUNCTION_ARGS)
1147
Numeric num1 = PG_GETARG_NUMERIC(0);
1148
Numeric num2 = PG_GETARG_NUMERIC(1);
1157
if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1158
PG_RETURN_NUMERIC(make_result(&const_nan));
1161
* Unpack the values, let sub_var() compute the result and return it.
1167
set_var_from_num(num1, &arg1);
1168
set_var_from_num(num2, &arg2);
1170
sub_var(&arg1, &arg2, &result);
1172
res = make_result(&result);
1178
PG_RETURN_NUMERIC(res);
1185
* Calculate the product of two numerics
1188
numeric_mul(PG_FUNCTION_ARGS)
1190
Numeric num1 = PG_GETARG_NUMERIC(0);
1191
Numeric num2 = PG_GETARG_NUMERIC(1);
1200
if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1201
PG_RETURN_NUMERIC(make_result(&const_nan));
1204
* Unpack the values, let mul_var() compute the result and return it.
1205
* Unlike add_var() and sub_var(), mul_var() will round its result. In
1206
* the case of numeric_mul(), which is invoked for the * operator on
1207
* numerics, we request exact representation for the product (rscale =
1208
* sum(dscale of arg1, dscale of arg2)).
1214
set_var_from_num(num1, &arg1);
1215
set_var_from_num(num2, &arg2);
1217
mul_var(&arg1, &arg2, &result, arg1.dscale + arg2.dscale);
1219
res = make_result(&result);
1225
PG_RETURN_NUMERIC(res);
1232
* Divide one numeric into another
1235
numeric_div(PG_FUNCTION_ARGS)
1237
Numeric num1 = PG_GETARG_NUMERIC(0);
1238
Numeric num2 = PG_GETARG_NUMERIC(1);
1248
if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1249
PG_RETURN_NUMERIC(make_result(&const_nan));
1252
* Unpack the arguments
1258
set_var_from_num(num1, &arg1);
1259
set_var_from_num(num2, &arg2);
1262
* Select scale for division result
1264
rscale = select_div_scale(&arg1, &arg2);
1267
* Do the divide and return the result
1269
div_var(&arg1, &arg2, &result, rscale);
1271
res = make_result(&result);
1277
PG_RETURN_NUMERIC(res);
1284
* Calculate the modulo of two numerics
1287
numeric_mod(PG_FUNCTION_ARGS)
1289
Numeric num1 = PG_GETARG_NUMERIC(0);
1290
Numeric num2 = PG_GETARG_NUMERIC(1);
1296
if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1297
PG_RETURN_NUMERIC(make_result(&const_nan));
1303
set_var_from_num(num1, &arg1);
1304
set_var_from_num(num2, &arg2);
1306
mod_var(&arg1, &arg2, &result);
1308
res = make_result(&result);
1314
PG_RETURN_NUMERIC(res);
1321
* Increment a number by one
1324
numeric_inc(PG_FUNCTION_ARGS)
1326
Numeric num = PG_GETARG_NUMERIC(0);
1333
if (NUMERIC_IS_NAN(num))
1334
PG_RETURN_NUMERIC(make_result(&const_nan));
1337
* Compute the result and return it
1341
set_var_from_num(num, &arg);
1343
add_var(&arg, &const_one, &arg);
1345
res = make_result(&arg);
1349
PG_RETURN_NUMERIC(res);
1354
* numeric_smaller() -
1356
* Return the smaller of two numbers
1359
numeric_smaller(PG_FUNCTION_ARGS)
1361
Numeric num1 = PG_GETARG_NUMERIC(0);
1362
Numeric num2 = PG_GETARG_NUMERIC(1);
1365
* Use cmp_numerics so that this will agree with the comparison
1366
* operators, particularly as regards comparisons involving NaN.
1368
if (cmp_numerics(num1, num2) < 0)
1369
PG_RETURN_NUMERIC(num1);
1371
PG_RETURN_NUMERIC(num2);
1376
* numeric_larger() -
1378
* Return the larger of two numbers
1381
numeric_larger(PG_FUNCTION_ARGS)
1383
Numeric num1 = PG_GETARG_NUMERIC(0);
1384
Numeric num2 = PG_GETARG_NUMERIC(1);
1387
* Use cmp_numerics so that this will agree with the comparison
1388
* operators, particularly as regards comparisons involving NaN.
1390
if (cmp_numerics(num1, num2) > 0)
1391
PG_RETURN_NUMERIC(num1);
1393
PG_RETURN_NUMERIC(num2);
1397
/* ----------------------------------------------------------------------
1399
* Advanced math functions
1401
* ----------------------------------------------------------------------
1410
numeric_fac(PG_FUNCTION_ARGS)
1412
int64 num = PG_GETARG_INT64(0);
1419
res = make_result(&const_one);
1420
PG_RETURN_NUMERIC(res);
1426
int8_to_numericvar(num, &result);
1428
for (num = num - 1; num > 1; num--)
1430
int8_to_numericvar(num, &fact);
1432
mul_var(&result, &fact, &result, 0);
1435
res = make_result(&result);
1440
PG_RETURN_NUMERIC(res);
1447
* Compute the square root of a numeric.
1450
numeric_sqrt(PG_FUNCTION_ARGS)
1452
Numeric num = PG_GETARG_NUMERIC(0);
1462
if (NUMERIC_IS_NAN(num))
1463
PG_RETURN_NUMERIC(make_result(&const_nan));
1466
* Unpack the argument and determine the result scale. We choose a
1467
* scale to give at least NUMERIC_MIN_SIG_DIGITS significant digits;
1468
* but in any case not less than the input's dscale.
1473
set_var_from_num(num, &arg);
1475
/* Assume the input was normalized, so arg.weight is accurate */
1476
sweight = (arg.weight + 1) * DEC_DIGITS / 2 - 1;
1478
rscale = NUMERIC_MIN_SIG_DIGITS - sweight;
1479
rscale = Max(rscale, arg.dscale);
1480
rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
1481
rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
1484
* Let sqrt_var() do the calculation and return the result.
1486
sqrt_var(&arg, &result, rscale);
1488
res = make_result(&result);
1493
PG_RETURN_NUMERIC(res);
1500
* Raise e to the power of x
1503
numeric_exp(PG_FUNCTION_ARGS)
1505
Numeric num = PG_GETARG_NUMERIC(0);
1515
if (NUMERIC_IS_NAN(num))
1516
PG_RETURN_NUMERIC(make_result(&const_nan));
1519
* Unpack the argument and determine the result scale. We choose a
1520
* scale to give at least NUMERIC_MIN_SIG_DIGITS significant digits;
1521
* but in any case not less than the input's dscale.
1526
set_var_from_num(num, &arg);
1528
/* convert input to float8, ignoring overflow */
1529
val = numericvar_to_double_no_overflow(&arg);
1532
* log10(result) = num * log10(e), so this is approximately the
1533
* decimal weight of the result:
1535
val *= 0.434294481903252;
1537
/* limit to something that won't cause integer overflow */
1538
val = Max(val, -NUMERIC_MAX_RESULT_SCALE);
1539
val = Min(val, NUMERIC_MAX_RESULT_SCALE);
1541
rscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
1542
rscale = Max(rscale, arg.dscale);
1543
rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
1544
rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
1547
* Let exp_var() do the calculation and return the result.
1549
exp_var(&arg, &result, rscale);
1551
res = make_result(&result);
1556
PG_RETURN_NUMERIC(res);
1563
* Compute the natural logarithm of x
1566
numeric_ln(PG_FUNCTION_ARGS)
1568
Numeric num = PG_GETARG_NUMERIC(0);
1578
if (NUMERIC_IS_NAN(num))
1579
PG_RETURN_NUMERIC(make_result(&const_nan));
1584
set_var_from_num(num, &arg);
1586
/* Approx decimal digits before decimal point */
1587
dec_digits = (arg.weight + 1) * DEC_DIGITS;
1590
rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(dec_digits - 1);
1591
else if (dec_digits < 1)
1592
rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(1 - dec_digits);
1594
rscale = NUMERIC_MIN_SIG_DIGITS;
1596
rscale = Max(rscale, arg.dscale);
1597
rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
1598
rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
1600
ln_var(&arg, &result, rscale);
1602
res = make_result(&result);
1607
PG_RETURN_NUMERIC(res);
1614
* Compute the logarithm of x in a given base
1617
numeric_log(PG_FUNCTION_ARGS)
1619
Numeric num1 = PG_GETARG_NUMERIC(0);
1620
Numeric num2 = PG_GETARG_NUMERIC(1);
1629
if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1630
PG_RETURN_NUMERIC(make_result(&const_nan));
1639
set_var_from_num(num1, &arg1);
1640
set_var_from_num(num2, &arg2);
1643
* Call log_var() to compute and return the result; note it handles
1644
* scale selection itself.
1646
log_var(&arg1, &arg2, &result);
1648
res = make_result(&result);
1654
PG_RETURN_NUMERIC(res);
1661
* Raise b to the power of x
1664
numeric_power(PG_FUNCTION_ARGS)
1666
Numeric num1 = PG_GETARG_NUMERIC(0);
1667
Numeric num2 = PG_GETARG_NUMERIC(1);
1671
NumericVar arg2_trunc;
1677
if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1678
PG_RETURN_NUMERIC(make_result(&const_nan));
1685
init_var(&arg2_trunc);
1688
set_var_from_num(num1, &arg1);
1689
set_var_from_num(num2, &arg2);
1690
set_var_from_var(&arg2, &arg2_trunc);
1692
trunc_var(&arg2_trunc, 0);
1695
* Return special SQLSTATE error codes for a few conditions mandated
1698
if ((cmp_var(&arg1, &const_zero) == 0 &&
1699
cmp_var(&arg2, &const_zero) < 0) ||
1700
(cmp_var(&arg1, &const_zero) < 0 &&
1701
cmp_var(&arg2, &arg2_trunc) != 0))
1703
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
1704
errmsg("invalid argument for power function")));
1707
* Call power_var() to compute and return the result; note it handles
1708
* scale selection itself.
1710
power_var(&arg1, &arg2, &result);
1712
res = make_result(&result);
1716
free_var(&arg2_trunc);
1719
PG_RETURN_NUMERIC(res);
1723
/* ----------------------------------------------------------------------
1725
* Type conversion functions
1727
* ----------------------------------------------------------------------
1732
int4_numeric(PG_FUNCTION_ARGS)
1734
int32 val = PG_GETARG_INT32(0);
1740
int8_to_numericvar((int64) val, &result);
1742
res = make_result(&result);
1746
PG_RETURN_NUMERIC(res);
1751
numeric_int4(PG_FUNCTION_ARGS)
1753
Numeric num = PG_GETARG_NUMERIC(0);
1757
/* XXX would it be better to return NULL? */
1758
if (NUMERIC_IS_NAN(num))
1760
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1761
errmsg("cannot convert NaN to integer")));
1763
/* Convert to variable format, then convert to int4 */
1765
set_var_from_num(num, &x);
1766
result = numericvar_to_int4(&x);
1768
PG_RETURN_INT32(result);
1772
* Given a NumericVar, convert it to an int32. If the NumericVar
1773
* exceeds the range of an int32, raise the appropriate error via
1774
* ereport(). The input NumericVar is *not* free'd.
1777
numericvar_to_int4(NumericVar *var)
1782
if (!numericvar_to_int8(var, &val))
1784
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1785
errmsg("integer out of range")));
1787
/* Down-convert to int4 */
1788
result = (int32) val;
1790
/* Test for overflow by reverse-conversion. */
1791
if ((int64) result != val)
1793
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1794
errmsg("integer out of range")));
1800
int8_numeric(PG_FUNCTION_ARGS)
1802
int64 val = PG_GETARG_INT64(0);
1808
int8_to_numericvar(val, &result);
1810
res = make_result(&result);
1814
PG_RETURN_NUMERIC(res);
1819
numeric_int8(PG_FUNCTION_ARGS)
1821
Numeric num = PG_GETARG_NUMERIC(0);
1825
/* XXX would it be better to return NULL? */
1826
if (NUMERIC_IS_NAN(num))
1828
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1829
errmsg("cannot convert NaN to bigint")));
1831
/* Convert to variable format and thence to int8 */
1833
set_var_from_num(num, &x);
1835
if (!numericvar_to_int8(&x, &result))
1837
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1838
errmsg("bigint out of range")));
1842
PG_RETURN_INT64(result);
1847
int2_numeric(PG_FUNCTION_ARGS)
1849
int16 val = PG_GETARG_INT16(0);
1855
int8_to_numericvar((int64) val, &result);
1857
res = make_result(&result);
1861
PG_RETURN_NUMERIC(res);
1866
numeric_int2(PG_FUNCTION_ARGS)
1868
Numeric num = PG_GETARG_NUMERIC(0);
1873
/* XXX would it be better to return NULL? */
1874
if (NUMERIC_IS_NAN(num))
1876
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1877
errmsg("cannot convert NaN to smallint")));
1879
/* Convert to variable format and thence to int8 */
1881
set_var_from_num(num, &x);
1883
if (!numericvar_to_int8(&x, &val))
1885
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1886
errmsg("smallint out of range")));
1890
/* Down-convert to int2 */
1891
result = (int16) val;
1893
/* Test for overflow by reverse-conversion. */
1894
if ((int64) result != val)
1896
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1897
errmsg("smallint out of range")));
1899
PG_RETURN_INT16(result);
1904
float8_numeric(PG_FUNCTION_ARGS)
1906
float8 val = PG_GETARG_FLOAT8(0);
1909
char buf[DBL_DIG + 100];
1912
PG_RETURN_NUMERIC(make_result(&const_nan));
1914
sprintf(buf, "%.*g", DBL_DIG, val);
1918
set_var_from_str(buf, &result);
1919
res = make_result(&result);
1923
PG_RETURN_NUMERIC(res);
1928
numeric_float8(PG_FUNCTION_ARGS)
1930
Numeric num = PG_GETARG_NUMERIC(0);
1934
if (NUMERIC_IS_NAN(num))
1935
PG_RETURN_FLOAT8(get_float8_nan());
1937
tmp = DatumGetCString(DirectFunctionCall1(numeric_out,
1938
NumericGetDatum(num)));
1940
result = DirectFunctionCall1(float8in, CStringGetDatum(tmp));
1944
PG_RETURN_DATUM(result);
1948
/* Convert numeric to float8; if out of range, return +/- HUGE_VAL */
1950
numeric_float8_no_overflow(PG_FUNCTION_ARGS)
1952
Numeric num = PG_GETARG_NUMERIC(0);
1955
if (NUMERIC_IS_NAN(num))
1956
PG_RETURN_FLOAT8(get_float8_nan());
1958
val = numeric_to_double_no_overflow(num);
1960
PG_RETURN_FLOAT8(val);
1964
float4_numeric(PG_FUNCTION_ARGS)
1966
float4 val = PG_GETARG_FLOAT4(0);
1969
char buf[FLT_DIG + 100];
1972
PG_RETURN_NUMERIC(make_result(&const_nan));
1974
sprintf(buf, "%.*g", FLT_DIG, val);
1978
set_var_from_str(buf, &result);
1979
res = make_result(&result);
1983
PG_RETURN_NUMERIC(res);
1988
numeric_float4(PG_FUNCTION_ARGS)
1990
Numeric num = PG_GETARG_NUMERIC(0);
1994
if (NUMERIC_IS_NAN(num))
1995
PG_RETURN_FLOAT4(get_float4_nan());
1997
tmp = DatumGetCString(DirectFunctionCall1(numeric_out,
1998
NumericGetDatum(num)));
2000
result = DirectFunctionCall1(float4in, CStringGetDatum(tmp));
2004
PG_RETURN_DATUM(result);
2009
text_numeric(PG_FUNCTION_ARGS)
2011
text *str = PG_GETARG_TEXT_P(0);
2016
len = (VARSIZE(str) - VARHDRSZ);
2017
s = palloc(len + 1);
2018
memcpy(s, VARDATA(str), len);
2021
result = DirectFunctionCall3(numeric_in, CStringGetDatum(s),
2022
ObjectIdGetDatum(0), Int32GetDatum(-1));
2030
numeric_text(PG_FUNCTION_ARGS)
2032
/* val is numeric, but easier to leave it as Datum */
2033
Datum val = PG_GETARG_DATUM(0);
2038
s = DatumGetCString(DirectFunctionCall1(numeric_out, val));
2041
result = (text *) palloc(VARHDRSZ + len);
2043
VARATT_SIZEP(result) = len + VARHDRSZ;
2044
memcpy(VARDATA(result), s, len);
2048
PG_RETURN_TEXT_P(result);
2052
/* ----------------------------------------------------------------------
2054
* Aggregate functions
2056
* The transition datatype for all these aggregates is a 3-element array
2057
* of Numeric, holding the values N, sum(X), sum(X*X) in that order.
2059
* We represent N as a numeric mainly to avoid having to build a special
2060
* datatype; it's unlikely it'd overflow an int4, but ...
2062
* ----------------------------------------------------------------------
2066
do_numeric_accum(ArrayType *transarray, Numeric newval)
2075
/* We assume the input is array of numeric */
2076
deconstruct_array(transarray,
2077
NUMERICOID, -1, false, 'i',
2078
&transdatums, &ndatums);
2080
elog(ERROR, "expected 3-element numeric array");
2082
sumX = transdatums[1];
2083
sumX2 = transdatums[2];
2085
N = DirectFunctionCall1(numeric_inc, N);
2086
sumX = DirectFunctionCall2(numeric_add, sumX,
2087
NumericGetDatum(newval));
2088
sumX2 = DirectFunctionCall2(numeric_add, sumX2,
2089
DirectFunctionCall2(numeric_mul,
2090
NumericGetDatum(newval),
2091
NumericGetDatum(newval)));
2094
transdatums[1] = sumX;
2095
transdatums[2] = sumX2;
2097
result = construct_array(transdatums, 3,
2098
NUMERICOID, -1, false, 'i');
2104
numeric_accum(PG_FUNCTION_ARGS)
2106
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
2107
Numeric newval = PG_GETARG_NUMERIC(1);
2109
PG_RETURN_ARRAYTYPE_P(do_numeric_accum(transarray, newval));
2113
* Integer data types all use Numeric accumulators to share code and
2114
* avoid risk of overflow. For int2 and int4 inputs, Numeric accumulation
2115
* is overkill for the N and sum(X) values, but definitely not overkill
2116
* for the sum(X*X) value. Hence, we use int2_accum and int4_accum only
2117
* for stddev/variance --- there are faster special-purpose accumulator
2118
* routines for SUM and AVG of these datatypes.
2122
int2_accum(PG_FUNCTION_ARGS)
2124
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
2125
Datum newval2 = PG_GETARG_DATUM(1);
2128
newval = DatumGetNumeric(DirectFunctionCall1(int2_numeric, newval2));
2130
PG_RETURN_ARRAYTYPE_P(do_numeric_accum(transarray, newval));
2134
int4_accum(PG_FUNCTION_ARGS)
2136
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
2137
Datum newval4 = PG_GETARG_DATUM(1);
2140
newval = DatumGetNumeric(DirectFunctionCall1(int4_numeric, newval4));
2142
PG_RETURN_ARRAYTYPE_P(do_numeric_accum(transarray, newval));
2146
int8_accum(PG_FUNCTION_ARGS)
2148
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
2149
Datum newval8 = PG_GETARG_DATUM(1);
2152
newval = DatumGetNumeric(DirectFunctionCall1(int8_numeric, newval8));
2154
PG_RETURN_ARRAYTYPE_P(do_numeric_accum(transarray, newval));
2158
numeric_avg(PG_FUNCTION_ARGS)
2160
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
2166
/* We assume the input is array of numeric */
2167
deconstruct_array(transarray,
2168
NUMERICOID, -1, false, 'i',
2169
&transdatums, &ndatums);
2171
elog(ERROR, "expected 3-element numeric array");
2172
N = DatumGetNumeric(transdatums[0]);
2173
sumX = DatumGetNumeric(transdatums[1]);
2176
/* SQL92 defines AVG of no values to be NULL */
2177
/* N is zero iff no digits (cf. numeric_uminus) */
2178
if (N->varlen == NUMERIC_HDRSZ)
2181
PG_RETURN_DATUM(DirectFunctionCall2(numeric_div,
2182
NumericGetDatum(sumX),
2183
NumericGetDatum(N)));
2187
numeric_variance(PG_FUNCTION_ARGS)
2189
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
2202
/* We assume the input is array of numeric */
2203
deconstruct_array(transarray,
2204
NUMERICOID, -1, false, 'i',
2205
&transdatums, &ndatums);
2207
elog(ERROR, "expected 3-element numeric array");
2208
N = DatumGetNumeric(transdatums[0]);
2209
sumX = DatumGetNumeric(transdatums[1]);
2210
sumX2 = DatumGetNumeric(transdatums[2]);
2212
if (NUMERIC_IS_NAN(N) || NUMERIC_IS_NAN(sumX) || NUMERIC_IS_NAN(sumX2))
2213
PG_RETURN_NUMERIC(make_result(&const_nan));
2215
/* Sample variance is undefined when N is 0 or 1, so return NULL */
2217
set_var_from_num(N, &vN);
2219
if (cmp_var(&vN, &const_one) <= 0)
2225
init_var(&vNminus1);
2226
sub_var(&vN, &const_one, &vNminus1);
2229
set_var_from_num(sumX, &vsumX);
2231
set_var_from_num(sumX2, &vsumX2);
2233
/* compute rscale for mul_var calls */
2234
rscale = vsumX.dscale * 2;
2236
mul_var(&vsumX, &vsumX, &vsumX, rscale); /* vsumX = sumX * sumX */
2237
mul_var(&vN, &vsumX2, &vsumX2, rscale); /* vsumX2 = N * sumX2 */
2238
sub_var(&vsumX2, &vsumX, &vsumX2); /* N * sumX2 - sumX * sumX */
2240
if (cmp_var(&vsumX2, &const_zero) <= 0)
2242
/* Watch out for roundoff error producing a negative numerator */
2243
res = make_result(&const_zero);
2247
mul_var(&vN, &vNminus1, &vNminus1, 0); /* N * (N - 1) */
2248
rscale = select_div_scale(&vsumX2, &vNminus1);
2249
div_var(&vsumX2, &vNminus1, &vsumX, rscale); /* variance */
2251
res = make_result(&vsumX);
2255
free_var(&vNminus1);
2259
PG_RETURN_NUMERIC(res);
2263
numeric_stddev(PG_FUNCTION_ARGS)
2265
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
2278
/* We assume the input is array of numeric */
2279
deconstruct_array(transarray,
2280
NUMERICOID, -1, false, 'i',
2281
&transdatums, &ndatums);
2283
elog(ERROR, "expected 3-element numeric array");
2284
N = DatumGetNumeric(transdatums[0]);
2285
sumX = DatumGetNumeric(transdatums[1]);
2286
sumX2 = DatumGetNumeric(transdatums[2]);
2288
if (NUMERIC_IS_NAN(N) || NUMERIC_IS_NAN(sumX) || NUMERIC_IS_NAN(sumX2))
2289
PG_RETURN_NUMERIC(make_result(&const_nan));
2291
/* Sample stddev is undefined when N is 0 or 1, so return NULL */
2293
set_var_from_num(N, &vN);
2295
if (cmp_var(&vN, &const_one) <= 0)
2301
init_var(&vNminus1);
2302
sub_var(&vN, &const_one, &vNminus1);
2305
set_var_from_num(sumX, &vsumX);
2307
set_var_from_num(sumX2, &vsumX2);
2309
/* compute rscale for mul_var calls */
2310
rscale = vsumX.dscale * 2;
2312
mul_var(&vsumX, &vsumX, &vsumX, rscale); /* vsumX = sumX * sumX */
2313
mul_var(&vN, &vsumX2, &vsumX2, rscale); /* vsumX2 = N * sumX2 */
2314
sub_var(&vsumX2, &vsumX, &vsumX2); /* N * sumX2 - sumX * sumX */
2316
if (cmp_var(&vsumX2, &const_zero) <= 0)
2318
/* Watch out for roundoff error producing a negative numerator */
2319
res = make_result(&const_zero);
2323
mul_var(&vN, &vNminus1, &vNminus1, 0); /* N * (N - 1) */
2324
rscale = select_div_scale(&vsumX2, &vNminus1);
2325
div_var(&vsumX2, &vNminus1, &vsumX, rscale); /* variance */
2326
sqrt_var(&vsumX, &vsumX, rscale); /* stddev */
2328
res = make_result(&vsumX);
2332
free_var(&vNminus1);
2336
PG_RETURN_NUMERIC(res);
2341
* SUM transition functions for integer datatypes.
2343
* To avoid overflow, we use accumulators wider than the input datatype.
2344
* A Numeric accumulator is needed for int8 input; for int4 and int2
2345
* inputs, we use int8 accumulators which should be sufficient for practical
2346
* purposes. (The latter two therefore don't really belong in this file,
2347
* but we keep them here anyway.)
2349
* Because SQL92 defines the SUM() of no values to be NULL, not zero,
2350
* the initial condition of the transition data value needs to be NULL. This
2351
* means we can't rely on ExecAgg to automatically insert the first non-null
2352
* data value into the transition data: it doesn't know how to do the type
2353
* conversion. The upshot is that these routines have to be marked non-strict
2354
* and handle substitution of the first non-null input themselves.
2358
int2_sum(PG_FUNCTION_ARGS)
2363
if (PG_ARGISNULL(0))
2365
/* No non-null input seen so far... */
2366
if (PG_ARGISNULL(1))
2367
PG_RETURN_NULL(); /* still no non-null */
2368
/* This is the first non-null input. */
2369
newval = (int64) PG_GETARG_INT16(1);
2370
PG_RETURN_INT64(newval);
2373
oldsum = PG_GETARG_INT64(0);
2375
/* Leave sum unchanged if new input is null. */
2376
if (PG_ARGISNULL(1))
2377
PG_RETURN_INT64(oldsum);
2379
/* OK to do the addition. */
2380
newval = oldsum + (int64) PG_GETARG_INT16(1);
2382
PG_RETURN_INT64(newval);
2386
int4_sum(PG_FUNCTION_ARGS)
2391
if (PG_ARGISNULL(0))
2393
/* No non-null input seen so far... */
2394
if (PG_ARGISNULL(1))
2395
PG_RETURN_NULL(); /* still no non-null */
2396
/* This is the first non-null input. */
2397
newval = (int64) PG_GETARG_INT32(1);
2398
PG_RETURN_INT64(newval);
2401
oldsum = PG_GETARG_INT64(0);
2403
/* Leave sum unchanged if new input is null. */
2404
if (PG_ARGISNULL(1))
2405
PG_RETURN_INT64(oldsum);
2407
/* OK to do the addition. */
2408
newval = oldsum + (int64) PG_GETARG_INT32(1);
2410
PG_RETURN_INT64(newval);
2414
int8_sum(PG_FUNCTION_ARGS)
2419
if (PG_ARGISNULL(0))
2421
/* No non-null input seen so far... */
2422
if (PG_ARGISNULL(1))
2423
PG_RETURN_NULL(); /* still no non-null */
2424
/* This is the first non-null input. */
2425
newval = DirectFunctionCall1(int8_numeric, PG_GETARG_DATUM(1));
2426
PG_RETURN_DATUM(newval);
2429
oldsum = PG_GETARG_NUMERIC(0);
2431
/* Leave sum unchanged if new input is null. */
2432
if (PG_ARGISNULL(1))
2433
PG_RETURN_NUMERIC(oldsum);
2435
/* OK to do the addition. */
2436
newval = DirectFunctionCall1(int8_numeric, PG_GETARG_DATUM(1));
2438
PG_RETURN_DATUM(DirectFunctionCall2(numeric_add,
2439
NumericGetDatum(oldsum), newval));
2444
* Routines for avg(int2) and avg(int4). The transition datatype
2445
* is a two-element int8 array, holding count and sum.
2448
typedef struct Int8TransTypeData
2450
#ifndef INT64_IS_BUSTED
2454
/* "int64" isn't really 64 bits, so fake up properly-aligned fields */
2460
} Int8TransTypeData;
2463
int2_avg_accum(PG_FUNCTION_ARGS)
2465
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P_COPY(0);
2466
int16 newval = PG_GETARG_INT16(1);
2467
Int8TransTypeData *transdata;
2470
* We copied the input array, so it's okay to scribble on it directly.
2472
if (ARR_SIZE(transarray) != ARR_OVERHEAD(1) + sizeof(Int8TransTypeData))
2473
elog(ERROR, "expected 2-element int8 array");
2474
transdata = (Int8TransTypeData *) ARR_DATA_PTR(transarray);
2477
transdata->sum += newval;
2479
PG_RETURN_ARRAYTYPE_P(transarray);
2483
int4_avg_accum(PG_FUNCTION_ARGS)
2485
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P_COPY(0);
2486
int32 newval = PG_GETARG_INT32(1);
2487
Int8TransTypeData *transdata;
2490
* We copied the input array, so it's okay to scribble on it directly.
2492
if (ARR_SIZE(transarray) != ARR_OVERHEAD(1) + sizeof(Int8TransTypeData))
2493
elog(ERROR, "expected 2-element int8 array");
2494
transdata = (Int8TransTypeData *) ARR_DATA_PTR(transarray);
2497
transdata->sum += newval;
2499
PG_RETURN_ARRAYTYPE_P(transarray);
2503
int8_avg(PG_FUNCTION_ARGS)
2505
ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0);
2506
Int8TransTypeData *transdata;
2510
if (ARR_SIZE(transarray) != ARR_OVERHEAD(1) + sizeof(Int8TransTypeData))
2511
elog(ERROR, "expected 2-element int8 array");
2512
transdata = (Int8TransTypeData *) ARR_DATA_PTR(transarray);
2514
/* SQL92 defines AVG of no values to be NULL */
2515
if (transdata->count == 0)
2518
countd = DirectFunctionCall1(int8_numeric,
2519
Int64GetDatumFast(transdata->count));
2520
sumd = DirectFunctionCall1(int8_numeric,
2521
Int64GetDatumFast(transdata->sum));
2523
PG_RETURN_DATUM(DirectFunctionCall2(numeric_div, sumd, countd));
2527
/* ----------------------------------------------------------------------
2531
* ----------------------------------------------------------------------
2534
#ifdef NUMERIC_DEBUG
2537
* dump_numeric() - Dump a value in the db storage format for debugging
2540
dump_numeric(const char *str, Numeric num)
2542
NumericDigit *digits = (NumericDigit *) num->n_data;
2546
ndigits = (num->varlen - NUMERIC_HDRSZ) / sizeof(NumericDigit);
2548
printf("%s: NUMERIC w=%d d=%d ", str, num->n_weight, NUMERIC_DSCALE(num));
2549
switch (NUMERIC_SIGN(num))
2561
printf("SIGN=0x%x", NUMERIC_SIGN(num));
2565
for (i = 0; i < ndigits; i++)
2566
printf(" %0*d", DEC_DIGITS, digits[i]);
2572
* dump_var() - Dump a value in the variable format for debugging
2575
dump_var(const char *str, NumericVar *var)
2579
printf("%s: VAR w=%d d=%d ", str, var->weight, var->dscale);
2592
printf("SIGN=0x%x", var->sign);
2596
for (i = 0; i < var->ndigits; i++)
2597
printf(" %0*d", DEC_DIGITS, var->digits[i]);
2601
#endif /* NUMERIC_DEBUG */
2604
/* ----------------------------------------------------------------------
2606
* Local functions follow
2608
* In general, these do not support NaNs --- callers must eliminate
2609
* the possibility of NaN first. (make_result() is an exception.)
2611
* ----------------------------------------------------------------------
2618
* Allocate a digit buffer of ndigits digits (plus a spare digit for rounding)
2621
alloc_var(NumericVar *var, int ndigits)
2623
digitbuf_free(var->buf);
2624
var->buf = digitbuf_alloc(ndigits + 1);
2625
var->buf[0] = 0; /* spare digit for rounding */
2626
var->digits = var->buf + 1;
2627
var->ndigits = ndigits;
2634
* Return the digit buffer of a variable to the free pool
2637
free_var(NumericVar *var)
2639
digitbuf_free(var->buf);
2642
var->sign = NUMERIC_NAN;
2649
* Set a variable to ZERO.
2650
* Note: its dscale is not touched.
2653
zero_var(NumericVar *var)
2655
digitbuf_free(var->buf);
2659
var->weight = 0; /* by convention; doesn't really matter */
2660
var->sign = NUMERIC_POS; /* anything but NAN... */
2665
* set_var_from_str()
2667
* Parse a string and put the number into a variable
2670
set_var_from_str(const char *str, NumericVar *dest)
2672
const char *cp = str;
2673
bool have_dp = FALSE;
2675
unsigned char *decdigits;
2676
int sign = NUMERIC_POS;
2683
NumericDigit *digits;
2686
* We first parse the string to extract decimal digits and determine
2687
* the correct decimal weight. Then convert to NBASE representation.
2690
/* skip leading spaces */
2693
if (!isspace((unsigned char) *cp))
2717
if (!isdigit((unsigned char) *cp))
2719
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2720
errmsg("invalid input syntax for type numeric: \"%s\"", str)));
2722
decdigits = (unsigned char *) palloc(strlen(cp) + DEC_DIGITS * 2);
2724
/* leading padding for digit alignment later */
2725
memset(decdigits, 0, DEC_DIGITS);
2730
if (isdigit((unsigned char) *cp))
2732
decdigits[i++] = *cp++ - '0';
2738
else if (*cp == '.')
2742
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2743
errmsg("invalid input syntax for type numeric: \"%s\"",
2752
ddigits = i - DEC_DIGITS;
2753
/* trailing padding for digit alignment later */
2754
memset(decdigits + i, 0, DEC_DIGITS - 1);
2756
/* Handle exponent, if any */
2757
if (*cp == 'e' || *cp == 'E')
2763
exponent = strtol(cp, &endptr, 10);
2766
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2767
errmsg("invalid input syntax for type numeric: \"%s\"",
2770
if (exponent > NUMERIC_MAX_PRECISION ||
2771
exponent < -NUMERIC_MAX_PRECISION)
2773
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2774
errmsg("invalid input syntax for type numeric: \"%s\"",
2776
dweight += (int) exponent;
2777
dscale -= (int) exponent;
2782
/* Should be nothing left but spaces */
2785
if (!isspace((unsigned char) *cp))
2787
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2788
errmsg("invalid input syntax for type numeric: \"%s\"",
2794
* Okay, convert pure-decimal representation to base NBASE. First we
2795
* need to determine the converted weight and ndigits. offset is the
2796
* number of decimal zeroes to insert before the first given digit to
2797
* have a correctly aligned first NBASE digit.
2800
weight = (dweight + 1 + DEC_DIGITS - 1) / DEC_DIGITS - 1;
2802
weight = -((-dweight - 1) / DEC_DIGITS + 1);
2803
offset = (weight + 1) * DEC_DIGITS - (dweight + 1);
2804
ndigits = (ddigits + offset + DEC_DIGITS - 1) / DEC_DIGITS;
2806
alloc_var(dest, ndigits);
2808
dest->weight = weight;
2809
dest->dscale = dscale;
2811
i = DEC_DIGITS - offset;
2812
digits = dest->digits;
2814
while (ndigits-- > 0)
2817
*digits++ = ((decdigits[i] * 10 + decdigits[i + 1]) * 10 +
2818
decdigits[i + 2]) * 10 + decdigits[i + 3];
2819
#elif DEC_DIGITS == 2
2820
*digits++ = decdigits[i] * 10 + decdigits[i + 1];
2821
#elif DEC_DIGITS == 1
2822
*digits++ = decdigits[i];
2824
#error unsupported NBASE
2831
/* Strip any leading/trailing zeroes, and normalize weight if zero */
2837
* set_var_from_num() -
2839
* Convert the packed db format into a variable
2842
set_var_from_num(Numeric num, NumericVar *dest)
2846
ndigits = (num->varlen - NUMERIC_HDRSZ) / sizeof(NumericDigit);
2848
alloc_var(dest, ndigits);
2850
dest->weight = num->n_weight;
2851
dest->sign = NUMERIC_SIGN(num);
2852
dest->dscale = NUMERIC_DSCALE(num);
2854
memcpy(dest->digits, num->n_data, ndigits * sizeof(NumericDigit));
2859
* set_var_from_var() -
2861
* Copy one variable into another
2864
set_var_from_var(NumericVar *value, NumericVar *dest)
2866
NumericDigit *newbuf;
2868
newbuf = digitbuf_alloc(value->ndigits + 1);
2869
newbuf[0] = 0; /* spare digit for rounding */
2870
memcpy(newbuf + 1, value->digits, value->ndigits * sizeof(NumericDigit));
2872
digitbuf_free(dest->buf);
2874
memmove(dest, value, sizeof(NumericVar));
2876
dest->digits = newbuf + 1;
2881
* get_str_from_var() -
2883
* Convert a var to text representation (guts of numeric_out).
2884
* CAUTION: var's contents may be modified by rounding!
2885
* Returns a palloc'd string.
2888
get_str_from_var(NumericVar *var, int dscale)
2905
* Check if we must round up before printing the value and do so.
2907
round_var(var, dscale);
2910
* Allocate space for the result.
2912
* i is set to to # of decimal digits before decimal point. dscale is the
2913
* # of decimal digits we will print after decimal point. We may
2914
* generate as many as DEC_DIGITS-1 excess digits at the end, and in
2915
* addition we need room for sign, decimal point, null terminator.
2917
i = (var->weight + 1) * DEC_DIGITS;
2921
str = palloc(i + dscale + DEC_DIGITS + 2);
2925
* Output a dash for negative values
2927
if (var->sign == NUMERIC_NEG)
2931
* Output all digits before the decimal point
2933
if (var->weight < 0)
2935
d = var->weight + 1;
2940
for (d = 0; d <= var->weight; d++)
2942
dig = (d < var->ndigits) ? var->digits[d] : 0;
2943
/* In the first digit, suppress extra leading decimal zeroes */
2946
bool putit = (d > 0);
2965
#elif DEC_DIGITS == 2
2968
if (d1 > 0 || d > 0)
2971
#elif DEC_DIGITS == 1
2974
#error unsupported NBASE
2980
* If requested, output a decimal point and all the digits that follow
2981
* it. We initially put out a multiple of DEC_DIGITS digits, then
2982
* truncate if needed.
2987
endcp = cp + dscale;
2988
for (i = 0; i < dscale; d++, i += DEC_DIGITS)
2990
dig = (d >= 0 && d < var->ndigits) ? var->digits[d] : 0;
3002
#elif DEC_DIGITS == 2
3007
#elif DEC_DIGITS == 1
3010
#error unsupported NBASE
3017
* terminate the string and return it
3027
* Create the packed db numeric format in palloc()'d memory from
3031
make_result(NumericVar *var)
3034
NumericDigit *digits = var->digits;
3035
int weight = var->weight;
3036
int sign = var->sign;
3040
if (sign == NUMERIC_NAN)
3042
result = (Numeric) palloc(NUMERIC_HDRSZ);
3044
result->varlen = NUMERIC_HDRSZ;
3045
result->n_weight = 0;
3046
result->n_sign_dscale = NUMERIC_NAN;
3048
dump_numeric("make_result()", result);
3054
/* truncate leading zeroes */
3055
while (n > 0 && *digits == 0)
3061
/* truncate trailing zeroes */
3062
while (n > 0 && digits[n - 1] == 0)
3065
/* If zero result, force to weight=0 and positive sign */
3072
/* Build the result */
3073
len = NUMERIC_HDRSZ + n * sizeof(NumericDigit);
3074
result = (Numeric) palloc(len);
3075
result->varlen = len;
3076
result->n_weight = weight;
3077
result->n_sign_dscale = sign | (var->dscale & NUMERIC_DSCALE_MASK);
3079
memcpy(result->n_data, digits, n * sizeof(NumericDigit));
3081
/* Check for overflow of int16 fields */
3082
if (result->n_weight != weight ||
3083
NUMERIC_DSCALE(result) != var->dscale)
3085
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
3086
errmsg("value overflows numeric format")));
3088
dump_numeric("make_result()", result);
3096
* Do bounds checking and rounding according to the attributes
3100
apply_typmod(NumericVar *var, int32 typmod)
3108
/* Do nothing if we have a default typmod (-1) */
3109
if (typmod < (int32) (VARHDRSZ))
3113
precision = (typmod >> 16) & 0xffff;
3114
scale = typmod & 0xffff;
3115
maxdigits = precision - scale;
3117
/* Round to target scale (and set var->dscale) */
3118
round_var(var, scale);
3121
* Check for overflow - note we can't do this before rounding, because
3122
* rounding could raise the weight. Also note that the var's weight
3123
* could be inflated by leading zeroes, which will be stripped before
3124
* storage but perhaps might not have been yet. In any case, we must
3125
* recognize a true zero, whose weight doesn't mean anything.
3127
ddigits = (var->weight + 1) * DEC_DIGITS;
3128
if (ddigits > maxdigits)
3130
/* Determine true weight; and check for all-zero result */
3131
for (i = 0; i < var->ndigits; i++)
3133
NumericDigit dig = var->digits[i];
3137
/* Adjust for any high-order decimal zero digits */
3143
else if (dig < 1000)
3145
#elif DEC_DIGITS == 2
3148
#elif DEC_DIGITS == 1
3151
#error unsupported NBASE
3153
if (ddigits > maxdigits)
3155
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
3156
errmsg("numeric field overflow"),
3157
errdetail("The absolute value is greater than or equal to 10^%d for field with precision %d, scale %d.",
3158
ddigits - 1, precision, scale)));
3161
ddigits -= DEC_DIGITS;
3167
* Convert numeric to int8, rounding if needed.
3169
* If overflow, return FALSE (no error is raised). Return TRUE if okay.
3171
* CAUTION: var's contents may be modified by rounding!
3174
numericvar_to_int8(NumericVar *var, int64 *result)
3176
NumericDigit *digits;
3184
/* Round to nearest integer */
3187
/* Check for zero input */
3189
ndigits = var->ndigits;
3197
* For input like 10000000000, we must treat stripped digits as real.
3198
* So the loop assumes there are weight+1 digits before the decimal
3201
weight = var->weight;
3202
Assert(weight >= 0 && ndigits <= weight + 1);
3204
/* Construct the result */
3205
digits = var->digits;
3206
neg = (var->sign == NUMERIC_NEG);
3208
for (i = 1; i <= weight; i++)
3216
* The overflow check is a bit tricky because we want to accept
3217
* INT64_MIN, which will overflow the positive accumulator. We
3218
* can detect this case easily though because INT64_MIN is the
3219
* only nonzero value for which -val == val (on a two's complement
3222
if ((val / NBASE) != oldval) /* possible overflow? */
3224
if (!neg || (-val) != val || val == 0 || oldval < 0)
3229
*result = neg ? -val : val;
3234
* Convert int8 value to numeric.
3237
int8_to_numericvar(int64 val, NumericVar *var)
3244
/* int8 can require at most 19 decimal digits; add one for safety */
3245
alloc_var(var, 20 / DEC_DIGITS);
3248
var->sign = NUMERIC_NEG;
3253
var->sign = NUMERIC_POS;
3263
ptr = var->digits + var->ndigits;
3269
newuval = uval / NBASE;
3270
*ptr = uval - newuval * NBASE;
3274
var->ndigits = ndigits;
3275
var->weight = ndigits - 1;
3279
* Convert numeric to float8; if out of range, return +/- HUGE_VAL
3282
numeric_to_double_no_overflow(Numeric num)
3288
tmp = DatumGetCString(DirectFunctionCall1(numeric_out,
3289
NumericGetDatum(num)));
3291
/* unlike float8in, we ignore ERANGE from strtod */
3292
val = strtod(tmp, &endptr);
3293
if (*endptr != '\0')
3295
/* shouldn't happen ... */
3297
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
3298
errmsg("invalid input syntax for type double precision: \"%s\"",
3307
/* As above, but work from a NumericVar */
3309
numericvar_to_double_no_overflow(NumericVar *var)
3315
tmp = get_str_from_var(var, var->dscale);
3317
/* unlike float8in, we ignore ERANGE from strtod */
3318
val = strtod(tmp, &endptr);
3319
if (*endptr != '\0')
3321
/* shouldn't happen ... */
3323
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
3324
errmsg("invalid input syntax for type double precision: \"%s\"",
3337
* Compare two values on variable level. We assume zeroes have been
3338
* truncated to no digits.
3341
cmp_var(NumericVar *var1, NumericVar *var2)
3343
if (var1->ndigits == 0)
3345
if (var2->ndigits == 0)
3347
if (var2->sign == NUMERIC_NEG)
3351
if (var2->ndigits == 0)
3353
if (var1->sign == NUMERIC_POS)
3358
if (var1->sign == NUMERIC_POS)
3360
if (var2->sign == NUMERIC_NEG)
3362
return cmp_abs(var1, var2);
3365
if (var2->sign == NUMERIC_POS)
3368
return cmp_abs(var2, var1);
3375
* Full version of add functionality on variable level (handling signs).
3376
* result might point to one of the operands too without danger.
3379
add_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
3382
* Decide on the signs of the two variables what to do
3384
if (var1->sign == NUMERIC_POS)
3386
if (var2->sign == NUMERIC_POS)
3389
* Both are positive result = +(ABS(var1) + ABS(var2))
3391
add_abs(var1, var2, result);
3392
result->sign = NUMERIC_POS;
3397
* var1 is positive, var2 is negative Must compare absolute
3400
switch (cmp_abs(var1, var2))
3404
* ABS(var1) == ABS(var2)
3409
result->dscale = Max(var1->dscale, var2->dscale);
3414
* ABS(var1) > ABS(var2)
3415
* result = +(ABS(var1) - ABS(var2))
3418
sub_abs(var1, var2, result);
3419
result->sign = NUMERIC_POS;
3424
* ABS(var1) < ABS(var2)
3425
* result = -(ABS(var2) - ABS(var1))
3428
sub_abs(var2, var1, result);
3429
result->sign = NUMERIC_NEG;
3436
if (var2->sign == NUMERIC_POS)
3439
* var1 is negative, var2 is positive
3440
* Must compare absolute values
3443
switch (cmp_abs(var1, var2))
3447
* ABS(var1) == ABS(var2)
3452
result->dscale = Max(var1->dscale, var2->dscale);
3457
* ABS(var1) > ABS(var2)
3458
* result = -(ABS(var1) - ABS(var2))
3461
sub_abs(var1, var2, result);
3462
result->sign = NUMERIC_NEG;
3467
* ABS(var1) < ABS(var2)
3468
* result = +(ABS(var2) - ABS(var1))
3471
sub_abs(var2, var1, result);
3472
result->sign = NUMERIC_POS;
3480
* result = -(ABS(var1) + ABS(var2))
3483
add_abs(var1, var2, result);
3484
result->sign = NUMERIC_NEG;
3493
* Full version of sub functionality on variable level (handling signs).
3494
* result might point to one of the operands too without danger.
3497
sub_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
3500
* Decide on the signs of the two variables what to do
3502
if (var1->sign == NUMERIC_POS)
3504
if (var2->sign == NUMERIC_NEG)
3507
* var1 is positive, var2 is negative
3508
* result = +(ABS(var1) + ABS(var2))
3511
add_abs(var1, var2, result);
3512
result->sign = NUMERIC_POS;
3518
* Must compare absolute values
3521
switch (cmp_abs(var1, var2))
3525
* ABS(var1) == ABS(var2)
3530
result->dscale = Max(var1->dscale, var2->dscale);
3535
* ABS(var1) > ABS(var2)
3536
* result = +(ABS(var1) - ABS(var2))
3539
sub_abs(var1, var2, result);
3540
result->sign = NUMERIC_POS;
3545
* ABS(var1) < ABS(var2)
3546
* result = -(ABS(var2) - ABS(var1))
3549
sub_abs(var2, var1, result);
3550
result->sign = NUMERIC_NEG;
3557
if (var2->sign == NUMERIC_NEG)
3561
* Must compare absolute values
3564
switch (cmp_abs(var1, var2))
3568
* ABS(var1) == ABS(var2)
3573
result->dscale = Max(var1->dscale, var2->dscale);
3578
* ABS(var1) > ABS(var2)
3579
* result = -(ABS(var1) - ABS(var2))
3582
sub_abs(var1, var2, result);
3583
result->sign = NUMERIC_NEG;
3588
* ABS(var1) < ABS(var2)
3589
* result = +(ABS(var2) - ABS(var1))
3592
sub_abs(var2, var1, result);
3593
result->sign = NUMERIC_POS;
3600
* var1 is negative, var2 is positive
3601
* result = -(ABS(var1) + ABS(var2))
3604
add_abs(var1, var2, result);
3605
result->sign = NUMERIC_NEG;
3614
* Multiplication on variable level. Product of var1 * var2 is stored
3615
* in result. Result is rounded to no more than rscale fractional digits.
3618
mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
3629
NumericDigit *res_digits;
3635
/* copy these values into local vars for speed in inner loop */
3636
int var1ndigits = var1->ndigits;
3637
int var2ndigits = var2->ndigits;
3638
NumericDigit *var1digits = var1->digits;
3639
NumericDigit *var2digits = var2->digits;
3641
if (var1ndigits == 0 || var2ndigits == 0)
3643
/* one or both inputs is zero; so is result */
3645
result->dscale = rscale;
3649
/* Determine result sign and (maximum possible) weight */
3650
if (var1->sign == var2->sign)
3651
res_sign = NUMERIC_POS;
3653
res_sign = NUMERIC_NEG;
3654
res_weight = var1->weight + var2->weight + 2;
3657
* Determine number of result digits to compute. If the exact result
3658
* would have more than rscale fractional digits, truncate the
3659
* computation with MUL_GUARD_DIGITS guard digits. We do that by
3660
* pretending that one or both inputs have fewer digits than they
3663
res_ndigits = var1ndigits + var2ndigits + 1;
3664
maxdigits = res_weight + 1 + (rscale * DEC_DIGITS) + MUL_GUARD_DIGITS;
3665
if (res_ndigits > maxdigits)
3669
/* no useful precision at all in the result... */
3671
result->dscale = rscale;
3674
/* force maxdigits odd so that input ndigits can be equal */
3675
if ((maxdigits & 1) == 0)
3677
if (var1ndigits > var2ndigits)
3679
var1ndigits -= res_ndigits - maxdigits;
3680
if (var1ndigits < var2ndigits)
3681
var1ndigits = var2ndigits = (var1ndigits + var2ndigits) / 2;
3685
var2ndigits -= res_ndigits - maxdigits;
3686
if (var2ndigits < var1ndigits)
3687
var1ndigits = var2ndigits = (var1ndigits + var2ndigits) / 2;
3689
res_ndigits = maxdigits;
3690
Assert(res_ndigits == var1ndigits + var2ndigits + 1);
3694
* We do the arithmetic in an array "dig[]" of signed int's. Since
3695
* INT_MAX is noticeably larger than NBASE*NBASE, this gives us
3696
* headroom to avoid normalizing carries immediately.
3698
* maxdig tracks the maximum possible value of any dig[] entry; when this
3699
* threatens to exceed INT_MAX, we take the time to propagate carries.
3700
* To avoid overflow in maxdig itself, it actually represents the max
3701
* possible value divided by NBASE-1.
3703
dig = (int *) palloc0(res_ndigits * sizeof(int));
3706
ri = res_ndigits - 1;
3707
for (i1 = var1ndigits - 1; i1 >= 0; ri--, i1--)
3709
int var1digit = var1digits[i1];
3714
/* Time to normalize? */
3715
maxdig += var1digit;
3716
if (maxdig > INT_MAX / (NBASE - 1))
3720
for (i = res_ndigits - 1; i >= 0; i--)
3722
newdig = dig[i] + carry;
3723
if (newdig >= NBASE)
3725
carry = newdig / NBASE;
3726
newdig -= carry * NBASE;
3733
/* Reset maxdig to indicate new worst-case */
3734
maxdig = 1 + var1digit;
3737
/* Add appropriate multiple of var2 into the accumulator */
3739
for (i2 = var2ndigits - 1; i2 >= 0; i2--)
3740
dig[i--] += var1digit * var2digits[i2];
3744
* Now we do a final carry propagation pass to normalize the result,
3745
* which we combine with storing the result digits into the output.
3746
* Note that this is still done at full precision w/guard digits.
3748
alloc_var(result, res_ndigits);
3749
res_digits = result->digits;
3751
for (i = res_ndigits - 1; i >= 0; i--)
3753
newdig = dig[i] + carry;
3754
if (newdig >= NBASE)
3756
carry = newdig / NBASE;
3757
newdig -= carry * NBASE;
3761
res_digits[i] = newdig;
3768
* Finally, round the result to the requested precision.
3770
result->weight = res_weight;
3771
result->sign = res_sign;
3773
/* Round to target rscale (and set result->dscale) */
3774
round_var(result, rscale);
3776
/* Strip leading and trailing zeroes */
3784
* Division on variable level. Quotient of var1 / var2 is stored
3785
* in result. Result is rounded to no more than rscale fractional digits.
3788
div_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
3799
NumericDigit *res_digits;
3807
/* copy these values into local vars for speed in inner loop */
3808
int var1ndigits = var1->ndigits;
3809
int var2ndigits = var2->ndigits;
3810
NumericDigit *var1digits = var1->digits;
3811
NumericDigit *var2digits = var2->digits;
3814
* First of all division by zero check; we must not be handed an
3815
* unnormalized divisor.
3817
if (var2ndigits == 0 || var2digits[0] == 0)
3819
(errcode(ERRCODE_DIVISION_BY_ZERO),
3820
errmsg("division by zero")));
3823
* Now result zero check
3825
if (var1ndigits == 0)
3828
result->dscale = rscale;
3833
* Determine the result sign, weight and number of digits to calculate
3835
if (var1->sign == var2->sign)
3836
res_sign = NUMERIC_POS;
3838
res_sign = NUMERIC_NEG;
3839
res_weight = var1->weight - var2->weight + 1;
3840
/* The number of accurate result digits we need to produce: */
3841
div_ndigits = res_weight + 1 + (rscale + DEC_DIGITS - 1) / DEC_DIGITS;
3842
/* Add guard digits for roundoff error */
3843
div_ndigits += DIV_GUARD_DIGITS;
3844
if (div_ndigits < DIV_GUARD_DIGITS)
3845
div_ndigits = DIV_GUARD_DIGITS;
3846
/* Must be at least var1ndigits, too, to simplify data-loading loop */
3847
if (div_ndigits < var1ndigits)
3848
div_ndigits = var1ndigits;
3851
* We do the arithmetic in an array "div[]" of signed int's. Since
3852
* INT_MAX is noticeably larger than NBASE*NBASE, this gives us
3853
* headroom to avoid normalizing carries immediately.
3855
* We start with div[] containing one zero digit followed by the
3856
* dividend's digits (plus appended zeroes to reach the desired
3857
* precision including guard digits). Each step of the main loop
3858
* computes an (approximate) quotient digit and stores it into div[],
3859
* removing one position of dividend space. A final pass of carry
3860
* propagation takes care of any mistaken quotient digits.
3862
div = (int *) palloc0((div_ndigits + 1) * sizeof(int));
3863
for (i = 0; i < var1ndigits; i++)
3864
div[i + 1] = var1digits[i];
3867
* We estimate each quotient digit using floating-point arithmetic,
3868
* taking the first four digits of the (current) dividend and divisor.
3869
* This must be float to avoid overflow.
3871
fdivisor = (double) var2digits[0];
3872
for (i = 1; i < 4; i++)
3875
if (i < var2ndigits)
3876
fdivisor += (double) var2digits[i];
3878
fdivisorinverse = 1.0 / fdivisor;
3881
* maxdiv tracks the maximum possible absolute value of any div[]
3882
* entry; when this threatens to exceed INT_MAX, we take the time to
3883
* propagate carries. To avoid overflow in maxdiv itself, it actually
3884
* represents the max possible abs. value divided by NBASE-1.
3889
* Outer loop computes next quotient digit, which will go into div[qi]
3891
for (qi = 0; qi < div_ndigits; qi++)
3893
/* Approximate the current dividend value */
3894
fdividend = (double) div[qi];
3895
for (i = 1; i < 4; i++)
3898
if (qi + i <= div_ndigits)
3899
fdividend += (double) div[qi + i];
3901
/* Compute the (approximate) quotient digit */
3902
fquotient = fdividend * fdivisorinverse;
3903
qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
3904
(((int) fquotient) - 1); /* truncate towards -infinity */
3908
/* Do we need to normalize now? */
3909
maxdiv += Abs(qdigit);
3910
if (maxdiv > INT_MAX / (NBASE - 1))
3914
for (i = div_ndigits; i > qi; i--)
3916
newdig = div[i] + carry;
3919
carry = -((-newdig - 1) / NBASE) - 1;
3920
newdig -= carry * NBASE;
3922
else if (newdig >= NBASE)
3924
carry = newdig / NBASE;
3925
newdig -= carry * NBASE;
3931
newdig = div[qi] + carry;
3935
* All the div[] digits except possibly div[qi] are now in
3936
* the range 0..NBASE-1.
3938
maxdiv = Abs(newdig) / (NBASE - 1);
3939
maxdiv = Max(maxdiv, 1);
3942
* Recompute the quotient digit since new info may have
3943
* propagated into the top four dividend digits
3945
fdividend = (double) div[qi];
3946
for (i = 1; i < 4; i++)
3949
if (qi + i <= div_ndigits)
3950
fdividend += (double) div[qi + i];
3952
/* Compute the (approximate) quotient digit */
3953
fquotient = fdividend * fdivisorinverse;
3954
qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
3955
(((int) fquotient) - 1); /* truncate towards
3957
maxdiv += Abs(qdigit);
3960
/* Subtract off the appropriate multiple of the divisor */
3963
int istop = Min(var2ndigits, div_ndigits - qi + 1);
3965
for (i = 0; i < istop; i++)
3966
div[qi + i] -= qdigit * var2digits[i];
3971
* The dividend digit we are about to replace might still be
3972
* nonzero. Fold it into the next digit position. We don't need
3973
* to worry about overflow here since this should nearly cancel
3974
* with the subtraction of the divisor.
3976
div[qi + 1] += div[qi] * NBASE;
3982
* Approximate and store the last quotient digit (div[div_ndigits])
3984
fdividend = (double) div[qi];
3985
for (i = 1; i < 4; i++)
3987
fquotient = fdividend * fdivisorinverse;
3988
qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
3989
(((int) fquotient) - 1); /* truncate towards -infinity */
3993
* Now we do a final carry propagation pass to normalize the result,
3994
* which we combine with storing the result digits into the output.
3995
* Note that this is still done at full precision w/guard digits.
3997
alloc_var(result, div_ndigits + 1);
3998
res_digits = result->digits;
4000
for (i = div_ndigits; i >= 0; i--)
4002
newdig = div[i] + carry;
4005
carry = -((-newdig - 1) / NBASE) - 1;
4006
newdig -= carry * NBASE;
4008
else if (newdig >= NBASE)
4010
carry = newdig / NBASE;
4011
newdig -= carry * NBASE;
4015
res_digits[i] = newdig;
4022
* Finally, round the result to the requested precision.
4024
result->weight = res_weight;
4025
result->sign = res_sign;
4027
/* Round to target rscale (and set result->dscale) */
4028
round_var(result, rscale);
4030
/* Strip leading and trailing zeroes */
4036
* Default scale selection for division
4038
* Returns the appropriate result scale for the division result.
4041
select_div_scale(NumericVar *var1, NumericVar *var2)
4047
NumericDigit firstdigit1,
4052
* The result scale of a division isn't specified in any SQL standard.
4053
* For PostgreSQL we select a result scale that will give at least
4054
* NUMERIC_MIN_SIG_DIGITS significant digits, so that numeric gives a
4055
* result no less accurate than float8; but use a scale not less than
4056
* either input's display scale.
4059
/* Get the actual (normalized) weight and first digit of each input */
4061
weight1 = 0; /* values to use if var1 is zero */
4063
for (i = 0; i < var1->ndigits; i++)
4065
firstdigit1 = var1->digits[i];
4066
if (firstdigit1 != 0)
4068
weight1 = var1->weight - i;
4073
weight2 = 0; /* values to use if var2 is zero */
4075
for (i = 0; i < var2->ndigits; i++)
4077
firstdigit2 = var2->digits[i];
4078
if (firstdigit2 != 0)
4080
weight2 = var2->weight - i;
4086
* Estimate weight of quotient. If the two first digits are equal, we
4087
* can't be sure, but assume that var1 is less than var2.
4089
qweight = weight1 - weight2;
4090
if (firstdigit1 <= firstdigit2)
4093
/* Select result scale */
4094
rscale = NUMERIC_MIN_SIG_DIGITS - qweight * DEC_DIGITS;
4095
rscale = Max(rscale, var1->dscale);
4096
rscale = Max(rscale, var2->dscale);
4097
rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
4098
rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
4107
* Calculate the modulo of two numerics at variable level
4110
mod_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
4118
* We do this using the equation
4119
* mod(x,y) = x - trunc(x/y)*y
4120
* We set rscale the same way numeric_div and numeric_mul do
4121
* to get the right answer from the equation. The final result,
4122
* however, need not be displayed to more precision than the inputs.
4125
rscale = select_div_scale(var1, var2);
4127
div_var(var1, var2, &tmp, rscale);
4131
mul_var(var2, &tmp, &tmp, var2->dscale + tmp.dscale);
4133
sub_var(var1, &tmp, result);
4135
round_var(result, Max(var1->dscale, var2->dscale));
4144
* Return the smallest integer greater than or equal to the argument
4148
ceil_var(NumericVar *var, NumericVar *result)
4153
set_var_from_var(var, &tmp);
4157
if (var->sign == NUMERIC_POS && cmp_var(var, &tmp) != 0)
4158
add_var(&tmp, &const_one, &tmp);
4160
set_var_from_var(&tmp, result);
4168
* Return the largest integer equal to or less than the argument
4172
floor_var(NumericVar *var, NumericVar *result)
4177
set_var_from_var(var, &tmp);
4181
if (var->sign == NUMERIC_NEG && cmp_var(var, &tmp) != 0)
4182
sub_var(&tmp, &const_one, &tmp);
4184
set_var_from_var(&tmp, result);
4192
* Compute the square root of x using Newton's algorithm
4195
sqrt_var(NumericVar *arg, NumericVar *result, int rscale)
4199
NumericVar last_val;
4203
local_rscale = rscale + 8;
4205
stat = cmp_var(arg, &const_zero);
4209
result->dscale = rscale;
4214
* SQL2003 defines sqrt() in terms of power, so we need to emit the
4215
* right SQLSTATE error code if the operand is negative.
4219
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
4220
errmsg("cannot take square root of a negative number")));
4224
init_var(&last_val);
4226
/* Copy arg in case it is the same var as result */
4227
set_var_from_var(arg, &tmp_arg);
4230
* Initialize the result to the first guess
4232
alloc_var(result, 1);
4233
result->digits[0] = tmp_arg.digits[0] / 2;
4234
if (result->digits[0] == 0)
4235
result->digits[0] = 1;
4236
result->weight = tmp_arg.weight / 2;
4237
result->sign = NUMERIC_POS;
4239
set_var_from_var(result, &last_val);
4243
div_var(&tmp_arg, result, &tmp_val, local_rscale);
4245
add_var(result, &tmp_val, result);
4246
mul_var(result, &const_zero_point_five, result, local_rscale);
4248
if (cmp_var(&last_val, result) == 0)
4250
set_var_from_var(result, &last_val);
4253
free_var(&last_val);
4257
/* Round to requested precision */
4258
round_var(result, rscale);
4265
* Raise e to the power of x
4268
exp_var(NumericVar *arg, NumericVar *result, int rscale)
4276
* We separate the integral and fraction parts of x, then compute
4277
* e^x = e^xint * e^xfrac
4278
* where e = exp(1) and e^xfrac = exp(xfrac) are computed by
4279
* exp_var_internal; the limited range of inputs allows that routine
4280
* to do a good job with a simple Taylor series. Raising e^xint is
4281
* done by repeated multiplications in power_var_int.
4286
set_var_from_var(arg, &x);
4288
if (x.sign == NUMERIC_NEG)
4291
x.sign = NUMERIC_POS;
4294
/* Extract the integer part, remove it from x */
4296
while (x.weight >= 0)
4301
xintval += x.digits[0];
4306
/* Guard against overflow */
4307
if (xintval >= NUMERIC_MAX_RESULT_SCALE * 3)
4309
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
4310
errmsg("argument for function \"exp\" too big")));
4313
/* Select an appropriate scale for internal calculation */
4314
local_rscale = rscale + MUL_GUARD_DIGITS * 2;
4316
/* Compute e^xfrac */
4317
exp_var_internal(&x, result, local_rscale);
4319
/* If there's an integer part, multiply by e^xint */
4325
exp_var_internal(&const_one, &e, local_rscale);
4326
power_var_int(&e, xintval, &e, local_rscale);
4327
mul_var(&e, result, result, local_rscale);
4331
/* Compensate for input sign, and round to requested rscale */
4333
div_var(&const_one, result, result, rscale);
4335
round_var(result, rscale);
4342
* exp_var_internal() -
4344
* Raise e to the power of x, where 0 <= x <= 1
4346
* NB: the result should be good to at least rscale digits, but it has
4347
* *not* been rounded off; the caller must do that if wanted.
4350
exp_var_internal(NumericVar *arg, NumericVar *result, int rscale)
4366
set_var_from_var(arg, &x);
4368
Assert(x.sign == NUMERIC_POS);
4370
local_rscale = rscale + 8;
4372
/* Reduce input into range 0 <= x <= 0.01 */
4373
while (cmp_var(&x, &const_zero_point_01) > 0)
4377
mul_var(&x, &const_zero_point_five, &x, x.dscale + 1);
4381
* Use the Taylor series
4383
* exp(x) = 1 + x + x^2/2! + x^3/3! + ...
4385
* Given the limited range of x, this should converge reasonably quickly.
4386
* We run the series until the terms fall below the local_rscale
4389
add_var(&const_one, &x, result);
4390
set_var_from_var(&x, &xpow);
4391
set_var_from_var(&const_one, &ifac);
4392
set_var_from_var(&const_one, &ni);
4396
add_var(&ni, &const_one, &ni);
4397
mul_var(&xpow, &x, &xpow, local_rscale);
4398
mul_var(&ifac, &ni, &ifac, 0);
4399
div_var(&xpow, &ifac, &elem, local_rscale);
4401
if (elem.ndigits == 0)
4404
add_var(result, &elem, result);
4407
/* Compensate for argument range reduction */
4409
mul_var(result, result, result, local_rscale);
4422
* Compute the natural log of x
4425
ln_var(NumericVar *arg, NumericVar *result, int rscale)
4435
cmp = cmp_var(arg, &const_zero);
4438
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
4439
errmsg("cannot take logarithm of zero")));
4442
(errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
4443
errmsg("cannot take logarithm of a negative number")));
4445
local_rscale = rscale + 8;
4453
set_var_from_var(arg, &x);
4454
set_var_from_var(&const_two, &fact);
4456
/* Reduce input into range 0.9 < x < 1.1 */
4457
while (cmp_var(&x, &const_zero_point_nine) <= 0)
4460
sqrt_var(&x, &x, local_rscale);
4461
mul_var(&fact, &const_two, &fact, 0);
4463
while (cmp_var(&x, &const_one_point_one) >= 0)
4466
sqrt_var(&x, &x, local_rscale);
4467
mul_var(&fact, &const_two, &fact, 0);
4471
* We use the Taylor series for 0.5 * ln((1+z)/(1-z)),
4473
* z + z^3/3 + z^5/5 + ...
4475
* where z = (x-1)/(x+1) is in the range (approximately) -0.053 .. 0.048
4476
* due to the above range-reduction of x.
4478
* The convergence of this is not as fast as one would like, but is
4479
* tolerable given that z is small.
4481
sub_var(&x, &const_one, result);
4482
add_var(&x, &const_one, &elem);
4483
div_var(result, &elem, result, local_rscale);
4484
set_var_from_var(result, &xx);
4485
mul_var(result, result, &x, local_rscale);
4487
set_var_from_var(&const_one, &ni);
4491
add_var(&ni, &const_two, &ni);
4492
mul_var(&xx, &x, &xx, local_rscale);
4493
div_var(&xx, &ni, &elem, local_rscale);
4495
if (elem.ndigits == 0)
4498
add_var(result, &elem, result);
4500
if (elem.weight < (result->weight - local_rscale * 2 / DEC_DIGITS))
4504
/* Compensate for argument range reduction, round to requested rscale */
4505
mul_var(result, &fact, result, rscale);
4518
* Compute the logarithm of num in a given base.
4520
* Note: this routine chooses dscale of the result.
4523
log_var(NumericVar *base, NumericVar *num, NumericVar *result)
4534
/* Set scale for ln() calculations --- compare numeric_ln() */
4536
/* Approx decimal digits before decimal point */
4537
dec_digits = (num->weight + 1) * DEC_DIGITS;
4540
rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(dec_digits - 1);
4541
else if (dec_digits < 1)
4542
rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(1 - dec_digits);
4544
rscale = NUMERIC_MIN_SIG_DIGITS;
4546
rscale = Max(rscale, base->dscale);
4547
rscale = Max(rscale, num->dscale);
4548
rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
4549
rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
4551
local_rscale = rscale + 8;
4553
/* Form natural logarithms */
4554
ln_var(base, &ln_base, local_rscale);
4555
ln_var(num, &ln_num, local_rscale);
4557
ln_base.dscale = rscale;
4558
ln_num.dscale = rscale;
4560
/* Select scale for division result */
4561
rscale = select_div_scale(&ln_num, &ln_base);
4563
div_var(&ln_num, &ln_base, result, rscale);
4573
* Raise base to the power of exp
4575
* Note: this routine chooses dscale of the result.
4578
power_var(NumericVar *base, NumericVar *exp, NumericVar *result)
4587
/* If exp can be represented as an integer, use power_var_int */
4588
if (exp->ndigits == 0 || exp->ndigits <= exp->weight + 1)
4590
/* exact integer, but does it fit in int? */
4594
/* must copy because numericvar_to_int8() scribbles on input */
4596
set_var_from_var(exp, &x);
4597
if (numericvar_to_int8(&x, &expval64))
4599
int expval = (int) expval64;
4601
/* Test for overflow by reverse-conversion. */
4602
if ((int64) expval == expval64)
4604
/* Okay, select rscale */
4605
rscale = NUMERIC_MIN_SIG_DIGITS;
4606
rscale = Max(rscale, base->dscale);
4607
rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
4608
rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
4610
power_var_int(base, expval, result, rscale);
4622
/* Set scale for ln() calculation --- need extra accuracy here */
4624
/* Approx decimal digits before decimal point */
4625
dec_digits = (base->weight + 1) * DEC_DIGITS;
4628
rscale = NUMERIC_MIN_SIG_DIGITS * 2 - (int) log10(dec_digits - 1);
4629
else if (dec_digits < 1)
4630
rscale = NUMERIC_MIN_SIG_DIGITS * 2 - (int) log10(1 - dec_digits);
4632
rscale = NUMERIC_MIN_SIG_DIGITS * 2;
4634
rscale = Max(rscale, base->dscale * 2);
4635
rscale = Max(rscale, exp->dscale * 2);
4636
rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE * 2);
4637
rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE * 2);
4639
local_rscale = rscale + 8;
4641
ln_var(base, &ln_base, local_rscale);
4643
mul_var(&ln_base, exp, &ln_num, local_rscale);
4645
/* Set scale for exp() -- compare numeric_exp() */
4647
/* convert input to float8, ignoring overflow */
4648
val = numericvar_to_double_no_overflow(&ln_num);
4651
* log10(result) = num * log10(e), so this is approximately the
4654
val *= 0.434294481903252;
4656
/* limit to something that won't cause integer overflow */
4657
val = Max(val, -NUMERIC_MAX_RESULT_SCALE);
4658
val = Min(val, NUMERIC_MAX_RESULT_SCALE);
4660
rscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
4661
rscale = Max(rscale, base->dscale);
4662
rscale = Max(rscale, exp->dscale);
4663
rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
4664
rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
4666
exp_var(&ln_num, result, rscale);
4675
* Raise base to the power of exp, where exp is an integer.
4678
power_var_int(NumericVar *base, int exp, NumericVar *result, int rscale)
4681
NumericVar base_prod;
4684
/* Detect some special cases, particularly 0^0. */
4689
if (base->ndigits == 0)
4691
(errcode(ERRCODE_FLOATING_POINT_EXCEPTION),
4692
errmsg("zero raised to zero is undefined")));
4693
set_var_from_var(&const_one, result);
4694
result->dscale = rscale; /* no need to round */
4697
set_var_from_var(base, result);
4698
round_var(result, rscale);
4701
div_var(&const_one, base, result, rscale);
4704
mul_var(base, base, result, rscale);
4711
* The general case repeatedly multiplies base according to the bit
4712
* pattern of exp. We do the multiplications with some extra
4718
local_rscale = rscale + MUL_GUARD_DIGITS * 2;
4720
init_var(&base_prod);
4721
set_var_from_var(base, &base_prod);
4724
set_var_from_var(base, result);
4726
set_var_from_var(&const_one, result);
4728
while ((exp >>= 1) > 0)
4730
mul_var(&base_prod, &base_prod, &base_prod, local_rscale);
4732
mul_var(&base_prod, result, result, local_rscale);
4735
free_var(&base_prod);
4737
/* Compensate for input sign, and round to requested rscale */
4739
div_var(&const_one, result, result, rscale);
4741
round_var(result, rscale);
4745
/* ----------------------------------------------------------------------
4747
* Following are the lowest level functions that operate unsigned
4748
* on the variable level
4750
* ----------------------------------------------------------------------
4757
* Compare the absolute values of var1 and var2
4758
* Returns: -1 for ABS(var1) < ABS(var2)
4759
* 0 for ABS(var1) == ABS(var2)
4760
* 1 for ABS(var1) > ABS(var2)
4764
cmp_abs(NumericVar *var1, NumericVar *var2)
4766
NumericDigit *var1digits = var1->digits;
4767
NumericDigit *var2digits = var2->digits;
4770
int w1 = var1->weight;
4771
int w2 = var2->weight;
4773
/* Check any digits before the first common digit */
4775
while (w1 > w2 && i1 < var1->ndigits)
4777
if (var1digits[i1++] != 0)
4781
while (w2 > w1 && i2 < var2->ndigits)
4783
if (var2digits[i2++] != 0)
4788
/* At this point, either w1 == w2 or we've run out of digits */
4792
while (i1 < var1->ndigits && i2 < var2->ndigits)
4794
int stat = var1digits[i1++] - var2digits[i2++];
4806
* At this point, we've run out of digits on one side or the other; so
4807
* any remaining nonzero digits imply that side is larger
4809
while (i1 < var1->ndigits)
4811
if (var1digits[i1++] != 0)
4814
while (i2 < var2->ndigits)
4816
if (var2digits[i2++] != 0)
4827
* Add the absolute values of two variables into result.
4828
* result might point to one of the operands without danger.
4831
add_abs(NumericVar *var1, NumericVar *var2, NumericVar *result)
4833
NumericDigit *res_buf;
4834
NumericDigit *res_digits;
4846
/* copy these values into local vars for speed in inner loop */
4847
int var1ndigits = var1->ndigits;
4848
int var2ndigits = var2->ndigits;
4849
NumericDigit *var1digits = var1->digits;
4850
NumericDigit *var2digits = var2->digits;
4852
res_weight = Max(var1->weight, var2->weight) + 1;
4854
res_dscale = Max(var1->dscale, var2->dscale);
4856
/* Note: here we are figuring rscale in base-NBASE digits */
4857
rscale1 = var1->ndigits - var1->weight - 1;
4858
rscale2 = var2->ndigits - var2->weight - 1;
4859
res_rscale = Max(rscale1, rscale2);
4861
res_ndigits = res_rscale + res_weight + 1;
4862
if (res_ndigits <= 0)
4865
res_buf = digitbuf_alloc(res_ndigits + 1);
4866
res_buf[0] = 0; /* spare digit for later rounding */
4867
res_digits = res_buf + 1;
4869
i1 = res_rscale + var1->weight + 1;
4870
i2 = res_rscale + var2->weight + 1;
4871
for (i = res_ndigits - 1; i >= 0; i--)
4875
if (i1 >= 0 && i1 < var1ndigits)
4876
carry += var1digits[i1];
4877
if (i2 >= 0 && i2 < var2ndigits)
4878
carry += var2digits[i2];
4882
res_digits[i] = carry - NBASE;
4887
res_digits[i] = carry;
4892
Assert(carry == 0); /* else we failed to allow for carry out */
4894
digitbuf_free(result->buf);
4895
result->ndigits = res_ndigits;
4896
result->buf = res_buf;
4897
result->digits = res_digits;
4898
result->weight = res_weight;
4899
result->dscale = res_dscale;
4901
/* Remove leading/trailing zeroes */
4909
* Subtract the absolute value of var2 from the absolute value of var1
4910
* and store in result. result might point to one of the operands
4913
* ABS(var1) MUST BE GREATER OR EQUAL ABS(var2) !!!
4916
sub_abs(NumericVar *var1, NumericVar *var2, NumericVar *result)
4918
NumericDigit *res_buf;
4919
NumericDigit *res_digits;
4931
/* copy these values into local vars for speed in inner loop */
4932
int var1ndigits = var1->ndigits;
4933
int var2ndigits = var2->ndigits;
4934
NumericDigit *var1digits = var1->digits;
4935
NumericDigit *var2digits = var2->digits;
4937
res_weight = var1->weight;
4939
res_dscale = Max(var1->dscale, var2->dscale);
4941
/* Note: here we are figuring rscale in base-NBASE digits */
4942
rscale1 = var1->ndigits - var1->weight - 1;
4943
rscale2 = var2->ndigits - var2->weight - 1;
4944
res_rscale = Max(rscale1, rscale2);
4946
res_ndigits = res_rscale + res_weight + 1;
4947
if (res_ndigits <= 0)
4950
res_buf = digitbuf_alloc(res_ndigits + 1);
4951
res_buf[0] = 0; /* spare digit for later rounding */
4952
res_digits = res_buf + 1;
4954
i1 = res_rscale + var1->weight + 1;
4955
i2 = res_rscale + var2->weight + 1;
4956
for (i = res_ndigits - 1; i >= 0; i--)
4960
if (i1 >= 0 && i1 < var1ndigits)
4961
borrow += var1digits[i1];
4962
if (i2 >= 0 && i2 < var2ndigits)
4963
borrow -= var2digits[i2];
4967
res_digits[i] = borrow + NBASE;
4972
res_digits[i] = borrow;
4977
Assert(borrow == 0); /* else caller gave us var1 < var2 */
4979
digitbuf_free(result->buf);
4980
result->ndigits = res_ndigits;
4981
result->buf = res_buf;
4982
result->digits = res_digits;
4983
result->weight = res_weight;
4984
result->dscale = res_dscale;
4986
/* Remove leading/trailing zeroes */
4993
* Round the value of a variable to no more than rscale decimal digits
4994
* after the decimal point. NOTE: we allow rscale < 0 here, implying
4995
* rounding before the decimal point.
4998
round_var(NumericVar *var, int rscale)
5000
NumericDigit *digits = var->digits;
5005
var->dscale = rscale;
5007
/* decimal digits wanted */
5008
di = (var->weight + 1) * DEC_DIGITS + rscale;
5011
* If di = 0, the value loses all digits, but could round up to 1 if
5012
* its first extra digit is >= 5. If di < 0 the result must be 0.
5018
var->sign = NUMERIC_POS;
5022
/* NBASE digits wanted */
5023
ndigits = (di + DEC_DIGITS - 1) / DEC_DIGITS;
5025
/* 0, or number of decimal digits to keep in last NBASE digit */
5028
if (ndigits < var->ndigits ||
5029
(ndigits == var->ndigits && di > 0))
5031
var->ndigits = ndigits;
5034
/* di must be zero */
5035
carry = (digits[ndigits] >= HALF_NBASE) ? 1 : 0;
5038
carry = (digits[ndigits] >= HALF_NBASE) ? 1 : 0;
5041
/* Must round within last NBASE digit */
5046
pow10 = round_powers[di];
5047
#elif DEC_DIGITS == 2
5050
#error unsupported NBASE
5052
extra = digits[--ndigits] % pow10;
5053
digits[ndigits] -= extra;
5055
if (extra >= pow10 / 2)
5057
pow10 += digits[ndigits];
5063
digits[ndigits] = pow10;
5068
/* Propagate carry if needed */
5071
carry += digits[--ndigits];
5074
digits[ndigits] = carry - NBASE;
5079
digits[ndigits] = carry;
5086
Assert(ndigits == -1); /* better not have added > 1 digit */
5087
Assert(var->digits > var->buf);
5099
* Truncate the value of a variable at rscale decimal digits after the
5100
* decimal point. NOTE: we allow rscale < 0 here, implying
5101
* truncation before the decimal point.
5104
trunc_var(NumericVar *var, int rscale)
5109
var->dscale = rscale;
5111
/* decimal digits wanted */
5112
di = (var->weight + 1) * DEC_DIGITS + rscale;
5115
* If di <= 0, the value loses all digits.
5121
var->sign = NUMERIC_POS;
5125
/* NBASE digits wanted */
5126
ndigits = (di + DEC_DIGITS - 1) / DEC_DIGITS;
5128
if (ndigits <= var->ndigits)
5130
var->ndigits = ndigits;
5133
/* no within-digit stuff to worry about */
5135
/* 0, or number of decimal digits to keep in last NBASE digit */
5140
/* Must truncate within last NBASE digit */
5141
NumericDigit *digits = var->digits;
5146
pow10 = round_powers[di];
5147
#elif DEC_DIGITS == 2
5150
#error unsupported NBASE
5152
extra = digits[--ndigits] % pow10;
5153
digits[ndigits] -= extra;
5163
* Strip any leading and trailing zeroes from a numeric variable
5166
strip_var(NumericVar *var)
5168
NumericDigit *digits = var->digits;
5169
int ndigits = var->ndigits;
5171
/* Strip leading zeroes */
5172
while (ndigits > 0 && *digits == 0)
5179
/* Strip trailing zeroes */
5180
while (ndigits > 0 && digits[ndigits - 1] == 0)
5183
/* If it's zero, normalize the sign and weight */
5186
var->sign = NUMERIC_POS;
5190
var->digits = digits;
5191
var->ndigits = ndigits;