granicus.if.org Git - postgresql/blob - src/backend/utils/adt/numeric.c

   1 /*-------------------------------------------------------------------------
   2  *
   3  * numeric.c
   4  *        An exact numeric data type for the Postgres database system
   5  *
   6  * Original coding 1998, Jan Wieck.  Heavily revised 2003, Tom Lane.
   7  *
   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,
  12  * pages 359-367.
  13  *
  14  * Copyright (c) 1998-2006, PostgreSQL Global Development Group
  15  *
  16  * IDENTIFICATION
  17  *        $PostgreSQL: pgsql/src/backend/utils/adt/numeric.c,v 1.95 2006/10/03 21:25:55 momjian Exp $
  18  *
  19  *-------------------------------------------------------------------------
  20  */
  21
  22 #include "postgres.h"
  23
  24 #include <ctype.h>
  25 #include <float.h>
  26 #include <limits.h>
  27 #include <math.h>
  28
  29 #include "catalog/pg_type.h"
  30 #include "libpq/pqformat.h"
  31 #include "utils/array.h"
  32 #include "utils/builtins.h"
  33 #include "utils/int8.h"
  34 #include "utils/numeric.h"
  35
  36 /* ----------
  37  * Uncomment the following to enable compilation of dump_numeric()
  38  * and dump_var() and to get a dump of any result produced by make_result().
  39  * ----------
  40 #define NUMERIC_DEBUG
  41  */
  42
  43
  44 /* ----------
  45  * Local data types
  46  *
  47  * Numeric values are represented in a base-NBASE floating point format.
  48  * Each "digit" ranges from 0 to NBASE-1.  The type NumericDigit is signed
  49  * and wide enough to store a digit.  We assume that NBASE*NBASE can fit in
  50  * an int.      Although the purely calculational routines could handle any even
  51  * NBASE that's less than sqrt(INT_MAX), in practice we are only interested
  52  * in NBASE a power of ten, so that I/O conversions and decimal rounding
  53  * are easy.  Also, it's actually more efficient if NBASE is rather less than
  54  * sqrt(INT_MAX), so that there is "headroom" for mul_var and div_var to
  55  * postpone processing carries.
  56  * ----------
  57  */
  58
  59 #if 0
  60 #define NBASE           10
  61 #define HALF_NBASE      5
  62 #define DEC_DIGITS      1                       /* decimal digits per NBASE digit */
  63 #define MUL_GUARD_DIGITS        4       /* these are measured in NBASE digits */
  64 #define DIV_GUARD_DIGITS        8
  65
  66 typedef signed char NumericDigit;
  67 #endif
  68
  69 #if 0
  70 #define NBASE           100
  71 #define HALF_NBASE      50
  72 #define DEC_DIGITS      2                       /* decimal digits per NBASE digit */
  73 #define MUL_GUARD_DIGITS        3       /* these are measured in NBASE digits */
  74 #define DIV_GUARD_DIGITS        6
  75
  76 typedef signed char NumericDigit;
  77 #endif
  78
  79 #if 1
  80 #define NBASE           10000
  81 #define HALF_NBASE      5000
  82 #define DEC_DIGITS      4                       /* decimal digits per NBASE digit */
  83 #define MUL_GUARD_DIGITS        2       /* these are measured in NBASE digits */
  84 #define DIV_GUARD_DIGITS        4
  85
  86 typedef int16 NumericDigit;
  87 #endif
  88
  89
  90 /* ----------
  91  * The value represented by a NumericVar is determined by the sign, weight,
  92  * ndigits, and digits[] array.
  93  * Note: the first digit of a NumericVar's value is assumed to be multiplied
  94  * by NBASE ** weight.  Another way to say it is that there are weight+1
  95  * digits before the decimal point.  It is possible to have weight < 0.
  96  *
  97  * buf points at the physical start of the palloc'd digit buffer for the
  98  * NumericVar.  digits points at the first digit in actual use (the one
  99  * with the specified weight).  We normally leave an unused digit or two
 100  * (preset to zeroes) between buf and digits, so that there is room to store
 101  * a carry out of the top digit without special pushups.  We just need to
 102  * decrement digits (and increment weight) to make room for the carry digit.
 103  * (There is no such extra space in a numeric value stored in the database,
 104  * only in a NumericVar in memory.)
 105  *
 106  * If buf is NULL then the digit buffer isn't actually palloc'd and should
 107  * not be freed --- see the constants below for an example.
 108  *
 109  * dscale, or display scale, is the nominal precision expressed as number
 110  * of digits after the decimal point (it must always be >= 0 at present).
 111  * dscale may be more than the number of physically stored fractional digits,
 112  * implying that we have suppressed storage of significant trailing zeroes.
 113  * It should never be less than the number of stored digits, since that would
 114  * imply hiding digits that are present.  NOTE that dscale is always expressed
 115  * in *decimal* digits, and so it may correspond to a fractional number of
 116  * base-NBASE digits --- divide by DEC_DIGITS to convert to NBASE digits.
 117  *
 118  * rscale, or result scale, is the target precision for a computation.
 119  * Like dscale it is expressed as number of *decimal* digits after the decimal
 120  * point, and is always >= 0 at present.
 121  * Note that rscale is not stored in variables --- it's figured on-the-fly
 122  * from the dscales of the inputs.
 123  *
 124  * NB: All the variable-level functions are written in a style that makes it
 125  * possible to give one and the same variable as argument and destination.
 126  * This is feasible because the digit buffer is separate from the variable.
 127  * ----------
 128  */
 129 typedef struct NumericVar
 130 {
 131         int                     ndigits;                /* # of digits in digits[] - can be 0! */
 132         int                     weight;                 /* weight of first digit */
 133         int                     sign;                   /* NUMERIC_POS, NUMERIC_NEG, or NUMERIC_NAN */
 134         int                     dscale;                 /* display scale */
 135         NumericDigit *buf;                      /* start of palloc'd space for digits[] */
 136         NumericDigit *digits;           /* base-NBASE digits */
 137 } NumericVar;
 138
 139
 140 /* ----------
 141  * Some preinitialized constants
 142  * ----------
 143  */
 144 static NumericDigit const_zero_data[1] = {0};
 145 static NumericVar const_zero =
 146 {0, 0, NUMERIC_POS, 0, NULL, const_zero_data};
 147
 148 static NumericDigit const_one_data[1] = {1};
 149 static NumericVar const_one =
 150 {1, 0, NUMERIC_POS, 0, NULL, const_one_data};
 151
 152 static NumericDigit const_two_data[1] = {2};
 153 static NumericVar const_two =
 154 {1, 0, NUMERIC_POS, 0, NULL, const_two_data};
 155
 156 #if DEC_DIGITS == 4
 157 static NumericDigit const_zero_point_five_data[1] = {5000};
 158 #elif DEC_DIGITS == 2
 159 static NumericDigit const_zero_point_five_data[1] = {50};
 160 #elif DEC_DIGITS == 1
 161 static NumericDigit const_zero_point_five_data[1] = {5};
 162 #endif
 163 static NumericVar const_zero_point_five =
 164 {1, -1, NUMERIC_POS, 1, NULL, const_zero_point_five_data};
 165
 166 #if DEC_DIGITS == 4
 167 static NumericDigit const_zero_point_nine_data[1] = {9000};
 168 #elif DEC_DIGITS == 2
 169 static NumericDigit const_zero_point_nine_data[1] = {90};
 170 #elif DEC_DIGITS == 1
 171 static NumericDigit const_zero_point_nine_data[1] = {9};
 172 #endif
 173 static NumericVar const_zero_point_nine =
 174 {1, -1, NUMERIC_POS, 1, NULL, const_zero_point_nine_data};
 175
 176 #if DEC_DIGITS == 4
 177 static NumericDigit const_zero_point_01_data[1] = {100};
 178 static NumericVar const_zero_point_01 =
 179 {1, -1, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
 180 #elif DEC_DIGITS == 2
 181 static NumericDigit const_zero_point_01_data[1] = {1};
 182 static NumericVar const_zero_point_01 =
 183 {1, -1, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
 184 #elif DEC_DIGITS == 1
 185 static NumericDigit const_zero_point_01_data[1] = {1};
 186 static NumericVar const_zero_point_01 =
 187 {1, -2, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
 188 #endif
 189
 190 #if DEC_DIGITS == 4
 191 static NumericDigit const_one_point_one_data[2] = {1, 1000};
 192 #elif DEC_DIGITS == 2
 193 static NumericDigit const_one_point_one_data[2] = {1, 10};
 194 #elif DEC_DIGITS == 1
 195 static NumericDigit const_one_point_one_data[2] = {1, 1};
 196 #endif
 197 static NumericVar const_one_point_one =
 198 {2, 0, NUMERIC_POS, 1, NULL, const_one_point_one_data};
 199
 200 static NumericVar const_nan =
 201 {0, 0, NUMERIC_NAN, 0, NULL, NULL};
 202
 203 #if DEC_DIGITS == 4
 204 static const int round_powers[4] = {0, 1000, 100, 10};
 205 #endif
 206
 207
 208 /* ----------
 209  * Local functions
 210  * ----------
 211  */
 212
 213 #ifdef NUMERIC_DEBUG
 214 static void dump_numeric(const char *str, Numeric num);
 215 static void dump_var(const char *str, NumericVar *var);
 216 #else
 217 #define dump_numeric(s,n)
 218 #define dump_var(s,v)
 219 #endif
 220
 221 #define digitbuf_alloc(ndigits)  \
 222         ((NumericDigit *) palloc((ndigits) * sizeof(NumericDigit)))
 223 #define digitbuf_free(buf)      \
 224         do { \
 225                  if ((buf) != NULL) \
 226                          pfree(buf); \
 227         } while (0)
 228
 229 #define init_var(v)             MemSetAligned(v, 0, sizeof(NumericVar))
 230
 231 #define NUMERIC_DIGITS(num) ((NumericDigit *)(num)->n_data)
 232 #define NUMERIC_NDIGITS(num) \
 233         (((num)->varlen - NUMERIC_HDRSZ) / sizeof(NumericDigit))
 234
 235 static void alloc_var(NumericVar *var, int ndigits);
 236 static void free_var(NumericVar *var);
 237 static void zero_var(NumericVar *var);
 238
 239 static void set_var_from_str(const char *str, NumericVar *dest);
 240 static void set_var_from_num(Numeric value, NumericVar *dest);
 241 static void set_var_from_var(NumericVar *value, NumericVar *dest);
 242 static char *get_str_from_var(NumericVar *var, int dscale);
 243
 244 static Numeric make_result(NumericVar *var);
 245
 246 static void apply_typmod(NumericVar *var, int32 typmod);
 247
 248 static int32 numericvar_to_int4(NumericVar *var);
 249 static bool numericvar_to_int8(NumericVar *var, int64 *result);
 250 static void int8_to_numericvar(int64 val, NumericVar *var);
 251 static double numeric_to_double_no_overflow(Numeric num);
 252 static double numericvar_to_double_no_overflow(NumericVar *var);
 253
 254 static int      cmp_numerics(Numeric num1, Numeric num2);
 255 static int      cmp_var(NumericVar *var1, NumericVar *var2);
 256 static int      cmp_var_common(const NumericDigit *var1digits, int var1ndigits,
 257                                                    int var1weight, int var1sign,
 258                                                    const NumericDigit *var2digits, int var2ndigits,
 259                                                    int var2weight, int var2sign);
 260 static void add_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
 261 static void sub_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
 262 static void mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 263                 int rscale);
 264 static void div_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 265                 int rscale, bool round);
 266 static int      select_div_scale(NumericVar *var1, NumericVar *var2);
 267 static void mod_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
 268 static void ceil_var(NumericVar *var, NumericVar *result);
 269 static void floor_var(NumericVar *var, NumericVar *result);
 270
 271 static void sqrt_var(NumericVar *arg, NumericVar *result, int rscale);
 272 static void exp_var(NumericVar *arg, NumericVar *result, int rscale);
 273 static void exp_var_internal(NumericVar *arg, NumericVar *result, int rscale);
 274 static void ln_var(NumericVar *arg, NumericVar *result, int rscale);
 275 static void log_var(NumericVar *base, NumericVar *num, NumericVar *result);
 276 static void power_var(NumericVar *base, NumericVar *exp, NumericVar *result);
 277 static void power_var_int(NumericVar *base, int exp, NumericVar *result,
 278                           int rscale);
 279
 280 static int      cmp_abs(NumericVar *var1, NumericVar *var2);
 281 static int      cmp_abs_common(const NumericDigit *var1digits, int var1ndigits,
 282                                                    int var1weight,
 283                                                    const NumericDigit *var2digits, int var2ndigits,
 284                                                    int var2weight);
 285 static void add_abs(NumericVar *var1, NumericVar *var2, NumericVar *result);
 286 static void sub_abs(NumericVar *var1, NumericVar *var2, NumericVar *result);
 287 static void round_var(NumericVar *var, int rscale);
 288 static void trunc_var(NumericVar *var, int rscale);
 289 static void strip_var(NumericVar *var);
 290 static void compute_bucket(Numeric operand, Numeric bound1, Numeric bound2,
 291                            NumericVar *count_var, NumericVar *result_var);
 292
 293
 294 /* ----------------------------------------------------------------------
 295  *
 296  * Input-, output- and rounding-functions
 297  *
 298  * ----------------------------------------------------------------------
 299  */
 300
 301
 302 /*
 303  * numeric_in() -
 304  *
 305  *      Input function for numeric data type
 306  */
 307 Datum
 308 numeric_in(PG_FUNCTION_ARGS)
 309 {
 310         char       *str = PG_GETARG_CSTRING(0);
 311
 312 #ifdef NOT_USED
 313         Oid                     typelem = PG_GETARG_OID(1);
 314 #endif
 315         int32           typmod = PG_GETARG_INT32(2);
 316         NumericVar      value;
 317         Numeric         res;
 318
 319         /*
 320          * Check for NaN
 321          */
 322         if (pg_strcasecmp(str, "NaN") == 0)
 323                 PG_RETURN_NUMERIC(make_result(&const_nan));
 324
 325         /*
 326          * Use set_var_from_str() to parse the input string and return it in the
 327          * packed DB storage format
 328          */
 329         init_var(&value);
 330         set_var_from_str(str, &value);
 331
 332         apply_typmod(&value, typmod);
 333
 334         res = make_result(&value);
 335         free_var(&value);
 336
 337         PG_RETURN_NUMERIC(res);
 338 }
 339
 340
 341 /*
 342  * numeric_out() -
 343  *
 344  *      Output function for numeric data type
 345  */
 346 Datum
 347 numeric_out(PG_FUNCTION_ARGS)
 348 {
 349         Numeric         num = PG_GETARG_NUMERIC(0);
 350         NumericVar      x;
 351         char       *str;
 352
 353         /*
 354          * Handle NaN
 355          */
 356         if (NUMERIC_IS_NAN(num))
 357                 PG_RETURN_CSTRING(pstrdup("NaN"));
 358
 359         /*
 360          * Get the number in the variable format.
 361          *
 362          * Even if we didn't need to change format, we'd still need to copy the
 363          * value to have a modifiable copy for rounding.  set_var_from_num() also
 364          * guarantees there is extra digit space in case we produce a carry out
 365          * from rounding.
 366          */
 367         init_var(&x);
 368         set_var_from_num(num, &x);
 369
 370         str = get_str_from_var(&x, x.dscale);
 371
 372         free_var(&x);
 373
 374         PG_RETURN_CSTRING(str);
 375 }
 376
 377 /*
 378  *              numeric_recv                    - converts external binary format to numeric
 379  *
 380  * External format is a sequence of int16's:
 381  * ndigits, weight, sign, dscale, NumericDigits.
 382  */
 383 Datum
 384 numeric_recv(PG_FUNCTION_ARGS)
 385 {
 386         StringInfo      buf = (StringInfo) PG_GETARG_POINTER(0);
 387
 388 #ifdef NOT_USED
 389         Oid                     typelem = PG_GETARG_OID(1);
 390 #endif
 391         int32           typmod = PG_GETARG_INT32(2);
 392         NumericVar      value;
 393         Numeric         res;
 394         int                     len,
 395                                 i;
 396
 397         init_var(&value);
 398
 399         len = (uint16) pq_getmsgint(buf, sizeof(uint16));
 400         if (len < 0 || len > NUMERIC_MAX_PRECISION + NUMERIC_MAX_RESULT_SCALE)
 401                 ereport(ERROR,
 402                                 (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
 403                                  errmsg("invalid length in external \"numeric\" value")));
 404
 405         alloc_var(&value, len);
 406
 407         value.weight = (int16) pq_getmsgint(buf, sizeof(int16));
 408         value.sign = (uint16) pq_getmsgint(buf, sizeof(uint16));
 409         if (!(value.sign == NUMERIC_POS ||
 410                   value.sign == NUMERIC_NEG ||
 411                   value.sign == NUMERIC_NAN))
 412                 ereport(ERROR,
 413                                 (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
 414                                  errmsg("invalid sign in external \"numeric\" value")));
 415
 416         value.dscale = (uint16) pq_getmsgint(buf, sizeof(uint16));
 417         for (i = 0; i < len; i++)
 418         {
 419                 NumericDigit d = pq_getmsgint(buf, sizeof(NumericDigit));
 420
 421                 if (d < 0 || d >= NBASE)
 422                         ereport(ERROR,
 423                                         (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
 424                                          errmsg("invalid digit in external \"numeric\" value")));
 425                 value.digits[i] = d;
 426         }
 427
 428         apply_typmod(&value, typmod);
 429
 430         res = make_result(&value);
 431         free_var(&value);
 432
 433         PG_RETURN_NUMERIC(res);
 434 }
 435
 436 /*
 437  *              numeric_send                    - converts numeric to binary format
 438  */
 439 Datum
 440 numeric_send(PG_FUNCTION_ARGS)
 441 {
 442         Numeric         num = PG_GETARG_NUMERIC(0);
 443         NumericVar      x;
 444         StringInfoData buf;
 445         int                     i;
 446
 447         init_var(&x);
 448         set_var_from_num(num, &x);
 449
 450         pq_begintypsend(&buf);
 451
 452         pq_sendint(&buf, x.ndigits, sizeof(int16));
 453         pq_sendint(&buf, x.weight, sizeof(int16));
 454         pq_sendint(&buf, x.sign, sizeof(int16));
 455         pq_sendint(&buf, x.dscale, sizeof(int16));
 456         for (i = 0; i < x.ndigits; i++)
 457                 pq_sendint(&buf, x.digits[i], sizeof(NumericDigit));
 458
 459         free_var(&x);
 460
 461         PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
 462 }
 463
 464
 465 /*
 466  * numeric() -
 467  *
 468  *      This is a special function called by the Postgres database system
 469  *      before a value is stored in a tuple's attribute. The precision and
 470  *      scale of the attribute have to be applied on the value.
 471  */
 472 Datum
 473 numeric         (PG_FUNCTION_ARGS)
 474 {
 475         Numeric         num = PG_GETARG_NUMERIC(0);
 476         int32           typmod = PG_GETARG_INT32(1);
 477         Numeric         new;
 478         int32           tmp_typmod;
 479         int                     precision;
 480         int                     scale;
 481         int                     ddigits;
 482         int                     maxdigits;
 483         NumericVar      var;
 484
 485         /*
 486          * Handle NaN
 487          */
 488         if (NUMERIC_IS_NAN(num))
 489                 PG_RETURN_NUMERIC(make_result(&const_nan));
 490
 491         /*
 492          * If the value isn't a valid type modifier, simply return a copy of the
 493          * input value
 494          */
 495         if (typmod < (int32) (VARHDRSZ))
 496         {
 497                 new = (Numeric) palloc(num->varlen);
 498                 memcpy(new, num, num->varlen);
 499                 PG_RETURN_NUMERIC(new);
 500         }
 501
 502         /*
 503          * Get the precision and scale out of the typmod value
 504          */
 505         tmp_typmod = typmod - VARHDRSZ;
 506         precision = (tmp_typmod >> 16) & 0xffff;
 507         scale = tmp_typmod & 0xffff;
 508         maxdigits = precision - scale;
 509
 510         /*
 511          * If the number is certainly in bounds and due to the target scale no
 512          * rounding could be necessary, just make a copy of the input and modify
 513          * its scale fields.  (Note we assume the existing dscale is honest...)
 514          */
 515         ddigits = (num->n_weight + 1) * DEC_DIGITS;
 516         if (ddigits <= maxdigits && scale >= NUMERIC_DSCALE(num))
 517         {
 518                 new = (Numeric) palloc(num->varlen);
 519                 memcpy(new, num, num->varlen);
 520                 new->n_sign_dscale = NUMERIC_SIGN(new) |
 521                         ((uint16) scale & NUMERIC_DSCALE_MASK);
 522                 PG_RETURN_NUMERIC(new);
 523         }
 524
 525         /*
 526          * We really need to fiddle with things - unpack the number into a
 527          * variable and let apply_typmod() do it.
 528          */
 529         init_var(&var);
 530
 531         set_var_from_num(num, &var);
 532         apply_typmod(&var, typmod);
 533         new = make_result(&var);
 534
 535         free_var(&var);
 536
 537         PG_RETURN_NUMERIC(new);
 538 }
 539
 540
 541 /* ----------------------------------------------------------------------
 542  *
 543  * Sign manipulation, rounding and the like
 544  *
 545  * ----------------------------------------------------------------------
 546  */
 547
 548 Datum
 549 numeric_abs(PG_FUNCTION_ARGS)
 550 {
 551         Numeric         num = PG_GETARG_NUMERIC(0);
 552         Numeric         res;
 553
 554         /*
 555          * Handle NaN
 556          */
 557         if (NUMERIC_IS_NAN(num))
 558                 PG_RETURN_NUMERIC(make_result(&const_nan));
 559
 560         /*
 561          * Do it the easy way directly on the packed format
 562          */
 563         res = (Numeric) palloc(num->varlen);
 564         memcpy(res, num, num->varlen);
 565
 566         res->n_sign_dscale = NUMERIC_POS | NUMERIC_DSCALE(num);
 567
 568         PG_RETURN_NUMERIC(res);
 569 }
 570
 571
 572 Datum
 573 numeric_uminus(PG_FUNCTION_ARGS)
 574 {
 575         Numeric         num = PG_GETARG_NUMERIC(0);
 576         Numeric         res;
 577
 578         /*
 579          * Handle NaN
 580          */
 581         if (NUMERIC_IS_NAN(num))
 582                 PG_RETURN_NUMERIC(make_result(&const_nan));
 583
 584         /*
 585          * Do it the easy way directly on the packed format
 586          */
 587         res = (Numeric) palloc(num->varlen);
 588         memcpy(res, num, num->varlen);
 589
 590         /*
 591          * The packed format is known to be totally zero digit trimmed always. So
 592          * we can identify a ZERO by the fact that there are no digits at all.  Do
 593          * nothing to a zero.
 594          */
 595         if (num->varlen != NUMERIC_HDRSZ)
 596         {
 597                 /* Else, flip the sign */
 598                 if (NUMERIC_SIGN(num) == NUMERIC_POS)
 599                         res->n_sign_dscale = NUMERIC_NEG | NUMERIC_DSCALE(num);
 600                 else
 601                         res->n_sign_dscale = NUMERIC_POS | NUMERIC_DSCALE(num);
 602         }
 603
 604         PG_RETURN_NUMERIC(res);
 605 }
 606
 607
 608 Datum
 609 numeric_uplus(PG_FUNCTION_ARGS)
 610 {
 611         Numeric         num = PG_GETARG_NUMERIC(0);
 612         Numeric         res;
 613
 614         res = (Numeric) palloc(num->varlen);
 615         memcpy(res, num, num->varlen);
 616
 617         PG_RETURN_NUMERIC(res);
 618 }
 619
 620 /*
 621  * numeric_sign() -
 622  *
 623  * returns -1 if the argument is less than 0, 0 if the argument is equal
 624  * to 0, and 1 if the argument is greater than zero.
 625  */
 626 Datum
 627 numeric_sign(PG_FUNCTION_ARGS)
 628 {
 629         Numeric         num = PG_GETARG_NUMERIC(0);
 630         Numeric         res;
 631         NumericVar      result;
 632
 633         /*
 634          * Handle NaN
 635          */
 636         if (NUMERIC_IS_NAN(num))
 637                 PG_RETURN_NUMERIC(make_result(&const_nan));
 638
 639         init_var(&result);
 640
 641         /*
 642          * The packed format is known to be totally zero digit trimmed always. So
 643          * we can identify a ZERO by the fact that there are no digits at all.
 644          */
 645         if (num->varlen == NUMERIC_HDRSZ)
 646                 set_var_from_var(&const_zero, &result);
 647         else
 648         {
 649                 /*
 650                  * And if there are some, we return a copy of ONE with the sign of our
 651                  * argument
 652                  */
 653                 set_var_from_var(&const_one, &result);
 654                 result.sign = NUMERIC_SIGN(num);
 655         }
 656
 657         res = make_result(&result);
 658         free_var(&result);
 659
 660         PG_RETURN_NUMERIC(res);
 661 }
 662
 663
 664 /*
 665  * numeric_round() -
 666  *
 667  *      Round a value to have 'scale' digits after the decimal point.
 668  *      We allow negative 'scale', implying rounding before the decimal
 669  *      point --- Oracle interprets rounding that way.
 670  */
 671 Datum
 672 numeric_round(PG_FUNCTION_ARGS)
 673 {
 674         Numeric         num = PG_GETARG_NUMERIC(0);
 675         int32           scale = PG_GETARG_INT32(1);
 676         Numeric         res;
 677         NumericVar      arg;
 678
 679         /*
 680          * Handle NaN
 681          */
 682         if (NUMERIC_IS_NAN(num))
 683                 PG_RETURN_NUMERIC(make_result(&const_nan));
 684
 685         /*
 686          * Limit the scale value to avoid possible overflow in calculations
 687          */
 688         scale = Max(scale, -NUMERIC_MAX_RESULT_SCALE);
 689         scale = Min(scale, NUMERIC_MAX_RESULT_SCALE);
 690
 691         /*
 692          * Unpack the argument and round it at the proper digit position
 693          */
 694         init_var(&arg);
 695         set_var_from_num(num, &arg);
 696
 697         round_var(&arg, scale);
 698
 699         /* We don't allow negative output dscale */
 700         if (scale < 0)
 701                 arg.dscale = 0;
 702
 703         /*
 704          * Return the rounded result
 705          */
 706         res = make_result(&arg);
 707
 708         free_var(&arg);
 709         PG_RETURN_NUMERIC(res);
 710 }
 711
 712
 713 /*
 714  * numeric_trunc() -
 715  *
 716  *      Truncate a value to have 'scale' digits after the decimal point.
 717  *      We allow negative 'scale', implying a truncation before the decimal
 718  *      point --- Oracle interprets truncation that way.
 719  */
 720 Datum
 721 numeric_trunc(PG_FUNCTION_ARGS)
 722 {
 723         Numeric         num = PG_GETARG_NUMERIC(0);
 724         int32           scale = PG_GETARG_INT32(1);
 725         Numeric         res;
 726         NumericVar      arg;
 727
 728         /*
 729          * Handle NaN
 730          */
 731         if (NUMERIC_IS_NAN(num))
 732                 PG_RETURN_NUMERIC(make_result(&const_nan));
 733
 734         /*
 735          * Limit the scale value to avoid possible overflow in calculations
 736          */
 737         scale = Max(scale, -NUMERIC_MAX_RESULT_SCALE);
 738         scale = Min(scale, NUMERIC_MAX_RESULT_SCALE);
 739
 740         /*
 741          * Unpack the argument and truncate it at the proper digit position
 742          */
 743         init_var(&arg);
 744         set_var_from_num(num, &arg);
 745
 746         trunc_var(&arg, scale);
 747
 748         /* We don't allow negative output dscale */
 749         if (scale < 0)
 750                 arg.dscale = 0;
 751
 752         /*
 753          * Return the truncated result
 754          */
 755         res = make_result(&arg);
 756
 757         free_var(&arg);
 758         PG_RETURN_NUMERIC(res);
 759 }
 760
 761
 762 /*
 763  * numeric_ceil() -
 764  *
 765  *      Return the smallest integer greater than or equal to the argument
 766  */
 767 Datum
 768 numeric_ceil(PG_FUNCTION_ARGS)
 769 {
 770         Numeric         num = PG_GETARG_NUMERIC(0);
 771         Numeric         res;
 772         NumericVar      result;
 773
 774         if (NUMERIC_IS_NAN(num))
 775                 PG_RETURN_NUMERIC(make_result(&const_nan));
 776
 777         init_var(&result);
 778
 779         set_var_from_num(num, &result);
 780         ceil_var(&result, &result);
 781
 782         res = make_result(&result);
 783         free_var(&result);
 784
 785         PG_RETURN_NUMERIC(res);
 786 }
 787
 788
 789 /*
 790  * numeric_floor() -
 791  *
 792  *      Return the largest integer equal to or less than the argument
 793  */
 794 Datum
 795 numeric_floor(PG_FUNCTION_ARGS)
 796 {
 797         Numeric         num = PG_GETARG_NUMERIC(0);
 798         Numeric         res;
 799         NumericVar      result;
 800
 801         if (NUMERIC_IS_NAN(num))
 802                 PG_RETURN_NUMERIC(make_result(&const_nan));
 803
 804         init_var(&result);
 805
 806         set_var_from_num(num, &result);
 807         floor_var(&result, &result);
 808
 809         res = make_result(&result);
 810         free_var(&result);
 811
 812         PG_RETURN_NUMERIC(res);
 813 }
 814
 815 /*
 816  * width_bucket_numeric() -
 817  *
 818  * 'bound1' and 'bound2' are the lower and upper bounds of the
 819  * histogram's range, respectively. 'count' is the number of buckets
 820  * in the histogram. width_bucket() returns an integer indicating the
 821  * bucket number that 'operand' belongs in for an equiwidth histogram
 822  * with the specified characteristics. An operand smaller than the
 823  * lower bound is assigned to bucket 0. An operand greater than the
 824  * upper bound is assigned to an additional bucket (with number
 825  * count+1).
 826  */
 827 Datum
 828 width_bucket_numeric(PG_FUNCTION_ARGS)
 829 {
 830         Numeric         operand = PG_GETARG_NUMERIC(0);
 831         Numeric         bound1 = PG_GETARG_NUMERIC(1);
 832         Numeric         bound2 = PG_GETARG_NUMERIC(2);
 833         int32           count = PG_GETARG_INT32(3);
 834         NumericVar      count_var;
 835         NumericVar      result_var;
 836         int32           result;
 837
 838         if (count <= 0)
 839                 ereport(ERROR,
 840                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
 841                                  errmsg("count must be greater than zero")));
 842
 843         init_var(&result_var);
 844         init_var(&count_var);
 845
 846         /* Convert 'count' to a numeric, for ease of use later */
 847         int8_to_numericvar((int64) count, &count_var);
 848
 849         switch (cmp_numerics(bound1, bound2))
 850         {
 851                 case 0:
 852                         ereport(ERROR,
 853                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
 854                                  errmsg("lower bound cannot equal upper bound")));
 855
 856                         /* bound1 < bound2 */
 857                 case -1:
 858                         if (cmp_numerics(operand, bound1) < 0)
 859                                 set_var_from_var(&const_zero, &result_var);
 860                         else if (cmp_numerics(operand, bound2) >= 0)
 861                                 add_var(&count_var, &const_one, &result_var);
 862                         else
 863                                 compute_bucket(operand, bound1, bound2,
 864                                                            &count_var, &result_var);
 865                         break;
 866
 867                         /* bound1 > bound2 */
 868                 case 1:
 869                         if (cmp_numerics(operand, bound1) > 0)
 870                                 set_var_from_var(&const_zero, &result_var);
 871                         else if (cmp_numerics(operand, bound2) <= 0)
 872                                 add_var(&count_var, &const_one, &result_var);
 873                         else
 874                                 compute_bucket(operand, bound1, bound2,
 875                                                            &count_var, &result_var);
 876                         break;
 877         }
 878
 879         result = numericvar_to_int4(&result_var);
 880
 881         free_var(&count_var);
 882         free_var(&result_var);
 883
 884         PG_RETURN_INT32(result);
 885 }
 886
 887 /*
 888  * compute_bucket() -
 889  *
 890  * If 'operand' is not outside the bucket range, determine the correct
 891  * bucket for it to go. The calculations performed by this function
 892  * are derived directly from the SQL2003 spec.
 893  */
 894 static void
 895 compute_bucket(Numeric operand, Numeric bound1, Numeric bound2,
 896                            NumericVar *count_var, NumericVar *result_var)
 897 {
 898         NumericVar      bound1_var;
 899         NumericVar      bound2_var;
 900         NumericVar      operand_var;
 901
 902         init_var(&bound1_var);
 903         init_var(&bound2_var);
 904         init_var(&operand_var);
 905
 906         set_var_from_num(bound1, &bound1_var);
 907         set_var_from_num(bound2, &bound2_var);
 908         set_var_from_num(operand, &operand_var);
 909
 910         if (cmp_var(&bound1_var, &bound2_var) < 0)
 911         {
 912                 sub_var(&operand_var, &bound1_var, &operand_var);
 913                 sub_var(&bound2_var, &bound1_var, &bound2_var);
 914                 div_var(&operand_var, &bound2_var, result_var,
 915                                 select_div_scale(&operand_var, &bound2_var), true);
 916         }
 917         else
 918         {
 919                 sub_var(&bound1_var, &operand_var, &operand_var);
 920                 sub_var(&bound1_var, &bound2_var, &bound1_var);
 921                 div_var(&operand_var, &bound1_var, result_var,
 922                                 select_div_scale(&operand_var, &bound1_var), true);
 923         }
 924
 925         mul_var(result_var, count_var, result_var,
 926                         result_var->dscale + count_var->dscale);
 927         add_var(result_var, &const_one, result_var);
 928         floor_var(result_var, result_var);
 929
 930         free_var(&bound1_var);
 931         free_var(&bound2_var);
 932         free_var(&operand_var);
 933 }
 934
 935 /* ----------------------------------------------------------------------
 936  *
 937  * Comparison functions
 938  *
 939  * Note: btree indexes need these routines not to leak memory; therefore,
 940  * be careful to free working copies of toasted datums.  Most places don't
 941  * need to be so careful.
 942  * ----------------------------------------------------------------------
 943  */
 944
 945
 946 Datum
 947 numeric_cmp(PG_FUNCTION_ARGS)
 948 {
 949         Numeric         num1 = PG_GETARG_NUMERIC(0);
 950         Numeric         num2 = PG_GETARG_NUMERIC(1);
 951         int                     result;
 952
 953         result = cmp_numerics(num1, num2);
 954
 955         PG_FREE_IF_COPY(num1, 0);
 956         PG_FREE_IF_COPY(num2, 1);
 957
 958         PG_RETURN_INT32(result);
 959 }
 960
 961
 962 Datum
 963 numeric_eq(PG_FUNCTION_ARGS)
 964 {
 965         Numeric         num1 = PG_GETARG_NUMERIC(0);
 966         Numeric         num2 = PG_GETARG_NUMERIC(1);
 967         bool            result;
 968
 969         result = cmp_numerics(num1, num2) == 0;
 970
 971         PG_FREE_IF_COPY(num1, 0);
 972         PG_FREE_IF_COPY(num2, 1);
 973
 974         PG_RETURN_BOOL(result);
 975 }
 976
 977 Datum
 978 numeric_ne(PG_FUNCTION_ARGS)
 979 {
 980         Numeric         num1 = PG_GETARG_NUMERIC(0);
 981         Numeric         num2 = PG_GETARG_NUMERIC(1);
 982         bool            result;
 983
 984         result = cmp_numerics(num1, num2) != 0;
 985
 986         PG_FREE_IF_COPY(num1, 0);
 987         PG_FREE_IF_COPY(num2, 1);
 988
 989         PG_RETURN_BOOL(result);
 990 }
 991
 992 Datum
 993 numeric_gt(PG_FUNCTION_ARGS)
 994 {
 995         Numeric         num1 = PG_GETARG_NUMERIC(0);
 996         Numeric         num2 = PG_GETARG_NUMERIC(1);
 997         bool            result;
 998
 999         result = cmp_numerics(num1, num2) > 0;
1000
1001         PG_FREE_IF_COPY(num1, 0);
1002         PG_FREE_IF_COPY(num2, 1);
1003
1004         PG_RETURN_BOOL(result);
1005 }
1006
1007 Datum
1008 numeric_ge(PG_FUNCTION_ARGS)
1009 {
1010         Numeric         num1 = PG_GETARG_NUMERIC(0);
1011         Numeric         num2 = PG_GETARG_NUMERIC(1);
1012         bool            result;
1013
1014         result = cmp_numerics(num1, num2) >= 0;
1015
1016         PG_FREE_IF_COPY(num1, 0);
1017         PG_FREE_IF_COPY(num2, 1);
1018
1019         PG_RETURN_BOOL(result);
1020 }
1021
1022 Datum
1023 numeric_lt(PG_FUNCTION_ARGS)
1024 {
1025         Numeric         num1 = PG_GETARG_NUMERIC(0);
1026         Numeric         num2 = PG_GETARG_NUMERIC(1);
1027         bool            result;
1028
1029         result = cmp_numerics(num1, num2) < 0;
1030
1031         PG_FREE_IF_COPY(num1, 0);
1032         PG_FREE_IF_COPY(num2, 1);
1033
1034         PG_RETURN_BOOL(result);
1035 }
1036
1037 Datum
1038 numeric_le(PG_FUNCTION_ARGS)
1039 {
1040         Numeric         num1 = PG_GETARG_NUMERIC(0);
1041         Numeric         num2 = PG_GETARG_NUMERIC(1);
1042         bool            result;
1043
1044         result = cmp_numerics(num1, num2) <= 0;
1045
1046         PG_FREE_IF_COPY(num1, 0);
1047         PG_FREE_IF_COPY(num2, 1);
1048
1049         PG_RETURN_BOOL(result);
1050 }
1051
1052 static int
1053 cmp_numerics(Numeric num1, Numeric num2)
1054 {
1055         int                     result;
1056
1057         /*
1058          * We consider all NANs to be equal and larger than any non-NAN. This is
1059          * somewhat arbitrary; the important thing is to have a consistent sort
1060          * order.
1061          */
1062         if (NUMERIC_IS_NAN(num1))
1063         {
1064                 if (NUMERIC_IS_NAN(num2))
1065                         result = 0;                     /* NAN = NAN */
1066                 else
1067                         result = 1;                     /* NAN > non-NAN */
1068         }
1069         else if (NUMERIC_IS_NAN(num2))
1070         {
1071                 result = -1;                    /* non-NAN < NAN */
1072         }
1073         else
1074         {
1075                 result = cmp_var_common(NUMERIC_DIGITS(num1), NUMERIC_NDIGITS(num1),
1076                                                                 num1->n_weight, NUMERIC_SIGN(num1),
1077                                                                 NUMERIC_DIGITS(num2), NUMERIC_NDIGITS(num2),
1078                                                                 num2->n_weight, NUMERIC_SIGN(num2));
1079         }
1080
1081         return result;
1082 }
1083
1084
1085 /* ----------------------------------------------------------------------
1086  *
1087  * Basic arithmetic functions
1088  *
1089  * ----------------------------------------------------------------------
1090  */
1091
1092
1093 /*
1094  * numeric_add() -
1095  *
1096  *      Add two numerics
1097  */
1098 Datum
1099 numeric_add(PG_FUNCTION_ARGS)
1100 {
1101         Numeric         num1 = PG_GETARG_NUMERIC(0);
1102         Numeric         num2 = PG_GETARG_NUMERIC(1);
1103         NumericVar      arg1;
1104         NumericVar      arg2;
1105         NumericVar      result;
1106         Numeric         res;
1107
1108         /*
1109          * Handle NaN
1110          */
1111         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1112                 PG_RETURN_NUMERIC(make_result(&const_nan));
1113
1114         /*
1115          * Unpack the values, let add_var() compute the result and return it.
1116          */
1117         init_var(&arg1);
1118         init_var(&arg2);
1119         init_var(&result);
1120
1121         set_var_from_num(num1, &arg1);
1122         set_var_from_num(num2, &arg2);
1123
1124         add_var(&arg1, &arg2, &result);
1125
1126         res = make_result(&result);
1127
1128         free_var(&arg1);
1129         free_var(&arg2);
1130         free_var(&result);
1131
1132         PG_RETURN_NUMERIC(res);
1133 }
1134
1135
1136 /*
1137  * numeric_sub() -
1138  *
1139  *      Subtract one numeric from another
1140  */
1141 Datum
1142 numeric_sub(PG_FUNCTION_ARGS)
1143 {
1144         Numeric         num1 = PG_GETARG_NUMERIC(0);
1145         Numeric         num2 = PG_GETARG_NUMERIC(1);
1146         NumericVar      arg1;
1147         NumericVar      arg2;
1148         NumericVar      result;
1149         Numeric         res;
1150
1151         /*
1152          * Handle NaN
1153          */
1154         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1155                 PG_RETURN_NUMERIC(make_result(&const_nan));
1156
1157         /*
1158          * Unpack the values, let sub_var() compute the result and return it.
1159          */
1160         init_var(&arg1);
1161         init_var(&arg2);
1162         init_var(&result);
1163
1164         set_var_from_num(num1, &arg1);
1165         set_var_from_num(num2, &arg2);
1166
1167         sub_var(&arg1, &arg2, &result);
1168
1169         res = make_result(&result);
1170
1171         free_var(&arg1);
1172         free_var(&arg2);
1173         free_var(&result);
1174
1175         PG_RETURN_NUMERIC(res);
1176 }
1177
1178
1179 /*
1180  * numeric_mul() -
1181  *
1182  *      Calculate the product of two numerics
1183  */
1184 Datum
1185 numeric_mul(PG_FUNCTION_ARGS)
1186 {
1187         Numeric         num1 = PG_GETARG_NUMERIC(0);
1188         Numeric         num2 = PG_GETARG_NUMERIC(1);
1189         NumericVar      arg1;
1190         NumericVar      arg2;
1191         NumericVar      result;
1192         Numeric         res;
1193
1194         /*
1195          * Handle NaN
1196          */
1197         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1198                 PG_RETURN_NUMERIC(make_result(&const_nan));
1199
1200         /*
1201          * Unpack the values, let mul_var() compute the result and return it.
1202          * Unlike add_var() and sub_var(), mul_var() will round its result. In the
1203          * case of numeric_mul(), which is invoked for the * operator on numerics,
1204          * we request exact representation for the product (rscale = sum(dscale of
1205          * arg1, dscale of arg2)).
1206          */
1207         init_var(&arg1);
1208         init_var(&arg2);
1209         init_var(&result);
1210
1211         set_var_from_num(num1, &arg1);
1212         set_var_from_num(num2, &arg2);
1213
1214         mul_var(&arg1, &arg2, &result, arg1.dscale + arg2.dscale);
1215
1216         res = make_result(&result);
1217
1218         free_var(&arg1);
1219         free_var(&arg2);
1220         free_var(&result);
1221
1222         PG_RETURN_NUMERIC(res);
1223 }
1224
1225
1226 /*
1227  * numeric_div() -
1228  *
1229  *      Divide one numeric into another
1230  */
1231 Datum
1232 numeric_div(PG_FUNCTION_ARGS)
1233 {
1234         Numeric         num1 = PG_GETARG_NUMERIC(0);
1235         Numeric         num2 = PG_GETARG_NUMERIC(1);
1236         NumericVar      arg1;
1237         NumericVar      arg2;
1238         NumericVar      result;
1239         Numeric         res;
1240         int                     rscale;
1241
1242         /*
1243          * Handle NaN
1244          */
1245         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1246                 PG_RETURN_NUMERIC(make_result(&const_nan));
1247
1248         /*
1249          * Unpack the arguments
1250          */
1251         init_var(&arg1);
1252         init_var(&arg2);
1253         init_var(&result);
1254
1255         set_var_from_num(num1, &arg1);
1256         set_var_from_num(num2, &arg2);
1257
1258         /*
1259          * Select scale for division result
1260          */
1261         rscale = select_div_scale(&arg1, &arg2);
1262
1263         /*
1264          * Do the divide and return the result
1265          */
1266         div_var(&arg1, &arg2, &result, rscale, true);
1267
1268         res = make_result(&result);
1269
1270         free_var(&arg1);
1271         free_var(&arg2);
1272         free_var(&result);
1273
1274         PG_RETURN_NUMERIC(res);
1275 }
1276
1277
1278 /*
1279  * numeric_mod() -
1280  *
1281  *      Calculate the modulo of two numerics
1282  */
1283 Datum
1284 numeric_mod(PG_FUNCTION_ARGS)
1285 {
1286         Numeric         num1 = PG_GETARG_NUMERIC(0);
1287         Numeric         num2 = PG_GETARG_NUMERIC(1);
1288         Numeric         res;
1289         NumericVar      arg1;
1290         NumericVar      arg2;
1291         NumericVar      result;
1292
1293         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1294                 PG_RETURN_NUMERIC(make_result(&const_nan));
1295
1296         init_var(&arg1);
1297         init_var(&arg2);
1298         init_var(&result);
1299
1300         set_var_from_num(num1, &arg1);
1301         set_var_from_num(num2, &arg2);
1302
1303         mod_var(&arg1, &arg2, &result);
1304
1305         res = make_result(&result);
1306
1307         free_var(&result);
1308         free_var(&arg2);
1309         free_var(&arg1);
1310
1311         PG_RETURN_NUMERIC(res);
1312 }
1313
1314
1315 /*
1316  * numeric_inc() -
1317  *
1318  *      Increment a number by one
1319  */
1320 Datum
1321 numeric_inc(PG_FUNCTION_ARGS)
1322 {
1323         Numeric         num = PG_GETARG_NUMERIC(0);
1324         NumericVar      arg;
1325         Numeric         res;
1326
1327         /*
1328          * Handle NaN
1329          */
1330         if (NUMERIC_IS_NAN(num))
1331                 PG_RETURN_NUMERIC(make_result(&const_nan));
1332
1333         /*
1334          * Compute the result and return it
1335          */
1336         init_var(&arg);
1337
1338         set_var_from_num(num, &arg);
1339
1340         add_var(&arg, &const_one, &arg);
1341
1342         res = make_result(&arg);
1343
1344         free_var(&arg);
1345
1346         PG_RETURN_NUMERIC(res);
1347 }
1348
1349
1350 /*
1351  * numeric_smaller() -
1352  *
1353  *      Return the smaller of two numbers
1354  */
1355 Datum
1356 numeric_smaller(PG_FUNCTION_ARGS)
1357 {
1358         Numeric         num1 = PG_GETARG_NUMERIC(0);
1359         Numeric         num2 = PG_GETARG_NUMERIC(1);
1360
1361         /*
1362          * Use cmp_numerics so that this will agree with the comparison operators,
1363          * particularly as regards comparisons involving NaN.
1364          */
1365         if (cmp_numerics(num1, num2) < 0)
1366                 PG_RETURN_NUMERIC(num1);
1367         else
1368                 PG_RETURN_NUMERIC(num2);
1369 }
1370
1371
1372 /*
1373  * numeric_larger() -
1374  *
1375  *      Return the larger of two numbers
1376  */
1377 Datum
1378 numeric_larger(PG_FUNCTION_ARGS)
1379 {
1380         Numeric         num1 = PG_GETARG_NUMERIC(0);
1381         Numeric         num2 = PG_GETARG_NUMERIC(1);
1382
1383         /*
1384          * Use cmp_numerics so that this will agree with the comparison operators,
1385          * particularly as regards comparisons involving NaN.
1386          */
1387         if (cmp_numerics(num1, num2) > 0)
1388                 PG_RETURN_NUMERIC(num1);
1389         else
1390                 PG_RETURN_NUMERIC(num2);
1391 }
1392
1393
1394 /* ----------------------------------------------------------------------
1395  *
1396  * Advanced math functions
1397  *
1398  * ----------------------------------------------------------------------
1399  */
1400
1401 /*
1402  * numeric_fac()
1403  *
1404  * Compute factorial
1405  */
1406 Datum
1407 numeric_fac(PG_FUNCTION_ARGS)
1408 {
1409         int64           num = PG_GETARG_INT64(0);
1410         Numeric         res;
1411         NumericVar      fact;
1412         NumericVar      result;
1413
1414         if (num <= 1)
1415         {
1416                 res = make_result(&const_one);
1417                 PG_RETURN_NUMERIC(res);
1418         }
1419
1420         init_var(&fact);
1421         init_var(&result);
1422
1423         int8_to_numericvar(num, &result);
1424
1425         for (num = num - 1; num > 1; num--)
1426         {
1427                 int8_to_numericvar(num, &fact);
1428
1429                 mul_var(&result, &fact, &result, 0);
1430         }
1431
1432         res = make_result(&result);
1433
1434         free_var(&fact);
1435         free_var(&result);
1436
1437         PG_RETURN_NUMERIC(res);
1438 }
1439
1440
1441 /*
1442  * numeric_sqrt() -
1443  *
1444  *      Compute the square root of a numeric.
1445  */
1446 Datum
1447 numeric_sqrt(PG_FUNCTION_ARGS)
1448 {
1449         Numeric         num = PG_GETARG_NUMERIC(0);
1450         Numeric         res;
1451         NumericVar      arg;
1452         NumericVar      result;
1453         int                     sweight;
1454         int                     rscale;
1455
1456         /*
1457          * Handle NaN
1458          */
1459         if (NUMERIC_IS_NAN(num))
1460                 PG_RETURN_NUMERIC(make_result(&const_nan));
1461
1462         /*
1463          * Unpack the argument and determine the result scale.  We choose a scale
1464          * to give at least NUMERIC_MIN_SIG_DIGITS significant digits; but in any
1465          * case not less than the input's dscale.
1466          */
1467         init_var(&arg);
1468         init_var(&result);
1469
1470         set_var_from_num(num, &arg);
1471
1472         /* Assume the input was normalized, so arg.weight is accurate */
1473         sweight = (arg.weight + 1) * DEC_DIGITS / 2 - 1;
1474
1475         rscale = NUMERIC_MIN_SIG_DIGITS - sweight;
1476         rscale = Max(rscale, arg.dscale);
1477         rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
1478         rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
1479
1480         /*
1481          * Let sqrt_var() do the calculation and return the result.
1482          */
1483         sqrt_var(&arg, &result, rscale);
1484
1485         res = make_result(&result);
1486
1487         free_var(&result);
1488         free_var(&arg);
1489
1490         PG_RETURN_NUMERIC(res);
1491 }
1492
1493
1494 /*
1495  * numeric_exp() -
1496  *
1497  *      Raise e to the power of x
1498  */
1499 Datum
1500 numeric_exp(PG_FUNCTION_ARGS)
1501 {
1502         Numeric         num = PG_GETARG_NUMERIC(0);
1503         Numeric         res;
1504         NumericVar      arg;
1505         NumericVar      result;
1506         int                     rscale;
1507         double          val;
1508
1509         /*
1510          * Handle NaN
1511          */
1512         if (NUMERIC_IS_NAN(num))
1513                 PG_RETURN_NUMERIC(make_result(&const_nan));
1514
1515         /*
1516          * Unpack the argument and determine the result scale.  We choose a scale
1517          * to give at least NUMERIC_MIN_SIG_DIGITS significant digits; but in any
1518          * case not less than the input's dscale.
1519          */
1520         init_var(&arg);
1521         init_var(&result);
1522
1523         set_var_from_num(num, &arg);
1524
1525         /* convert input to float8, ignoring overflow */
1526         val = numericvar_to_double_no_overflow(&arg);
1527
1528         /*
1529          * log10(result) = num * log10(e), so this is approximately the decimal
1530          * weight of the result:
1531          */
1532         val *= 0.434294481903252;
1533
1534         /* limit to something that won't cause integer overflow */
1535         val = Max(val, -NUMERIC_MAX_RESULT_SCALE);
1536         val = Min(val, NUMERIC_MAX_RESULT_SCALE);
1537
1538         rscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
1539         rscale = Max(rscale, arg.dscale);
1540         rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
1541         rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
1542
1543         /*
1544          * Let exp_var() do the calculation and return the result.
1545          */
1546         exp_var(&arg, &result, rscale);
1547
1548         res = make_result(&result);
1549
1550         free_var(&result);
1551         free_var(&arg);
1552
1553         PG_RETURN_NUMERIC(res);
1554 }
1555
1556
1557 /*
1558  * numeric_ln() -
1559  *
1560  *      Compute the natural logarithm of x
1561  */
1562 Datum
1563 numeric_ln(PG_FUNCTION_ARGS)
1564 {
1565         Numeric         num = PG_GETARG_NUMERIC(0);
1566         Numeric         res;
1567         NumericVar      arg;
1568         NumericVar      result;
1569         int                     dec_digits;
1570         int                     rscale;
1571
1572         /*
1573          * Handle NaN
1574          */
1575         if (NUMERIC_IS_NAN(num))
1576                 PG_RETURN_NUMERIC(make_result(&const_nan));
1577
1578         init_var(&arg);
1579         init_var(&result);
1580
1581         set_var_from_num(num, &arg);
1582
1583         /* Approx decimal digits before decimal point */
1584         dec_digits = (arg.weight + 1) * DEC_DIGITS;
1585
1586         if (dec_digits > 1)
1587                 rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(dec_digits - 1);
1588         else if (dec_digits < 1)
1589                 rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(1 - dec_digits);
1590         else
1591                 rscale = NUMERIC_MIN_SIG_DIGITS;
1592
1593         rscale = Max(rscale, arg.dscale);
1594         rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
1595         rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
1596
1597         ln_var(&arg, &result, rscale);
1598
1599         res = make_result(&result);
1600
1601         free_var(&result);
1602         free_var(&arg);
1603
1604         PG_RETURN_NUMERIC(res);
1605 }
1606
1607
1608 /*
1609  * numeric_log() -
1610  *
1611  *      Compute the logarithm of x in a given base
1612  */
1613 Datum
1614 numeric_log(PG_FUNCTION_ARGS)
1615 {
1616         Numeric         num1 = PG_GETARG_NUMERIC(0);
1617         Numeric         num2 = PG_GETARG_NUMERIC(1);
1618         Numeric         res;
1619         NumericVar      arg1;
1620         NumericVar      arg2;
1621         NumericVar      result;
1622
1623         /*
1624          * Handle NaN
1625          */
1626         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1627                 PG_RETURN_NUMERIC(make_result(&const_nan));
1628
1629         /*
1630          * Initialize things
1631          */
1632         init_var(&arg1);
1633         init_var(&arg2);
1634         init_var(&result);
1635
1636         set_var_from_num(num1, &arg1);
1637         set_var_from_num(num2, &arg2);
1638
1639         /*
1640          * Call log_var() to compute and return the result; note it handles scale
1641          * selection itself.
1642          */
1643         log_var(&arg1, &arg2, &result);
1644
1645         res = make_result(&result);
1646
1647         free_var(&result);
1648         free_var(&arg2);
1649         free_var(&arg1);
1650
1651         PG_RETURN_NUMERIC(res);
1652 }
1653
1654
1655 /*
1656  * numeric_power() -
1657  *
1658  *      Raise b to the power of x
1659  */
1660 Datum
1661 numeric_power(PG_FUNCTION_ARGS)
1662 {
1663         Numeric         num1 = PG_GETARG_NUMERIC(0);
1664         Numeric         num2 = PG_GETARG_NUMERIC(1);
1665         Numeric         res;
1666         NumericVar      arg1;
1667         NumericVar      arg2;
1668         NumericVar      arg2_trunc;
1669         NumericVar      result;
1670
1671         /*
1672          * Handle NaN
1673          */
1674         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1675                 PG_RETURN_NUMERIC(make_result(&const_nan));
1676
1677         /*
1678          * Initialize things
1679          */
1680         init_var(&arg1);
1681         init_var(&arg2);
1682         init_var(&arg2_trunc);
1683         init_var(&result);
1684
1685         set_var_from_num(num1, &arg1);
1686         set_var_from_num(num2, &arg2);
1687         set_var_from_var(&arg2, &arg2_trunc);
1688
1689         trunc_var(&arg2_trunc, 0);
1690
1691         /*
1692          * Return special SQLSTATE error codes for a few conditions mandated by
1693          * the standard.
1694          */
1695         if ((cmp_var(&arg1, &const_zero) == 0 &&
1696                  cmp_var(&arg2, &const_zero) < 0) ||
1697                 (cmp_var(&arg1, &const_zero) < 0 &&
1698                  cmp_var(&arg2, &arg2_trunc) != 0))
1699                 ereport(ERROR,
1700                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
1701                                  errmsg("invalid argument for power function")));
1702
1703         /*
1704          * Call power_var() to compute and return the result; note it handles
1705          * scale selection itself.
1706          */
1707         power_var(&arg1, &arg2, &result);
1708
1709         res = make_result(&result);
1710
1711         free_var(&result);
1712         free_var(&arg2);
1713         free_var(&arg2_trunc);
1714         free_var(&arg1);
1715
1716         PG_RETURN_NUMERIC(res);
1717 }
1718
1719
1720 /* ----------------------------------------------------------------------
1721  *
1722  * Type conversion functions
1723  *
1724  * ----------------------------------------------------------------------
1725  */
1726
1727
1728 Datum
1729 int4_numeric(PG_FUNCTION_ARGS)
1730 {
1731         int32           val = PG_GETARG_INT32(0);
1732         Numeric         res;
1733         NumericVar      result;
1734
1735         init_var(&result);
1736
1737         int8_to_numericvar((int64) val, &result);
1738
1739         res = make_result(&result);
1740
1741         free_var(&result);
1742
1743         PG_RETURN_NUMERIC(res);
1744 }
1745
1746
1747 Datum
1748 numeric_int4(PG_FUNCTION_ARGS)
1749 {
1750         Numeric         num = PG_GETARG_NUMERIC(0);
1751         NumericVar      x;
1752         int32           result;
1753
1754         /* XXX would it be better to return NULL? */
1755         if (NUMERIC_IS_NAN(num))
1756                 ereport(ERROR,
1757                                 (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1758                                  errmsg("cannot convert NaN to integer")));
1759
1760         /* Convert to variable format, then convert to int4 */
1761         init_var(&x);
1762         set_var_from_num(num, &x);
1763         result = numericvar_to_int4(&x);
1764         free_var(&x);
1765         PG_RETURN_INT32(result);
1766 }
1767
1768 /*
1769  * Given a NumericVar, convert it to an int32. If the NumericVar
1770  * exceeds the range of an int32, raise the appropriate error via
1771  * ereport(). The input NumericVar is *not* free'd.
1772  */
1773 static int32
1774 numericvar_to_int4(NumericVar *var)
1775 {
1776         int32           result;
1777         int64           val;
1778
1779         if (!numericvar_to_int8(var, &val))
1780                 ereport(ERROR,
1781                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1782                                  errmsg("integer out of range")));
1783
1784         /* Down-convert to int4 */
1785         result = (int32) val;
1786
1787         /* Test for overflow by reverse-conversion. */
1788         if ((int64) result != val)
1789                 ereport(ERROR,
1790                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1791                                  errmsg("integer out of range")));
1792
1793         return result;
1794 }
1795
1796 Datum
1797 int8_numeric(PG_FUNCTION_ARGS)
1798 {
1799         int64           val = PG_GETARG_INT64(0);
1800         Numeric         res;
1801         NumericVar      result;
1802
1803         init_var(&result);
1804
1805         int8_to_numericvar(val, &result);
1806
1807         res = make_result(&result);
1808
1809         free_var(&result);
1810
1811         PG_RETURN_NUMERIC(res);
1812 }
1813
1814
1815 Datum
1816 numeric_int8(PG_FUNCTION_ARGS)
1817 {
1818         Numeric         num = PG_GETARG_NUMERIC(0);
1819         NumericVar      x;
1820         int64           result;
1821
1822         /* XXX would it be better to return NULL? */
1823         if (NUMERIC_IS_NAN(num))
1824                 ereport(ERROR,
1825                                 (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1826                                  errmsg("cannot convert NaN to bigint")));
1827
1828         /* Convert to variable format and thence to int8 */
1829         init_var(&x);
1830         set_var_from_num(num, &x);
1831
1832         if (!numericvar_to_int8(&x, &result))
1833                 ereport(ERROR,
1834                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1835                                  errmsg("bigint out of range")));
1836
1837         free_var(&x);
1838
1839         PG_RETURN_INT64(result);
1840 }
1841
1842
1843 Datum
1844 int2_numeric(PG_FUNCTION_ARGS)
1845 {
1846         int16           val = PG_GETARG_INT16(0);
1847         Numeric         res;
1848         NumericVar      result;
1849
1850         init_var(&result);
1851
1852         int8_to_numericvar((int64) val, &result);
1853
1854         res = make_result(&result);
1855
1856         free_var(&result);
1857
1858         PG_RETURN_NUMERIC(res);
1859 }
1860
1861
1862 Datum
1863 numeric_int2(PG_FUNCTION_ARGS)
1864 {
1865         Numeric         num = PG_GETARG_NUMERIC(0);
1866         NumericVar      x;
1867         int64           val;
1868         int16           result;
1869
1870         /* XXX would it be better to return NULL? */
1871         if (NUMERIC_IS_NAN(num))
1872                 ereport(ERROR,
1873                                 (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1874                                  errmsg("cannot convert NaN to smallint")));
1875
1876         /* Convert to variable format and thence to int8 */
1877         init_var(&x);
1878         set_var_from_num(num, &x);
1879
1880         if (!numericvar_to_int8(&x, &val))
1881                 ereport(ERROR,
1882                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1883                                  errmsg("smallint out of range")));
1884
1885         free_var(&x);
1886
1887         /* Down-convert to int2 */
1888         result = (int16) val;
1889
1890         /* Test for overflow by reverse-conversion. */
1891         if ((int64) result != val)
1892                 ereport(ERROR,
1893                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1894                                  errmsg("smallint out of range")));
1895
1896         PG_RETURN_INT16(result);
1897 }
1898
1899
1900 Datum
1901 float8_numeric(PG_FUNCTION_ARGS)
1902 {
1903         float8          val = PG_GETARG_FLOAT8(0);
1904         Numeric         res;
1905         NumericVar      result;
1906         char            buf[DBL_DIG + 100];
1907
1908         if (isnan(val))
1909                 PG_RETURN_NUMERIC(make_result(&const_nan));
1910
1911         sprintf(buf, "%.*g", DBL_DIG, val);
1912
1913         init_var(&result);
1914
1915         set_var_from_str(buf, &result);
1916         res = make_result(&result);
1917
1918         free_var(&result);
1919
1920         PG_RETURN_NUMERIC(res);
1921 }
1922
1923
1924 Datum
1925 numeric_float8(PG_FUNCTION_ARGS)
1926 {
1927         Numeric         num = PG_GETARG_NUMERIC(0);
1928         char       *tmp;
1929         Datum           result;
1930
1931         if (NUMERIC_IS_NAN(num))
1932                 PG_RETURN_FLOAT8(get_float8_nan());
1933
1934         tmp = DatumGetCString(DirectFunctionCall1(numeric_out,
1935                                                                                           NumericGetDatum(num)));
1936
1937         result = DirectFunctionCall1(float8in, CStringGetDatum(tmp));
1938
1939         pfree(tmp);
1940
1941         PG_RETURN_DATUM(result);
1942 }
1943
1944
1945 /* Convert numeric to float8; if out of range, return +/- HUGE_VAL */
1946 Datum
1947 numeric_float8_no_overflow(PG_FUNCTION_ARGS)
1948 {
1949         Numeric         num = PG_GETARG_NUMERIC(0);
1950         double          val;
1951
1952         if (NUMERIC_IS_NAN(num))
1953                 PG_RETURN_FLOAT8(get_float8_nan());
1954
1955         val = numeric_to_double_no_overflow(num);
1956
1957         PG_RETURN_FLOAT8(val);
1958 }
1959
1960 Datum
1961 float4_numeric(PG_FUNCTION_ARGS)
1962 {
1963         float4          val = PG_GETARG_FLOAT4(0);
1964         Numeric         res;
1965         NumericVar      result;
1966         char            buf[FLT_DIG + 100];
1967
1968         if (isnan(val))
1969                 PG_RETURN_NUMERIC(make_result(&const_nan));
1970
1971         sprintf(buf, "%.*g", FLT_DIG, val);
1972
1973         init_var(&result);
1974
1975         set_var_from_str(buf, &result);
1976         res = make_result(&result);
1977
1978         free_var(&result);
1979
1980         PG_RETURN_NUMERIC(res);
1981 }
1982
1983
1984 Datum
1985 numeric_float4(PG_FUNCTION_ARGS)
1986 {
1987         Numeric         num = PG_GETARG_NUMERIC(0);
1988         char       *tmp;
1989         Datum           result;
1990
1991         if (NUMERIC_IS_NAN(num))
1992                 PG_RETURN_FLOAT4(get_float4_nan());
1993
1994         tmp = DatumGetCString(DirectFunctionCall1(numeric_out,
1995                                                                                           NumericGetDatum(num)));
1996
1997         result = DirectFunctionCall1(float4in, CStringGetDatum(tmp));
1998
1999         pfree(tmp);
2000
2001         PG_RETURN_DATUM(result);
2002 }
2003
2004
2005 Datum
2006 text_numeric(PG_FUNCTION_ARGS)
2007 {
2008         text       *str = PG_GETARG_TEXT_P(0);
2009         int                     len;
2010         char       *s;
2011         Datum           result;
2012
2013         len = (VARSIZE(str) - VARHDRSZ);
2014         s = palloc(len + 1);
2015         memcpy(s, VARDATA(str), len);
2016         *(s + len) = '\0';
2017
2018         result = DirectFunctionCall3(numeric_in, CStringGetDatum(s),
2019                                                                  ObjectIdGetDatum(0), Int32GetDatum(-1));
2020
2021         pfree(s);
2022
2023         return result;
2024 }
2025
2026 Datum
2027 numeric_text(PG_FUNCTION_ARGS)
2028 {
2029         /* val is numeric, but easier to leave it as Datum */
2030         Datum           val = PG_GETARG_DATUM(0);
2031         char       *s;
2032         int                     len;
2033         text       *result;
2034
2035         s = DatumGetCString(DirectFunctionCall1(numeric_out, val));
2036         len = strlen(s);
2037
2038         result = (text *) palloc(VARHDRSZ + len);
2039
2040         VARATT_SIZEP(result) = len + VARHDRSZ;
2041         memcpy(VARDATA(result), s, len);
2042
2043         pfree(s);
2044
2045         PG_RETURN_TEXT_P(result);
2046 }
2047
2048
2049 /* ----------------------------------------------------------------------
2050  *
2051  * Aggregate functions
2052  *
2053  * The transition datatype for all these aggregates is a 3-element array
2054  * of Numeric, holding the values N, sum(X), sum(X*X) in that order.
2055  *
2056  * We represent N as a numeric mainly to avoid having to build a special
2057  * datatype; it's unlikely it'd overflow an int4, but ...
2058  *
2059  * ----------------------------------------------------------------------
2060  */
2061
2062 static ArrayType *
2063 do_numeric_accum(ArrayType *transarray, Numeric newval)
2064 {
2065         Datum      *transdatums;
2066         int                     ndatums;
2067         Datum           N,
2068                                 sumX,
2069                                 sumX2;
2070         ArrayType  *result;
2071
2072         /* We assume the input is array of numeric */
2073         deconstruct_array(transarray,
2074                                           NUMERICOID, -1, false, 'i',
2075                                           &transdatums, NULL, &ndatums);
2076         if (ndatums != 3)
2077                 elog(ERROR, "expected 3-element numeric array");
2078         N = transdatums[0];
2079         sumX = transdatums[1];
2080         sumX2 = transdatums[2];
2081
2082         N = DirectFunctionCall1(numeric_inc, N);
2083         sumX = DirectFunctionCall2(numeric_add, sumX,
2084                                                            NumericGetDatum(newval));
2085         sumX2 = DirectFunctionCall2(numeric_add, sumX2,
2086                                                                 DirectFunctionCall2(numeric_mul,
2087                                                                                                         NumericGetDatum(newval),
2088                                                                                                         NumericGetDatum(newval)));
2089
2090         transdatums[0] = N;
2091         transdatums[1] = sumX;
2092         transdatums[2] = sumX2;
2093
2094         result = construct_array(transdatums, 3,
2095                                                          NUMERICOID, -1, false, 'i');
2096
2097         return result;
2098 }
2099
2100 Datum
2101 numeric_accum(PG_FUNCTION_ARGS)
2102 {
2103         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2104         Numeric         newval = PG_GETARG_NUMERIC(1);
2105
2106         PG_RETURN_ARRAYTYPE_P(do_numeric_accum(transarray, newval));
2107 }
2108
2109 /*
2110  * Integer data types all use Numeric accumulators to share code and
2111  * avoid risk of overflow.      For int2 and int4 inputs, Numeric accumulation
2112  * is overkill for the N and sum(X) values, but definitely not overkill
2113  * for the sum(X*X) value.      Hence, we use int2_accum and int4_accum only
2114  * for stddev/variance --- there are faster special-purpose accumulator
2115  * routines for SUM and AVG of these datatypes.
2116  */
2117
2118 Datum
2119 int2_accum(PG_FUNCTION_ARGS)
2120 {
2121         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2122         Datum           newval2 = PG_GETARG_DATUM(1);
2123         Numeric         newval;
2124
2125         newval = DatumGetNumeric(DirectFunctionCall1(int2_numeric, newval2));
2126
2127         PG_RETURN_ARRAYTYPE_P(do_numeric_accum(transarray, newval));
2128 }
2129
2130 Datum
2131 int4_accum(PG_FUNCTION_ARGS)
2132 {
2133         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2134         Datum           newval4 = PG_GETARG_DATUM(1);
2135         Numeric         newval;
2136
2137         newval = DatumGetNumeric(DirectFunctionCall1(int4_numeric, newval4));
2138
2139         PG_RETURN_ARRAYTYPE_P(do_numeric_accum(transarray, newval));
2140 }
2141
2142 Datum
2143 int8_accum(PG_FUNCTION_ARGS)
2144 {
2145         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2146         Datum           newval8 = PG_GETARG_DATUM(1);
2147         Numeric         newval;
2148
2149         newval = DatumGetNumeric(DirectFunctionCall1(int8_numeric, newval8));
2150
2151         PG_RETURN_ARRAYTYPE_P(do_numeric_accum(transarray, newval));
2152 }
2153
2154 Datum
2155 numeric_avg(PG_FUNCTION_ARGS)
2156 {
2157         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2158         Datum      *transdatums;
2159         int                     ndatums;
2160         Numeric         N,
2161                                 sumX;
2162
2163         /* We assume the input is array of numeric */
2164         deconstruct_array(transarray,
2165                                           NUMERICOID, -1, false, 'i',
2166                                           &transdatums, NULL, &ndatums);
2167         if (ndatums != 3)
2168                 elog(ERROR, "expected 3-element numeric array");
2169         N = DatumGetNumeric(transdatums[0]);
2170         sumX = DatumGetNumeric(transdatums[1]);
2171         /* ignore sumX2 */
2172
2173         /* SQL92 defines AVG of no values to be NULL */
2174         /* N is zero iff no digits (cf. numeric_uminus) */
2175         if (N->varlen == NUMERIC_HDRSZ)
2176                 PG_RETURN_NULL();
2177
2178         PG_RETURN_DATUM(DirectFunctionCall2(numeric_div,
2179                                                                                 NumericGetDatum(sumX),
2180                                                                                 NumericGetDatum(N)));
2181 }
2182
2183 /*
2184  * Workhorse routine for the standard deviance and variance
2185  * aggregates. 'transarray' is the aggregate's transition
2186  * array. 'variance' specifies whether we should calculate the
2187  * variance or the standard deviation. 'sample' indicates whether the
2188  * caller is interested in the sample or the population
2189  * variance/stddev.
2190  *
2191  * If appropriate variance statistic is undefined for the input,
2192  * *is_null is set to true and NULL is returned.
2193  */
2194 static Numeric
2195 numeric_stddev_internal(ArrayType *transarray,
2196                                                 bool variance, bool sample,
2197                                                 bool *is_null)
2198 {
2199         Datum      *transdatums;
2200         int                     ndatums;
2201         Numeric         N,
2202                                 sumX,
2203                                 sumX2,
2204                                 res;
2205         NumericVar      vN,
2206                                 vsumX,
2207                                 vsumX2,
2208                                 vNminus1;
2209         NumericVar *comp;
2210         int                     rscale;
2211
2212         *is_null = false;
2213
2214         /* We assume the input is array of numeric */
2215         deconstruct_array(transarray,
2216                                           NUMERICOID, -1, false, 'i',
2217                                           &transdatums, NULL, &ndatums);
2218         if (ndatums != 3)
2219                 elog(ERROR, "expected 3-element numeric array");
2220         N = DatumGetNumeric(transdatums[0]);
2221         sumX = DatumGetNumeric(transdatums[1]);
2222         sumX2 = DatumGetNumeric(transdatums[2]);
2223
2224         if (NUMERIC_IS_NAN(N) || NUMERIC_IS_NAN(sumX) || NUMERIC_IS_NAN(sumX2))
2225                 return make_result(&const_nan);
2226
2227         init_var(&vN);
2228         set_var_from_num(N, &vN);
2229
2230         /*
2231          * Sample stddev and variance are undefined when N <= 1;
2232          * population stddev is undefined when N == 0. Return NULL in
2233          * either case.
2234          */
2235         if (sample)
2236                 comp = &const_one;
2237         else
2238                 comp = &const_zero;
2239
2240         if (cmp_var(&vN, comp) <= 0)
2241         {
2242                 free_var(&vN);
2243                 *is_null = true;
2244                 return NULL;
2245         }
2246
2247         init_var(&vNminus1);
2248         sub_var(&vN, &const_one, &vNminus1);
2249
2250         init_var(&vsumX);
2251         set_var_from_num(sumX, &vsumX);
2252         init_var(&vsumX2);
2253         set_var_from_num(sumX2, &vsumX2);
2254
2255         /* compute rscale for mul_var calls */
2256         rscale = vsumX.dscale * 2;
2257
2258         mul_var(&vsumX, &vsumX, &vsumX, rscale);        /* vsumX = sumX * sumX */
2259         mul_var(&vN, &vsumX2, &vsumX2, rscale);         /* vsumX2 = N * sumX2 */
2260         sub_var(&vsumX2, &vsumX, &vsumX2);                      /* N * sumX2 - sumX * sumX */
2261
2262         if (cmp_var(&vsumX2, &const_zero) <= 0)
2263         {
2264                 /* Watch out for roundoff error producing a negative numerator */
2265                 res = make_result(&const_zero);
2266         }
2267         else
2268         {
2269                 mul_var(&vN, &vNminus1, &vNminus1, 0);                          /* N * (N - 1) */
2270                 rscale = select_div_scale(&vsumX2, &vNminus1);
2271                 div_var(&vsumX2, &vNminus1, &vsumX, rscale, true);      /* variance */
2272                 if (!variance)
2273                         sqrt_var(&vsumX, &vsumX, rscale);                               /* stddev */
2274
2275                 res = make_result(&vsumX);
2276         }
2277
2278         free_var(&vN);
2279         free_var(&vNminus1);
2280         free_var(&vsumX);
2281         free_var(&vsumX2);
2282
2283         return res;
2284 }
2285
2286 Datum
2287 numeric_var_samp(PG_FUNCTION_ARGS)
2288 {
2289         Numeric res;
2290         bool is_null;
2291
2292         res = numeric_stddev_internal(PG_GETARG_ARRAYTYPE_P(0),
2293                                                                   true, true, &is_null);
2294
2295         if (is_null)
2296                 PG_RETURN_NULL();
2297         else
2298                 PG_RETURN_NUMERIC(res);
2299 }
2300
2301 Datum
2302 numeric_stddev_samp(PG_FUNCTION_ARGS)
2303 {
2304         Numeric res;
2305         bool is_null;
2306
2307         res = numeric_stddev_internal(PG_GETARG_ARRAYTYPE_P(0),
2308                                                                   false, true, &is_null);
2309
2310         if (is_null)
2311                 PG_RETURN_NULL();
2312         else
2313                 PG_RETURN_NUMERIC(res);
2314 }
2315
2316 Datum
2317 numeric_var_pop(PG_FUNCTION_ARGS)
2318 {
2319         Numeric res;
2320         bool is_null;
2321
2322         res = numeric_stddev_internal(PG_GETARG_ARRAYTYPE_P(0),
2323                                                                   true, false, &is_null);
2324
2325         if (is_null)
2326                 PG_RETURN_NULL();
2327         else
2328                 PG_RETURN_NUMERIC(res);
2329 }
2330
2331 Datum
2332 numeric_stddev_pop(PG_FUNCTION_ARGS)
2333 {
2334         Numeric res;
2335         bool is_null;
2336
2337         res = numeric_stddev_internal(PG_GETARG_ARRAYTYPE_P(0),
2338                                                                   false, false, &is_null);
2339
2340         if (is_null)
2341                 PG_RETURN_NULL();
2342         else
2343                 PG_RETURN_NUMERIC(res);
2344 }
2345
2346 /*
2347  * SUM transition functions for integer datatypes.
2348  *
2349  * To avoid overflow, we use accumulators wider than the input datatype.
2350  * A Numeric accumulator is needed for int8 input; for int4 and int2
2351  * inputs, we use int8 accumulators which should be sufficient for practical
2352  * purposes.  (The latter two therefore don't really belong in this file,
2353  * but we keep them here anyway.)
2354  *
2355  * Because SQL92 defines the SUM() of no values to be NULL, not zero,
2356  * the initial condition of the transition data value needs to be NULL. This
2357  * means we can't rely on ExecAgg to automatically insert the first non-null
2358  * data value into the transition data: it doesn't know how to do the type
2359  * conversion.  The upshot is that these routines have to be marked non-strict
2360  * and handle substitution of the first non-null input themselves.
2361  */
2362
2363 Datum
2364 int2_sum(PG_FUNCTION_ARGS)
2365 {
2366         int64           newval;
2367
2368         if (PG_ARGISNULL(0))
2369         {
2370                 /* No non-null input seen so far... */
2371                 if (PG_ARGISNULL(1))
2372                         PG_RETURN_NULL();       /* still no non-null */
2373                 /* This is the first non-null input. */
2374                 newval = (int64) PG_GETARG_INT16(1);
2375                 PG_RETURN_INT64(newval);
2376         }
2377
2378         /*
2379          * If we're invoked by nodeAgg, we can cheat and modify out first
2380          * parameter in-place to avoid palloc overhead. If not, we need to return
2381          * the new value of the transition variable.
2382          */
2383         if (fcinfo->context && IsA(fcinfo->context, AggState))
2384         {
2385                 int64      *oldsum = (int64 *) PG_GETARG_POINTER(0);
2386
2387                 /* Leave the running sum unchanged in the new input is null */
2388                 if (!PG_ARGISNULL(1))
2389                         *oldsum = *oldsum + (int64) PG_GETARG_INT16(1);
2390
2391                 PG_RETURN_POINTER(oldsum);
2392         }
2393         else
2394         {
2395                 int64           oldsum = PG_GETARG_INT64(0);
2396
2397                 /* Leave sum unchanged if new input is null. */
2398                 if (PG_ARGISNULL(1))
2399                         PG_RETURN_INT64(oldsum);
2400
2401                 /* OK to do the addition. */
2402                 newval = oldsum + (int64) PG_GETARG_INT16(1);
2403
2404                 PG_RETURN_INT64(newval);
2405         }
2406 }
2407
2408 Datum
2409 int4_sum(PG_FUNCTION_ARGS)
2410 {
2411         int64           newval;
2412
2413         if (PG_ARGISNULL(0))
2414         {
2415                 /* No non-null input seen so far... */
2416                 if (PG_ARGISNULL(1))
2417                         PG_RETURN_NULL();       /* still no non-null */
2418                 /* This is the first non-null input. */
2419                 newval = (int64) PG_GETARG_INT32(1);
2420                 PG_RETURN_INT64(newval);
2421         }
2422
2423         /*
2424          * If we're invoked by nodeAgg, we can cheat and modify out first
2425          * parameter in-place to avoid palloc overhead. If not, we need to return
2426          * the new value of the transition variable.
2427          */
2428         if (fcinfo->context && IsA(fcinfo->context, AggState))
2429         {
2430                 int64      *oldsum = (int64 *) PG_GETARG_POINTER(0);
2431
2432                 /* Leave the running sum unchanged in the new input is null */
2433                 if (!PG_ARGISNULL(1))
2434                         *oldsum = *oldsum + (int64) PG_GETARG_INT32(1);
2435
2436                 PG_RETURN_POINTER(oldsum);
2437         }
2438         else
2439         {
2440                 int64           oldsum = PG_GETARG_INT64(0);
2441
2442                 /* Leave sum unchanged if new input is null. */
2443                 if (PG_ARGISNULL(1))
2444                         PG_RETURN_INT64(oldsum);
2445
2446                 /* OK to do the addition. */
2447                 newval = oldsum + (int64) PG_GETARG_INT32(1);
2448
2449                 PG_RETURN_INT64(newval);
2450         }
2451 }
2452
2453 Datum
2454 int8_sum(PG_FUNCTION_ARGS)
2455 {
2456         Numeric         oldsum;
2457         Datum           newval;
2458
2459         if (PG_ARGISNULL(0))
2460         {
2461                 /* No non-null input seen so far... */
2462                 if (PG_ARGISNULL(1))
2463                         PG_RETURN_NULL();       /* still no non-null */
2464                 /* This is the first non-null input. */
2465                 newval = DirectFunctionCall1(int8_numeric, PG_GETARG_DATUM(1));
2466                 PG_RETURN_DATUM(newval);
2467         }
2468
2469         /*
2470          * Note that we cannot special-case the nodeAgg case here, as we do for
2471          * int2_sum and int4_sum: numeric is of variable size, so we cannot modify
2472          * our first parameter in-place.
2473          */
2474
2475         oldsum = PG_GETARG_NUMERIC(0);
2476
2477         /* Leave sum unchanged if new input is null. */
2478         if (PG_ARGISNULL(1))
2479                 PG_RETURN_NUMERIC(oldsum);
2480
2481         /* OK to do the addition. */
2482         newval = DirectFunctionCall1(int8_numeric, PG_GETARG_DATUM(1));
2483
2484         PG_RETURN_DATUM(DirectFunctionCall2(numeric_add,
2485                                                                                 NumericGetDatum(oldsum), newval));
2486 }
2487
2488
2489 /*
2490  * Routines for avg(int2) and avg(int4).  The transition datatype
2491  * is a two-element int8 array, holding count and sum.
2492  */
2493
2494 typedef struct Int8TransTypeData
2495 {
2496 #ifndef INT64_IS_BUSTED
2497         int64           count;
2498         int64           sum;
2499 #else
2500         /* "int64" isn't really 64 bits, so fake up properly-aligned fields */
2501         int32           count;
2502         int32           pad1;
2503         int32           sum;
2504         int32           pad2;
2505 #endif
2506 } Int8TransTypeData;
2507
2508 Datum
2509 int2_avg_accum(PG_FUNCTION_ARGS)
2510 {
2511         ArrayType  *transarray;
2512         int16           newval = PG_GETARG_INT16(1);
2513         Int8TransTypeData *transdata;
2514
2515         /*
2516          * If we're invoked by nodeAgg, we can cheat and modify our first
2517          * parameter in-place to reduce palloc overhead. Otherwise we need to make
2518          * a copy of it before scribbling on it.
2519          */
2520         if (fcinfo->context && IsA(fcinfo->context, AggState))
2521                 transarray = PG_GETARG_ARRAYTYPE_P(0);
2522         else
2523                 transarray = PG_GETARG_ARRAYTYPE_P_COPY(0);
2524
2525         if (ARR_HASNULL(transarray) ||
2526                 ARR_SIZE(transarray) != ARR_OVERHEAD_NONULLS(1) + sizeof(Int8TransTypeData))
2527                 elog(ERROR, "expected 2-element int8 array");
2528
2529         transdata = (Int8TransTypeData *) ARR_DATA_PTR(transarray);
2530         transdata->count++;
2531         transdata->sum += newval;
2532
2533         PG_RETURN_ARRAYTYPE_P(transarray);
2534 }
2535
2536 Datum
2537 int4_avg_accum(PG_FUNCTION_ARGS)
2538 {
2539         ArrayType  *transarray;
2540         int32           newval = PG_GETARG_INT32(1);
2541         Int8TransTypeData *transdata;
2542
2543         /*
2544          * If we're invoked by nodeAgg, we can cheat and modify our first
2545          * parameter in-place to reduce palloc overhead. Otherwise we need to make
2546          * a copy of it before scribbling on it.
2547          */
2548         if (fcinfo->context && IsA(fcinfo->context, AggState))
2549                 transarray = PG_GETARG_ARRAYTYPE_P(0);
2550         else
2551                 transarray = PG_GETARG_ARRAYTYPE_P_COPY(0);
2552
2553         if (ARR_HASNULL(transarray) ||
2554                 ARR_SIZE(transarray) != ARR_OVERHEAD_NONULLS(1) + sizeof(Int8TransTypeData))
2555                 elog(ERROR, "expected 2-element int8 array");
2556
2557         transdata = (Int8TransTypeData *) ARR_DATA_PTR(transarray);
2558         transdata->count++;
2559         transdata->sum += newval;
2560
2561         PG_RETURN_ARRAYTYPE_P(transarray);
2562 }
2563
2564 Datum
2565 int8_avg(PG_FUNCTION_ARGS)
2566 {
2567         ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
2568         Int8TransTypeData *transdata;
2569         Datum           countd,
2570                                 sumd;
2571
2572         if (ARR_HASNULL(transarray) ||
2573                 ARR_SIZE(transarray) != ARR_OVERHEAD_NONULLS(1) + sizeof(Int8TransTypeData))
2574                 elog(ERROR, "expected 2-element int8 array");
2575         transdata = (Int8TransTypeData *) ARR_DATA_PTR(transarray);
2576
2577         /* SQL92 defines AVG of no values to be NULL */
2578         if (transdata->count == 0)
2579                 PG_RETURN_NULL();
2580
2581         countd = DirectFunctionCall1(int8_numeric,
2582                                                                  Int64GetDatumFast(transdata->count));
2583         sumd = DirectFunctionCall1(int8_numeric,
2584                                                            Int64GetDatumFast(transdata->sum));
2585
2586         PG_RETURN_DATUM(DirectFunctionCall2(numeric_div, sumd, countd));
2587 }
2588
2589
2590 /* ----------------------------------------------------------------------
2591  *
2592  * Debug support
2593  *
2594  * ----------------------------------------------------------------------
2595  */
2596
2597 #ifdef NUMERIC_DEBUG
2598
2599 /*
2600  * dump_numeric() - Dump a value in the db storage format for debugging
2601  */
2602 static void
2603 dump_numeric(const char *str, Numeric num)
2604 {
2605         NumericDigit *digits = NUMERIC_DIGITS(num);
2606         int                     ndigits;
2607         int                     i;
2608
2609         ndigits = NUMERIC_NDIGITS(num);
2610
2611         printf("%s: NUMERIC w=%d d=%d ", str, num->n_weight, NUMERIC_DSCALE(num));
2612         switch (NUMERIC_SIGN(num))
2613         {
2614                 case NUMERIC_POS:
2615                         printf("POS");
2616                         break;
2617                 case NUMERIC_NEG:
2618                         printf("NEG");
2619                         break;
2620                 case NUMERIC_NAN:
2621                         printf("NaN");
2622                         break;
2623                 default:
2624                         printf("SIGN=0x%x", NUMERIC_SIGN(num));
2625                         break;
2626         }
2627
2628         for (i = 0; i < ndigits; i++)
2629                 printf(" %0*d", DEC_DIGITS, digits[i]);
2630         printf("\n");
2631 }
2632
2633
2634 /*
2635  * dump_var() - Dump a value in the variable format for debugging
2636  */
2637 static void
2638 dump_var(const char *str, NumericVar *var)
2639 {
2640         int                     i;
2641
2642         printf("%s: VAR w=%d d=%d ", str, var->weight, var->dscale);
2643         switch (var->sign)
2644         {
2645                 case NUMERIC_POS:
2646                         printf("POS");
2647                         break;
2648                 case NUMERIC_NEG:
2649                         printf("NEG");
2650                         break;
2651                 case NUMERIC_NAN:
2652                         printf("NaN");
2653                         break;
2654                 default:
2655                         printf("SIGN=0x%x", var->sign);
2656                         break;
2657         }
2658
2659         for (i = 0; i < var->ndigits; i++)
2660                 printf(" %0*d", DEC_DIGITS, var->digits[i]);
2661
2662         printf("\n");
2663 }
2664 #endif   /* NUMERIC_DEBUG */
2665
2666
2667 /* ----------------------------------------------------------------------
2668  *
2669  * Local functions follow
2670  *
2671  * In general, these do not support NaNs --- callers must eliminate
2672  * the possibility of NaN first.  (make_result() is an exception.)
2673  *
2674  * ----------------------------------------------------------------------
2675  */
2676
2677
2678 /*
2679  * alloc_var() -
2680  *
2681  *      Allocate a digit buffer of ndigits digits (plus a spare digit for rounding)
2682  */
2683 static void
2684 alloc_var(NumericVar *var, int ndigits)
2685 {
2686         digitbuf_free(var->buf);
2687         var->buf = digitbuf_alloc(ndigits + 1);
2688         var->buf[0] = 0;                        /* spare digit for rounding */
2689         var->digits = var->buf + 1;
2690         var->ndigits = ndigits;
2691 }
2692
2693
2694 /*
2695  * free_var() -
2696  *
2697  *      Return the digit buffer of a variable to the free pool
2698  */
2699 static void
2700 free_var(NumericVar *var)
2701 {
2702         digitbuf_free(var->buf);
2703         var->buf = NULL;
2704         var->digits = NULL;
2705         var->sign = NUMERIC_NAN;
2706 }
2707
2708
2709 /*
2710  * zero_var() -
2711  *
2712  *      Set a variable to ZERO.
2713  *      Note: its dscale is not touched.
2714  */
2715 static void
2716 zero_var(NumericVar *var)
2717 {
2718         digitbuf_free(var->buf);
2719         var->buf = NULL;
2720         var->digits = NULL;
2721         var->ndigits = 0;
2722         var->weight = 0;                        /* by convention; doesn't really matter */
2723         var->sign = NUMERIC_POS;        /* anything but NAN... */
2724 }
2725
2726
2727 /*
2728  * set_var_from_str()
2729  *
2730  *      Parse a string and put the number into a variable
2731  */
2732 static void
2733 set_var_from_str(const char *str, NumericVar *dest)
2734 {
2735         const char *cp = str;
2736         bool            have_dp = FALSE;
2737         int                     i;
2738         unsigned char *decdigits;
2739         int                     sign = NUMERIC_POS;
2740         int                     dweight = -1;
2741         int                     ddigits;
2742         int                     dscale = 0;
2743         int                     weight;
2744         int                     ndigits;
2745         int                     offset;
2746         NumericDigit *digits;
2747
2748         /*
2749          * We first parse the string to extract decimal digits and determine the
2750          * correct decimal weight.      Then convert to NBASE representation.
2751          */
2752
2753         /* skip leading spaces */
2754         while (*cp)
2755         {
2756                 if (!isspace((unsigned char) *cp))
2757                         break;
2758                 cp++;
2759         }
2760
2761         switch (*cp)
2762         {
2763                 case '+':
2764                         sign = NUMERIC_POS;
2765                         cp++;
2766                         break;
2767
2768                 case '-':
2769                         sign = NUMERIC_NEG;
2770                         cp++;
2771                         break;
2772         }
2773
2774         if (*cp == '.')
2775         {
2776                 have_dp = TRUE;
2777                 cp++;
2778         }
2779
2780         if (!isdigit((unsigned char) *cp))
2781                 ereport(ERROR,
2782                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2783                           errmsg("invalid input syntax for type numeric: \"%s\"", str)));
2784
2785         decdigits = (unsigned char *) palloc(strlen(cp) + DEC_DIGITS * 2);
2786
2787         /* leading padding for digit alignment later */
2788         memset(decdigits, 0, DEC_DIGITS);
2789         i = DEC_DIGITS;
2790
2791         while (*cp)
2792         {
2793                 if (isdigit((unsigned char) *cp))
2794                 {
2795                         decdigits[i++] = *cp++ - '0';
2796                         if (!have_dp)
2797                                 dweight++;
2798                         else
2799                                 dscale++;
2800                 }
2801                 else if (*cp == '.')
2802                 {
2803                         if (have_dp)
2804                                 ereport(ERROR,
2805                                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2806                                           errmsg("invalid input syntax for type numeric: \"%s\"",
2807                                                          str)));
2808                         have_dp = TRUE;
2809                         cp++;
2810                 }
2811                 else
2812                         break;
2813         }
2814
2815         ddigits = i - DEC_DIGITS;
2816         /* trailing padding for digit alignment later */
2817         memset(decdigits + i, 0, DEC_DIGITS - 1);
2818
2819         /* Handle exponent, if any */
2820         if (*cp == 'e' || *cp == 'E')
2821         {
2822                 long            exponent;
2823                 char       *endptr;
2824
2825                 cp++;
2826                 exponent = strtol(cp, &endptr, 10);
2827                 if (endptr == cp)
2828                         ereport(ERROR,
2829                                         (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2830                                          errmsg("invalid input syntax for type numeric: \"%s\"",
2831                                                         str)));
2832                 cp = endptr;
2833                 if (exponent > NUMERIC_MAX_PRECISION ||
2834                         exponent < -NUMERIC_MAX_PRECISION)
2835                         ereport(ERROR,
2836                                         (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2837                                          errmsg("invalid input syntax for type numeric: \"%s\"",
2838                                                         str)));
2839                 dweight += (int) exponent;
2840                 dscale -= (int) exponent;
2841                 if (dscale < 0)
2842                         dscale = 0;
2843         }
2844
2845         /* Should be nothing left but spaces */
2846         while (*cp)
2847         {
2848                 if (!isspace((unsigned char) *cp))
2849                         ereport(ERROR,
2850                                         (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2851                                          errmsg("invalid input syntax for type numeric: \"%s\"",
2852                                                         str)));
2853                 cp++;
2854         }
2855
2856         /*
2857          * Okay, convert pure-decimal representation to base NBASE.  First we need
2858          * to determine the converted weight and ndigits.  offset is the number of
2859          * decimal zeroes to insert before the first given digit to have a
2860          * correctly aligned first NBASE digit.
2861          */
2862         if (dweight >= 0)
2863                 weight = (dweight + 1 + DEC_DIGITS - 1) / DEC_DIGITS - 1;
2864         else
2865                 weight = -((-dweight - 1) / DEC_DIGITS + 1);
2866         offset = (weight + 1) * DEC_DIGITS - (dweight + 1);
2867         ndigits = (ddigits + offset + DEC_DIGITS - 1) / DEC_DIGITS;
2868
2869         alloc_var(dest, ndigits);
2870         dest->sign = sign;
2871         dest->weight = weight;
2872         dest->dscale = dscale;
2873
2874         i = DEC_DIGITS - offset;
2875         digits = dest->digits;
2876
2877         while (ndigits-- > 0)
2878         {
2879 #if DEC_DIGITS == 4
2880                 *digits++ = ((decdigits[i] * 10 + decdigits[i + 1]) * 10 +
2881                                          decdigits[i + 2]) * 10 + decdigits[i + 3];
2882 #elif DEC_DIGITS == 2
2883                 *digits++ = decdigits[i] * 10 + decdigits[i + 1];
2884 #elif DEC_DIGITS == 1
2885                 *digits++ = decdigits[i];
2886 #else
2887 #error unsupported NBASE
2888 #endif
2889                 i += DEC_DIGITS;
2890         }
2891
2892         pfree(decdigits);
2893
2894         /* Strip any leading/trailing zeroes, and normalize weight if zero */
2895         strip_var(dest);
2896 }
2897
2898
2899 /*
2900  * set_var_from_num() -
2901  *
2902  *      Convert the packed db format into a variable
2903  */
2904 static void
2905 set_var_from_num(Numeric num, NumericVar *dest)
2906 {
2907         int                     ndigits;
2908
2909         ndigits = NUMERIC_NDIGITS(num);
2910
2911         alloc_var(dest, ndigits);
2912
2913         dest->weight = num->n_weight;
2914         dest->sign = NUMERIC_SIGN(num);
2915         dest->dscale = NUMERIC_DSCALE(num);
2916
2917         memcpy(dest->digits, num->n_data, ndigits * sizeof(NumericDigit));
2918 }
2919
2920
2921 /*
2922  * set_var_from_var() -
2923  *
2924  *      Copy one variable into another
2925  */
2926 static void
2927 set_var_from_var(NumericVar *value, NumericVar *dest)
2928 {
2929         NumericDigit *newbuf;
2930
2931         newbuf = digitbuf_alloc(value->ndigits + 1);
2932         newbuf[0] = 0;                          /* spare digit for rounding */
2933         memcpy(newbuf + 1, value->digits, value->ndigits * sizeof(NumericDigit));
2934
2935         digitbuf_free(dest->buf);
2936
2937         memmove(dest, value, sizeof(NumericVar));
2938         dest->buf = newbuf;
2939         dest->digits = newbuf + 1;
2940 }
2941
2942
2943 /*
2944  * get_str_from_var() -
2945  *
2946  *      Convert a var to text representation (guts of numeric_out).
2947  *      CAUTION: var's contents may be modified by rounding!
2948  *      Returns a palloc'd string.
2949  */
2950 static char *
2951 get_str_from_var(NumericVar *var, int dscale)
2952 {
2953         char       *str;
2954         char       *cp;
2955         char       *endcp;
2956         int                     i;
2957         int                     d;
2958         NumericDigit dig;
2959
2960 #if DEC_DIGITS > 1
2961         NumericDigit d1;
2962 #endif
2963
2964         if (dscale < 0)
2965                 dscale = 0;
2966
2967         /*
2968          * Check if we must round up before printing the value and do so.
2969          */
2970         round_var(var, dscale);
2971
2972         /*
2973          * Allocate space for the result.
2974          *
2975          * i is set to to # of decimal digits before decimal point. dscale is the
2976          * # of decimal digits we will print after decimal point. We may generate
2977          * as many as DEC_DIGITS-1 excess digits at the end, and in addition we
2978          * need room for sign, decimal point, null terminator.
2979          */
2980         i = (var->weight + 1) * DEC_DIGITS;
2981         if (i <= 0)
2982                 i = 1;
2983
2984         str = palloc(i + dscale + DEC_DIGITS + 2);
2985         cp = str;
2986
2987         /*
2988          * Output a dash for negative values
2989          */
2990         if (var->sign == NUMERIC_NEG)
2991                 *cp++ = '-';
2992
2993         /*
2994          * Output all digits before the decimal point
2995          */
2996         if (var->weight < 0)
2997         {
2998                 d = var->weight + 1;
2999                 *cp++ = '0';
3000         }
3001         else
3002         {
3003                 for (d = 0; d <= var->weight; d++)
3004                 {
3005                         dig = (d < var->ndigits) ? var->digits[d] : 0;
3006                         /* In the first digit, suppress extra leading decimal zeroes */
3007 #if DEC_DIGITS == 4
3008                         {
3009                                 bool            putit = (d > 0);
3010
3011                                 d1 = dig / 1000;
3012                                 dig -= d1 * 1000;
3013                                 putit |= (d1 > 0);
3014                                 if (putit)
3015                                         *cp++ = d1 + '0';
3016                                 d1 = dig / 100;
3017                                 dig -= d1 * 100;
3018                                 putit |= (d1 > 0);
3019                                 if (putit)
3020                                         *cp++ = d1 + '0';
3021                                 d1 = dig / 10;
3022                                 dig -= d1 * 10;
3023                                 putit |= (d1 > 0);
3024                                 if (putit)
3025                                         *cp++ = d1 + '0';
3026                                 *cp++ = dig + '0';
3027                         }
3028 #elif DEC_DIGITS == 2
3029                         d1 = dig / 10;
3030                         dig -= d1 * 10;
3031                         if (d1 > 0 || d > 0)
3032                                 *cp++ = d1 + '0';
3033                         *cp++ = dig + '0';
3034 #elif DEC_DIGITS == 1
3035                         *cp++ = dig + '0';
3036 #else
3037 #error unsupported NBASE
3038 #endif
3039                 }
3040         }
3041
3042         /*
3043          * If requested, output a decimal point and all the digits that follow it.
3044          * We initially put out a multiple of DEC_DIGITS digits, then truncate if
3045          * needed.
3046          */
3047         if (dscale > 0)
3048         {
3049                 *cp++ = '.';
3050                 endcp = cp + dscale;
3051                 for (i = 0; i < dscale; d++, i += DEC_DIGITS)
3052                 {
3053                         dig = (d >= 0 && d < var->ndigits) ? var->digits[d] : 0;
3054 #if DEC_DIGITS == 4
3055                         d1 = dig / 1000;
3056                         dig -= d1 * 1000;
3057                         *cp++ = d1 + '0';
3058                         d1 = dig / 100;
3059                         dig -= d1 * 100;
3060                         *cp++ = d1 + '0';
3061                         d1 = dig / 10;
3062                         dig -= d1 * 10;
3063                         *cp++ = d1 + '0';
3064                         *cp++ = dig + '0';
3065 #elif DEC_DIGITS == 2
3066                         d1 = dig / 10;
3067                         dig -= d1 * 10;
3068                         *cp++ = d1 + '0';
3069                         *cp++ = dig + '0';
3070 #elif DEC_DIGITS == 1
3071                         *cp++ = dig + '0';
3072 #else
3073 #error unsupported NBASE
3074 #endif
3075                 }
3076                 cp = endcp;
3077         }
3078
3079         /*
3080          * terminate the string and return it
3081          */
3082         *cp = '\0';
3083         return str;
3084 }
3085
3086
3087 /*
3088  * make_result() -
3089  *
3090  *      Create the packed db numeric format in palloc()'d memory from
3091  *      a variable.
3092  */
3093 static Numeric
3094 make_result(NumericVar *var)
3095 {
3096         Numeric         result;
3097         NumericDigit *digits = var->digits;
3098         int                     weight = var->weight;
3099         int                     sign = var->sign;
3100         int                     n;
3101         Size            len;
3102
3103         if (sign == NUMERIC_NAN)
3104         {
3105                 result = (Numeric) palloc(NUMERIC_HDRSZ);
3106
3107                 result->varlen = NUMERIC_HDRSZ;
3108                 result->n_weight = 0;
3109                 result->n_sign_dscale = NUMERIC_NAN;
3110
3111                 dump_numeric("make_result()", result);
3112                 return result;
3113         }
3114
3115         n = var->ndigits;
3116
3117         /* truncate leading zeroes */
3118         while (n > 0 && *digits == 0)
3119         {
3120                 digits++;
3121                 weight--;
3122                 n--;
3123         }
3124         /* truncate trailing zeroes */
3125         while (n > 0 && digits[n - 1] == 0)
3126                 n--;
3127
3128         /* If zero result, force to weight=0 and positive sign */
3129         if (n == 0)
3130         {
3131                 weight = 0;
3132                 sign = NUMERIC_POS;
3133         }
3134
3135         /* Build the result */
3136         len = NUMERIC_HDRSZ + n * sizeof(NumericDigit);
3137         result = (Numeric) palloc(len);
3138         result->varlen = len;
3139         result->n_weight = weight;
3140         result->n_sign_dscale = sign | (var->dscale & NUMERIC_DSCALE_MASK);
3141
3142         memcpy(result->n_data, digits, n * sizeof(NumericDigit));
3143
3144         /* Check for overflow of int16 fields */
3145         if (result->n_weight != weight ||
3146                 NUMERIC_DSCALE(result) != var->dscale)
3147                 ereport(ERROR,
3148                                 (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
3149                                  errmsg("value overflows numeric format")));
3150
3151         dump_numeric("make_result()", result);
3152         return result;
3153 }
3154
3155
3156 /*
3157  * apply_typmod() -
3158  *
3159  *      Do bounds checking and rounding according to the attributes
3160  *      typmod field.
3161  */
3162 static void
3163 apply_typmod(NumericVar *var, int32 typmod)
3164 {
3165         int                     precision;
3166         int                     scale;
3167         int                     maxdigits;
3168         int                     ddigits;
3169         int                     i;
3170
3171         /* Do nothing if we have a default typmod (-1) */
3172         if (typmod < (int32) (VARHDRSZ))
3173                 return;
3174
3175         typmod -= VARHDRSZ;
3176         precision = (typmod >> 16) & 0xffff;
3177         scale = typmod & 0xffff;
3178         maxdigits = precision - scale;
3179
3180         /* Round to target scale (and set var->dscale) */
3181         round_var(var, scale);
3182
3183         /*
3184          * Check for overflow - note we can't do this before rounding, because
3185          * rounding could raise the weight.  Also note that the var's weight could
3186          * be inflated by leading zeroes, which will be stripped before storage
3187          * but perhaps might not have been yet. In any case, we must recognize a
3188          * true zero, whose weight doesn't mean anything.
3189          */
3190         ddigits = (var->weight + 1) * DEC_DIGITS;
3191         if (ddigits > maxdigits)
3192         {
3193                 /* Determine true weight; and check for all-zero result */
3194                 for (i = 0; i < var->ndigits; i++)
3195                 {
3196                         NumericDigit dig = var->digits[i];
3197
3198                         if (dig)
3199                         {
3200                                 /* Adjust for any high-order decimal zero digits */
3201 #if DEC_DIGITS == 4
3202                                 if (dig < 10)
3203                                         ddigits -= 3;
3204                                 else if (dig < 100)
3205                                         ddigits -= 2;
3206                                 else if (dig < 1000)
3207                                         ddigits -= 1;
3208 #elif DEC_DIGITS == 2
3209                                 if (dig < 10)
3210                                         ddigits -= 1;
3211 #elif DEC_DIGITS == 1
3212                                 /* no adjustment */
3213 #else
3214 #error unsupported NBASE
3215 #endif
3216                                 if (ddigits > maxdigits)
3217                                         ereport(ERROR,
3218                                                         (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
3219                                                          errmsg("numeric field overflow"),
3220                                                          errdetail("A field with precision %d, scale %d must round to an absolute value less than %s%d.",
3221                                                                            precision, scale,
3222                                                                            /* Display 10^0 as 1 */
3223                                                                            maxdigits ? "10^" : "",
3224                                                                            maxdigits ? maxdigits : 1
3225                                                                            )));
3226                                 break;
3227                         }
3228                         ddigits -= DEC_DIGITS;
3229                 }
3230         }
3231 }
3232
3233 /*
3234  * Convert numeric to int8, rounding if needed.
3235  *
3236  * If overflow, return FALSE (no error is raised).      Return TRUE if okay.
3237  *
3238  *      CAUTION: var's contents may be modified by rounding!
3239  */
3240 static bool
3241 numericvar_to_int8(NumericVar *var, int64 *result)
3242 {
3243         NumericDigit *digits;
3244         int                     ndigits;
3245         int                     weight;
3246         int                     i;
3247         int64           val,
3248                                 oldval;
3249         bool            neg;
3250
3251         /* Round to nearest integer */
3252         round_var(var, 0);
3253
3254         /* Check for zero input */
3255         strip_var(var);
3256         ndigits = var->ndigits;
3257         if (ndigits == 0)
3258         {
3259                 *result = 0;
3260                 return true;
3261         }
3262
3263         /*
3264          * For input like 10000000000, we must treat stripped digits as real. So
3265          * the loop assumes there are weight+1 digits before the decimal point.
3266          */
3267         weight = var->weight;
3268         Assert(weight >= 0 && ndigits <= weight + 1);
3269
3270         /* Construct the result */
3271         digits = var->digits;
3272         neg = (var->sign == NUMERIC_NEG);
3273         val = digits[0];
3274         for (i = 1; i <= weight; i++)
3275         {
3276                 oldval = val;
3277                 val *= NBASE;
3278                 if (i < ndigits)
3279                         val += digits[i];
3280
3281                 /*
3282                  * The overflow check is a bit tricky because we want to accept
3283                  * INT64_MIN, which will overflow the positive accumulator.  We can
3284                  * detect this case easily though because INT64_MIN is the only
3285                  * nonzero value for which -val == val (on a two's complement machine,
3286                  * anyway).
3287                  */
3288                 if ((val / NBASE) != oldval)    /* possible overflow? */
3289                 {
3290                         if (!neg || (-val) != val || val == 0 || oldval < 0)
3291                                 return false;
3292                 }
3293         }
3294
3295         *result = neg ? -val : val;
3296         return true;
3297 }
3298
3299 /*
3300  * Convert int8 value to numeric.
3301  */
3302 static void
3303 int8_to_numericvar(int64 val, NumericVar *var)
3304 {
3305         uint64          uval,
3306                                 newuval;
3307         NumericDigit *ptr;
3308         int                     ndigits;
3309
3310         /* int8 can require at most 19 decimal digits; add one for safety */
3311         alloc_var(var, 20 / DEC_DIGITS);
3312         if (val < 0)
3313         {
3314                 var->sign = NUMERIC_NEG;
3315                 uval = -val;
3316         }
3317         else
3318         {
3319                 var->sign = NUMERIC_POS;
3320                 uval = val;
3321         }
3322         var->dscale = 0;
3323         if (val == 0)
3324         {
3325                 var->ndigits = 0;
3326                 var->weight = 0;
3327                 return;
3328         }
3329         ptr = var->digits + var->ndigits;
3330         ndigits = 0;
3331         do
3332         {
3333                 ptr--;
3334                 ndigits++;
3335                 newuval = uval / NBASE;
3336                 *ptr = uval - newuval * NBASE;
3337                 uval = newuval;
3338         } while (uval);
3339         var->digits = ptr;
3340         var->ndigits = ndigits;
3341         var->weight = ndigits - 1;
3342 }
3343
3344 /*
3345  * Convert numeric to float8; if out of range, return +/- HUGE_VAL
3346  */
3347 static double
3348 numeric_to_double_no_overflow(Numeric num)
3349 {
3350         char       *tmp;
3351         double          val;
3352         char       *endptr;
3353
3354         tmp = DatumGetCString(DirectFunctionCall1(numeric_out,
3355                                                                                           NumericGetDatum(num)));
3356
3357         /* unlike float8in, we ignore ERANGE from strtod */
3358         val = strtod(tmp, &endptr);
3359         if (*endptr != '\0')
3360         {
3361                 /* shouldn't happen ... */
3362                 ereport(ERROR,
3363                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
3364                          errmsg("invalid input syntax for type double precision: \"%s\"",
3365                                         tmp)));
3366         }
3367
3368         pfree(tmp);
3369
3370         return val;
3371 }
3372
3373 /* As above, but work from a NumericVar */
3374 static double
3375 numericvar_to_double_no_overflow(NumericVar *var)
3376 {
3377         char       *tmp;
3378         double          val;
3379         char       *endptr;
3380
3381         tmp = get_str_from_var(var, var->dscale);
3382
3383         /* unlike float8in, we ignore ERANGE from strtod */
3384         val = strtod(tmp, &endptr);
3385         if (*endptr != '\0')
3386         {
3387                 /* shouldn't happen ... */
3388                 ereport(ERROR,
3389                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
3390                          errmsg("invalid input syntax for type double precision: \"%s\"",
3391                                         tmp)));
3392         }
3393
3394         pfree(tmp);
3395
3396         return val;
3397 }
3398
3399
3400 /*
3401  * cmp_var() -
3402  *
3403  *      Compare two values on variable level.  We assume zeroes have been
3404  *      truncated to no digits.
3405  */
3406 static int
3407 cmp_var(NumericVar *var1, NumericVar *var2)
3408 {
3409         return cmp_var_common(var1->digits, var1->ndigits,
3410                                                   var1->weight, var1->sign,
3411                                                   var2->digits, var2->ndigits,
3412                                                   var2->weight, var2->sign);
3413 }
3414
3415 /*
3416  * cmp_var_common() -
3417  *
3418  *  Main routine of cmp_var(). This function can be used by both
3419  *  NumericVar and Numeric.
3420  */
3421 static int
3422 cmp_var_common(const NumericDigit *var1digits, int var1ndigits,
3423                            int var1weight, int var1sign,
3424                            const NumericDigit *var2digits, int var2ndigits,
3425                            int var2weight, int var2sign)
3426 {
3427         if (var1ndigits == 0)
3428         {
3429                 if (var2ndigits == 0)
3430                         return 0;
3431                 if (var2sign == NUMERIC_NEG)
3432                         return 1;
3433                 return -1;
3434         }
3435         if (var2ndigits == 0)
3436         {
3437                 if (var1sign == NUMERIC_POS)
3438                         return 1;
3439                 return -1;
3440         }
3441
3442         if (var1sign == NUMERIC_POS)
3443         {
3444                 if (var2sign == NUMERIC_NEG)
3445                         return 1;
3446                 return cmp_abs_common(var1digits, var1ndigits, var1weight,
3447                                                           var2digits, var2ndigits, var2weight);
3448         }
3449
3450         if (var2sign == NUMERIC_POS)
3451                 return -1;
3452
3453         return cmp_abs_common(var2digits, var2ndigits, var2weight,
3454                                                   var1digits, var1ndigits, var1weight);
3455 }
3456
3457
3458 /*
3459  * add_var() -
3460  *
3461  *      Full version of add functionality on variable level (handling signs).
3462  *      result might point to one of the operands too without danger.
3463  */
3464 static void
3465 add_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
3466 {
3467         /*
3468          * Decide on the signs of the two variables what to do
3469          */
3470         if (var1->sign == NUMERIC_POS)
3471         {
3472                 if (var2->sign == NUMERIC_POS)
3473                 {
3474                         /*
3475                          * Both are positive result = +(ABS(var1) + ABS(var2))
3476                          */
3477                         add_abs(var1, var2, result);
3478                         result->sign = NUMERIC_POS;
3479                 }
3480                 else
3481                 {
3482                         /*
3483                          * var1 is positive, var2 is negative Must compare absolute values
3484                          */
3485                         switch (cmp_abs(var1, var2))
3486                         {
3487                                 case 0:
3488                                         /* ----------
3489                                          * ABS(var1) == ABS(var2)
3490                                          * result = ZERO
3491                                          * ----------
3492                                          */
3493                                         zero_var(result);
3494                                         result->dscale = Max(var1->dscale, var2->dscale);
3495                                         break;
3496
3497                                 case 1:
3498                                         /* ----------
3499                                          * ABS(var1) > ABS(var2)
3500                                          * result = +(ABS(var1) - ABS(var2))
3501                                          * ----------
3502                                          */
3503                                         sub_abs(var1, var2, result);
3504                                         result->sign = NUMERIC_POS;
3505                                         break;
3506
3507                                 case -1:
3508                                         /* ----------
3509                                          * ABS(var1) < ABS(var2)
3510                                          * result = -(ABS(var2) - ABS(var1))
3511                                          * ----------
3512                                          */
3513                                         sub_abs(var2, var1, result);
3514                                         result->sign = NUMERIC_NEG;
3515                                         break;
3516                         }
3517                 }
3518         }
3519         else
3520         {
3521                 if (var2->sign == NUMERIC_POS)
3522                 {
3523                         /* ----------
3524                          * var1 is negative, var2 is positive
3525                          * Must compare absolute values
3526                          * ----------
3527                          */
3528                         switch (cmp_abs(var1, var2))
3529                         {
3530                                 case 0:
3531                                         /* ----------
3532                                          * ABS(var1) == ABS(var2)
3533                                          * result = ZERO
3534                                          * ----------
3535                                          */
3536                                         zero_var(result);
3537                                         result->dscale = Max(var1->dscale, var2->dscale);
3538                                         break;
3539
3540                                 case 1:
3541                                         /* ----------
3542                                          * ABS(var1) > ABS(var2)
3543                                          * result = -(ABS(var1) - ABS(var2))
3544                                          * ----------
3545                                          */
3546                                         sub_abs(var1, var2, result);
3547                                         result->sign = NUMERIC_NEG;
3548                                         break;
3549
3550                                 case -1:
3551                                         /* ----------
3552                                          * ABS(var1) < ABS(var2)
3553                                          * result = +(ABS(var2) - ABS(var1))
3554                                          * ----------
3555                                          */
3556                                         sub_abs(var2, var1, result);
3557                                         result->sign = NUMERIC_POS;
3558                                         break;
3559                         }
3560                 }
3561                 else
3562                 {
3563                         /* ----------
3564                          * Both are negative
3565                          * result = -(ABS(var1) + ABS(var2))
3566                          * ----------
3567                          */
3568                         add_abs(var1, var2, result);
3569                         result->sign = NUMERIC_NEG;
3570                 }
3571         }
3572 }
3573
3574
3575 /*
3576  * sub_var() -
3577  *
3578  *      Full version of sub functionality on variable level (handling signs).
3579  *      result might point to one of the operands too without danger.
3580  */
3581 static void
3582 sub_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
3583 {
3584         /*
3585          * Decide on the signs of the two variables what to do
3586          */
3587         if (var1->sign == NUMERIC_POS)
3588         {
3589                 if (var2->sign == NUMERIC_NEG)
3590                 {
3591                         /* ----------
3592                          * var1 is positive, var2 is negative
3593                          * result = +(ABS(var1) + ABS(var2))
3594                          * ----------
3595                          */
3596                         add_abs(var1, var2, result);
3597                         result->sign = NUMERIC_POS;
3598                 }
3599                 else
3600                 {
3601                         /* ----------
3602                          * Both are positive
3603                          * Must compare absolute values
3604                          * ----------
3605                          */
3606                         switch (cmp_abs(var1, var2))
3607                         {
3608                                 case 0:
3609                                         /* ----------
3610                                          * ABS(var1) == ABS(var2)
3611                                          * result = ZERO
3612                                          * ----------
3613                                          */
3614                                         zero_var(result);
3615                                         result->dscale = Max(var1->dscale, var2->dscale);
3616                                         break;
3617
3618                                 case 1:
3619                                         /* ----------
3620                                          * ABS(var1) > ABS(var2)
3621                                          * result = +(ABS(var1) - ABS(var2))
3622                                          * ----------
3623                                          */
3624                                         sub_abs(var1, var2, result);
3625                                         result->sign = NUMERIC_POS;
3626                                         break;
3627
3628                                 case -1:
3629                                         /* ----------
3630                                          * ABS(var1) < ABS(var2)
3631                                          * result = -(ABS(var2) - ABS(var1))
3632                                          * ----------
3633                                          */
3634                                         sub_abs(var2, var1, result);
3635                                         result->sign = NUMERIC_NEG;
3636                                         break;
3637                         }
3638                 }
3639         }
3640         else
3641         {
3642                 if (var2->sign == NUMERIC_NEG)
3643                 {
3644                         /* ----------
3645                          * Both are negative
3646                          * Must compare absolute values
3647                          * ----------
3648                          */
3649                         switch (cmp_abs(var1, var2))
3650                         {
3651                                 case 0:
3652                                         /* ----------
3653                                          * ABS(var1) == ABS(var2)
3654                                          * result = ZERO
3655                                          * ----------
3656                                          */
3657                                         zero_var(result);
3658                                         result->dscale = Max(var1->dscale, var2->dscale);
3659                                         break;
3660
3661                                 case 1:
3662                                         /* ----------
3663                                          * ABS(var1) > ABS(var2)
3664                                          * result = -(ABS(var1) - ABS(var2))
3665                                          * ----------
3666                                          */
3667                                         sub_abs(var1, var2, result);
3668                                         result->sign = NUMERIC_NEG;
3669                                         break;
3670
3671                                 case -1:
3672                                         /* ----------
3673                                          * ABS(var1) < ABS(var2)
3674                                          * result = +(ABS(var2) - ABS(var1))
3675                                          * ----------
3676                                          */
3677                                         sub_abs(var2, var1, result);
3678                                         result->sign = NUMERIC_POS;
3679                                         break;
3680                         }
3681                 }
3682                 else
3683                 {
3684                         /* ----------
3685                          * var1 is negative, var2 is positive
3686                          * result = -(ABS(var1) + ABS(var2))
3687                          * ----------
3688                          */
3689                         add_abs(var1, var2, result);
3690                         result->sign = NUMERIC_NEG;
3691                 }
3692         }
3693 }
3694
3695
3696 /*
3697  * mul_var() -
3698  *
3699  *      Multiplication on variable level. Product of var1 * var2 is stored
3700  *      in result.      Result is rounded to no more than rscale fractional digits.
3701  */
3702 static void
3703 mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
3704                 int rscale)
3705 {
3706         int                     res_ndigits;
3707         int                     res_sign;
3708         int                     res_weight;
3709         int                     maxdigits;
3710         int                *dig;
3711         int                     carry;
3712         int                     maxdig;
3713         int                     newdig;
3714         NumericDigit *res_digits;
3715         int                     i,
3716                                 ri,
3717                                 i1,
3718                                 i2;
3719
3720         /* copy these values into local vars for speed in inner loop */
3721         int                     var1ndigits = var1->ndigits;
3722         int                     var2ndigits = var2->ndigits;
3723         NumericDigit *var1digits = var1->digits;
3724         NumericDigit *var2digits = var2->digits;
3725
3726         if (var1ndigits == 0 || var2ndigits == 0)
3727         {
3728                 /* one or both inputs is zero; so is result */
3729                 zero_var(result);
3730                 result->dscale = rscale;
3731                 return;
3732         }
3733
3734         /* Determine result sign and (maximum possible) weight */
3735         if (var1->sign == var2->sign)
3736                 res_sign = NUMERIC_POS;
3737         else
3738                 res_sign = NUMERIC_NEG;
3739         res_weight = var1->weight + var2->weight + 2;
3740
3741         /*
3742          * Determine number of result digits to compute.  If the exact result
3743          * would have more than rscale fractional digits, truncate the computation
3744          * with MUL_GUARD_DIGITS guard digits.  We do that by pretending that one
3745          * or both inputs have fewer digits than they really do.
3746          */
3747         res_ndigits = var1ndigits + var2ndigits + 1;
3748         maxdigits = res_weight + 1 + (rscale * DEC_DIGITS) + MUL_GUARD_DIGITS;
3749         if (res_ndigits > maxdigits)
3750         {
3751                 if (maxdigits < 3)
3752                 {
3753                         /* no useful precision at all in the result... */
3754                         zero_var(result);
3755                         result->dscale = rscale;
3756                         return;
3757                 }
3758                 /* force maxdigits odd so that input ndigits can be equal */
3759                 if ((maxdigits & 1) == 0)
3760                         maxdigits++;
3761                 if (var1ndigits > var2ndigits)
3762                 {
3763                         var1ndigits -= res_ndigits - maxdigits;
3764                         if (var1ndigits < var2ndigits)
3765                                 var1ndigits = var2ndigits = (var1ndigits + var2ndigits) / 2;
3766                 }
3767                 else
3768                 {
3769                         var2ndigits -= res_ndigits - maxdigits;
3770                         if (var2ndigits < var1ndigits)
3771                                 var1ndigits = var2ndigits = (var1ndigits + var2ndigits) / 2;
3772                 }
3773                 res_ndigits = maxdigits;
3774                 Assert(res_ndigits == var1ndigits + var2ndigits + 1);
3775         }
3776
3777         /*
3778          * We do the arithmetic in an array "dig[]" of signed int's.  Since
3779          * INT_MAX is noticeably larger than NBASE*NBASE, this gives us headroom
3780          * to avoid normalizing carries immediately.
3781          *
3782          * maxdig tracks the maximum possible value of any dig[] entry; when this
3783          * threatens to exceed INT_MAX, we take the time to propagate carries. To
3784          * avoid overflow in maxdig itself, it actually represents the max
3785          * possible value divided by NBASE-1.
3786          */
3787         dig = (int *) palloc0(res_ndigits * sizeof(int));
3788         maxdig = 0;
3789
3790         ri = res_ndigits - 1;
3791         for (i1 = var1ndigits - 1; i1 >= 0; ri--, i1--)
3792         {
3793                 int                     var1digit = var1digits[i1];
3794
3795                 if (var1digit == 0)
3796                         continue;
3797
3798                 /* Time to normalize? */
3799                 maxdig += var1digit;
3800                 if (maxdig > INT_MAX / (NBASE - 1))
3801                 {
3802                         /* Yes, do it */
3803                         carry = 0;
3804                         for (i = res_ndigits - 1; i >= 0; i--)
3805                         {
3806                                 newdig = dig[i] + carry;
3807                                 if (newdig >= NBASE)
3808                                 {
3809                                         carry = newdig / NBASE;
3810                                         newdig -= carry * NBASE;
3811                                 }
3812                                 else
3813                                         carry = 0;
3814                                 dig[i] = newdig;
3815                         }
3816                         Assert(carry == 0);
3817                         /* Reset maxdig to indicate new worst-case */
3818                         maxdig = 1 + var1digit;
3819                 }
3820
3821                 /* Add appropriate multiple of var2 into the accumulator */
3822                 i = ri;
3823                 for (i2 = var2ndigits - 1; i2 >= 0; i2--)
3824                         dig[i--] += var1digit * var2digits[i2];
3825         }
3826
3827         /*
3828          * Now we do a final carry propagation pass to normalize the result, which
3829          * we combine with storing the result digits into the output. Note that
3830          * this is still done at full precision w/guard digits.
3831          */
3832         alloc_var(result, res_ndigits);
3833         res_digits = result->digits;
3834         carry = 0;
3835         for (i = res_ndigits - 1; i >= 0; i--)
3836         {
3837                 newdig = dig[i] + carry;
3838                 if (newdig >= NBASE)
3839                 {
3840                         carry = newdig / NBASE;
3841                         newdig -= carry * NBASE;
3842                 }
3843                 else
3844                         carry = 0;
3845                 res_digits[i] = newdig;
3846         }
3847         Assert(carry == 0);
3848
3849         pfree(dig);
3850
3851         /*
3852          * Finally, round the result to the requested precision.
3853          */
3854         result->weight = res_weight;
3855         result->sign = res_sign;
3856
3857         /* Round to target rscale (and set result->dscale) */
3858         round_var(result, rscale);
3859
3860         /* Strip leading and trailing zeroes */
3861         strip_var(result);
3862 }
3863
3864
3865 /*
3866  * div_var() -
3867  *
3868  *      Division on variable level. Quotient of var1 / var2 is stored
3869  *      in result.      Result is rounded to no more than rscale fractional digits.
3870  */
3871 static void
3872 div_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
3873                 int rscale, bool round)
3874 {
3875         int                     div_ndigits;
3876         int                     res_sign;
3877         int                     res_weight;
3878         int                *div;
3879         int                     qdigit;
3880         int                     carry;
3881         int                     maxdiv;
3882         int                     newdig;
3883         NumericDigit *res_digits;
3884         double          fdividend,
3885                                 fdivisor,
3886                                 fdivisorinverse,
3887                                 fquotient;
3888         int                     qi;
3889         int                     i;
3890
3891         /* copy these values into local vars for speed in inner loop */
3892         int                     var1ndigits = var1->ndigits;
3893         int                     var2ndigits = var2->ndigits;
3894         NumericDigit *var1digits = var1->digits;
3895         NumericDigit *var2digits = var2->digits;
3896
3897         /*
3898          * First of all division by zero check; we must not be handed an
3899          * unnormalized divisor.
3900          */
3901         if (var2ndigits == 0 || var2digits[0] == 0)
3902                 ereport(ERROR,
3903                                 (errcode(ERRCODE_DIVISION_BY_ZERO),
3904                                  errmsg("division by zero")));
3905
3906         /*
3907          * Now result zero check
3908          */
3909         if (var1ndigits == 0)
3910         {
3911                 zero_var(result);
3912                 result->dscale = rscale;
3913                 return;
3914         }
3915
3916         /*
3917          * Determine the result sign, weight and number of digits to calculate
3918          */
3919         if (var1->sign == var2->sign)
3920                 res_sign = NUMERIC_POS;
3921         else
3922                 res_sign = NUMERIC_NEG;
3923         res_weight = var1->weight - var2->weight + 1;
3924         /* The number of accurate result digits we need to produce: */
3925         div_ndigits = res_weight + 1 + (rscale + DEC_DIGITS - 1) / DEC_DIGITS;
3926         /* Add guard digits for roundoff error */
3927         div_ndigits += DIV_GUARD_DIGITS;
3928         if (div_ndigits < DIV_GUARD_DIGITS)
3929                 div_ndigits = DIV_GUARD_DIGITS;
3930         /* Must be at least var1ndigits, too, to simplify data-loading loop */
3931         if (div_ndigits < var1ndigits)
3932                 div_ndigits = var1ndigits;
3933
3934         /*
3935          * We do the arithmetic in an array "div[]" of signed int's.  Since
3936          * INT_MAX is noticeably larger than NBASE*NBASE, this gives us headroom
3937          * to avoid normalizing carries immediately.
3938          *
3939          * We start with div[] containing one zero digit followed by the
3940          * dividend's digits (plus appended zeroes to reach the desired precision
3941          * including guard digits).  Each step of the main loop computes an
3942          * (approximate) quotient digit and stores it into div[], removing one
3943          * position of dividend space.  A final pass of carry propagation takes
3944          * care of any mistaken quotient digits.
3945          */
3946         div = (int *) palloc0((div_ndigits + 1) * sizeof(int));
3947         for (i = 0; i < var1ndigits; i++)
3948                 div[i + 1] = var1digits[i];
3949
3950         /*
3951          * We estimate each quotient digit using floating-point arithmetic, taking
3952          * the first four digits of the (current) dividend and divisor. This must
3953          * be float to avoid overflow.
3954          */
3955         fdivisor = (double) var2digits[0];
3956         for (i = 1; i < 4; i++)
3957         {
3958                 fdivisor *= NBASE;
3959                 if (i < var2ndigits)
3960                         fdivisor += (double) var2digits[i];
3961         }
3962         fdivisorinverse = 1.0 / fdivisor;
3963
3964         /*
3965          * maxdiv tracks the maximum possible absolute value of any div[] entry;
3966          * when this threatens to exceed INT_MAX, we take the time to propagate
3967          * carries.  To avoid overflow in maxdiv itself, it actually represents
3968          * the max possible abs. value divided by NBASE-1.
3969          */
3970         maxdiv = 1;
3971
3972         /*
3973          * Outer loop computes next quotient digit, which will go into div[qi]
3974          */
3975         for (qi = 0; qi < div_ndigits; qi++)
3976         {
3977                 /* Approximate the current dividend value */
3978                 fdividend = (double) div[qi];
3979                 for (i = 1; i < 4; i++)
3980                 {
3981                         fdividend *= NBASE;
3982                         if (qi + i <= div_ndigits)
3983                                 fdividend += (double) div[qi + i];
3984                 }
3985                 /* Compute the (approximate) quotient digit */
3986                 fquotient = fdividend * fdivisorinverse;
3987                 qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
3988                         (((int) fquotient) - 1);        /* truncate towards -infinity */
3989
3990                 if (qdigit != 0)
3991                 {
3992                         /* Do we need to normalize now? */
3993                         maxdiv += Abs(qdigit);
3994                         if (maxdiv > INT_MAX / (NBASE - 1))
3995                         {
3996                                 /* Yes, do it */
3997                                 carry = 0;
3998                                 for (i = div_ndigits; i > qi; i--)
3999                                 {
4000                                         newdig = div[i] + carry;
4001                                         if (newdig < 0)
4002                                         {
4003                                                 carry = -((-newdig - 1) / NBASE) - 1;
4004                                                 newdig -= carry * NBASE;
4005                                         }
4006                                         else if (newdig >= NBASE)
4007                                         {
4008                                                 carry = newdig / NBASE;
4009                                                 newdig -= carry * NBASE;
4010                                         }
4011                                         else
4012                                                 carry = 0;
4013                                         div[i] = newdig;
4014                                 }
4015                                 newdig = div[qi] + carry;
4016                                 div[qi] = newdig;
4017
4018                                 /*
4019                                  * All the div[] digits except possibly div[qi] are now in the
4020                                  * range 0..NBASE-1.
4021                                  */
4022                                 maxdiv = Abs(newdig) / (NBASE - 1);
4023                                 maxdiv = Max(maxdiv, 1);
4024
4025                                 /*
4026                                  * Recompute the quotient digit since new info may have
4027                                  * propagated into the top four dividend digits
4028                                  */
4029                                 fdividend = (double) div[qi];
4030                                 for (i = 1; i < 4; i++)
4031                                 {
4032                                         fdividend *= NBASE;
4033                                         if (qi + i <= div_ndigits)
4034                                                 fdividend += (double) div[qi + i];
4035                                 }
4036                                 /* Compute the (approximate) quotient digit */
4037                                 fquotient = fdividend * fdivisorinverse;
4038                                 qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
4039                                         (((int) fquotient) - 1);        /* truncate towards -infinity */
4040                                 maxdiv += Abs(qdigit);
4041                         }
4042
4043                         /* Subtract off the appropriate multiple of the divisor */
4044                         if (qdigit != 0)
4045                         {
4046                                 int                     istop = Min(var2ndigits, div_ndigits - qi + 1);
4047
4048                                 for (i = 0; i < istop; i++)
4049                                         div[qi + i] -= qdigit * var2digits[i];
4050                         }
4051                 }
4052
4053                 /*
4054                  * The dividend digit we are about to replace might still be nonzero.
4055                  * Fold it into the next digit position.  We don't need to worry about
4056                  * overflow here since this should nearly cancel with the subtraction
4057                  * of the divisor.
4058                  */
4059                 div[qi + 1] += div[qi] * NBASE;
4060
4061                 div[qi] = qdigit;
4062         }
4063
4064         /*
4065          * Approximate and store the last quotient digit (div[div_ndigits])
4066          */
4067         fdividend = (double) div[qi];
4068         for (i = 1; i < 4; i++)
4069                 fdividend *= NBASE;
4070         fquotient = fdividend * fdivisorinverse;
4071         qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
4072                 (((int) fquotient) - 1);        /* truncate towards -infinity */
4073         div[qi] = qdigit;
4074
4075         /*
4076          * Now we do a final carry propagation pass to normalize the result, which
4077          * we combine with storing the result digits into the output. Note that
4078          * this is still done at full precision w/guard digits.
4079          */
4080         alloc_var(result, div_ndigits + 1);
4081         res_digits = result->digits;
4082         carry = 0;
4083         for (i = div_ndigits; i >= 0; i--)
4084         {
4085                 newdig = div[i] + carry;
4086                 if (newdig < 0)
4087                 {
4088                         carry = -((-newdig - 1) / NBASE) - 1;
4089                         newdig -= carry * NBASE;
4090                 }
4091                 else if (newdig >= NBASE)
4092                 {
4093                         carry = newdig / NBASE;
4094                         newdig -= carry * NBASE;
4095                 }
4096                 else
4097                         carry = 0;
4098                 res_digits[i] = newdig;
4099         }
4100         Assert(carry == 0);
4101
4102         pfree(div);
4103
4104         /*
4105          * Finally, round the result to the requested precision.
4106          */
4107         result->weight = res_weight;
4108         result->sign = res_sign;
4109
4110         /* Round to target rscale (and set result->dscale) */
4111         if (round)
4112                 round_var(result, rscale);
4113         else
4114                 trunc_var(result, rscale);
4115
4116         /* Strip leading and trailing zeroes */
4117         strip_var(result);
4118 }
4119
4120
4121 /*
4122  * Default scale selection for division
4123  *
4124  * Returns the appropriate result scale for the division result.
4125  */
4126 static int
4127 select_div_scale(NumericVar *var1, NumericVar *var2)
4128 {
4129         int                     weight1,
4130                                 weight2,
4131                                 qweight,
4132                                 i;
4133         NumericDigit firstdigit1,
4134                                 firstdigit2;
4135         int                     rscale;
4136
4137         /*
4138          * The result scale of a division isn't specified in any SQL standard. For
4139          * PostgreSQL we select a result scale that will give at least
4140          * NUMERIC_MIN_SIG_DIGITS significant digits, so that numeric gives a
4141          * result no less accurate than float8; but use a scale not less than
4142          * either input's display scale.
4143          */
4144
4145         /* Get the actual (normalized) weight and first digit of each input */
4146
4147         weight1 = 0;                            /* values to use if var1 is zero */
4148         firstdigit1 = 0;
4149         for (i = 0; i < var1->ndigits; i++)
4150         {
4151                 firstdigit1 = var1->digits[i];
4152                 if (firstdigit1 != 0)
4153                 {
4154                         weight1 = var1->weight - i;
4155                         break;
4156                 }
4157         }
4158
4159         weight2 = 0;                            /* values to use if var2 is zero */
4160         firstdigit2 = 0;
4161         for (i = 0; i < var2->ndigits; i++)
4162         {
4163                 firstdigit2 = var2->digits[i];
4164                 if (firstdigit2 != 0)
4165                 {
4166                         weight2 = var2->weight - i;
4167                         break;
4168                 }
4169         }
4170
4171         /*
4172          * Estimate weight of quotient.  If the two first digits are equal, we
4173          * can't be sure, but assume that var1 is less than var2.
4174          */
4175         qweight = weight1 - weight2;
4176         if (firstdigit1 <= firstdigit2)
4177                 qweight--;
4178
4179         /* Select result scale */
4180         rscale = NUMERIC_MIN_SIG_DIGITS - qweight * DEC_DIGITS;
4181         rscale = Max(rscale, var1->dscale);
4182         rscale = Max(rscale, var2->dscale);
4183         rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
4184         rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
4185
4186         return rscale;
4187 }
4188
4189
4190 /*
4191  * mod_var() -
4192  *
4193  *      Calculate the modulo of two numerics at variable level
4194  */
4195 static void
4196 mod_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
4197 {
4198         NumericVar      tmp;
4199         int                     rscale;
4200
4201         init_var(&tmp);
4202
4203         /* ---------
4204          * We do this using the equation
4205          *              mod(x,y) = x - trunc(x/y)*y
4206          * We set rscale the same way numeric_div and numeric_mul do
4207          * to get the right answer from the equation.  The final result,
4208          * however, need not be displayed to more precision than the inputs.
4209          * ----------
4210          */
4211         rscale = select_div_scale(var1, var2);
4212
4213         div_var(var1, var2, &tmp, rscale, false);
4214
4215         trunc_var(&tmp, 0);
4216
4217         mul_var(var2, &tmp, &tmp, var2->dscale + tmp.dscale);
4218
4219         sub_var(var1, &tmp, result);
4220
4221         round_var(result, Max(var1->dscale, var2->dscale));
4222
4223         free_var(&tmp);
4224 }
4225
4226
4227 /*
4228  * ceil_var() -
4229  *
4230  *      Return the smallest integer greater than or equal to the argument
4231  *      on variable level
4232  */
4233 static void
4234 ceil_var(NumericVar *var, NumericVar *result)
4235 {
4236         NumericVar      tmp;
4237
4238         init_var(&tmp);
4239         set_var_from_var(var, &tmp);
4240
4241         trunc_var(&tmp, 0);
4242
4243         if (var->sign == NUMERIC_POS && cmp_var(var, &tmp) != 0)
4244                 add_var(&tmp, &const_one, &tmp);
4245
4246         set_var_from_var(&tmp, result);
4247         free_var(&tmp);
4248 }
4249
4250
4251 /*
4252  * floor_var() -
4253  *
4254  *      Return the largest integer equal to or less than the argument
4255  *      on variable level
4256  */
4257 static void
4258 floor_var(NumericVar *var, NumericVar *result)
4259 {
4260         NumericVar      tmp;
4261
4262         init_var(&tmp);
4263         set_var_from_var(var, &tmp);
4264
4265         trunc_var(&tmp, 0);
4266
4267         if (var->sign == NUMERIC_NEG && cmp_var(var, &tmp) != 0)
4268                 sub_var(&tmp, &const_one, &tmp);
4269
4270         set_var_from_var(&tmp, result);
4271         free_var(&tmp);
4272 }
4273
4274
4275 /*
4276  * sqrt_var() -
4277  *
4278  *      Compute the square root of x using Newton's algorithm
4279  */
4280 static void
4281 sqrt_var(NumericVar *arg, NumericVar *result, int rscale)
4282 {
4283         NumericVar      tmp_arg;
4284         NumericVar      tmp_val;
4285         NumericVar      last_val;
4286         int                     local_rscale;
4287         int                     stat;
4288
4289         local_rscale = rscale + 8;
4290
4291         stat = cmp_var(arg, &const_zero);
4292         if (stat == 0)
4293         {
4294                 zero_var(result);
4295                 result->dscale = rscale;
4296                 return;
4297         }
4298
4299         /*
4300          * SQL2003 defines sqrt() in terms of power, so we need to emit the right
4301          * SQLSTATE error code if the operand is negative.
4302          */
4303         if (stat < 0)
4304                 ereport(ERROR,
4305                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
4306                                  errmsg("cannot take square root of a negative number")));
4307
4308         init_var(&tmp_arg);
4309         init_var(&tmp_val);
4310         init_var(&last_val);
4311
4312         /* Copy arg in case it is the same var as result */
4313         set_var_from_var(arg, &tmp_arg);
4314
4315         /*
4316          * Initialize the result to the first guess
4317          */
4318         alloc_var(result, 1);
4319         result->digits[0] = tmp_arg.digits[0] / 2;
4320         if (result->digits[0] == 0)
4321                 result->digits[0] = 1;
4322         result->weight = tmp_arg.weight / 2;
4323         result->sign = NUMERIC_POS;
4324
4325         set_var_from_var(result, &last_val);
4326
4327         for (;;)
4328         {
4329                 div_var(&tmp_arg, result, &tmp_val, local_rscale, true);
4330
4331                 add_var(result, &tmp_val, result);
4332                 mul_var(result, &const_zero_point_five, result, local_rscale);
4333
4334                 if (cmp_var(&last_val, result) == 0)
4335                         break;
4336                 set_var_from_var(result, &last_val);
4337         }
4338
4339         free_var(&last_val);
4340         free_var(&tmp_val);
4341         free_var(&tmp_arg);
4342
4343         /* Round to requested precision */
4344         round_var(result, rscale);
4345 }
4346
4347
4348 /*
4349  * exp_var() -
4350  *
4351  *      Raise e to the power of x
4352  */
4353 static void
4354 exp_var(NumericVar *arg, NumericVar *result, int rscale)
4355 {
4356         NumericVar      x;
4357         int                     xintval;
4358         bool            xneg = FALSE;
4359         int                     local_rscale;
4360
4361         /*----------
4362          * We separate the integral and fraction parts of x, then compute
4363          *              e^x = e^xint * e^xfrac
4364          * where e = exp(1) and e^xfrac = exp(xfrac) are computed by
4365          * exp_var_internal; the limited range of inputs allows that routine
4366          * to do a good job with a simple Taylor series.  Raising e^xint is
4367          * done by repeated multiplications in power_var_int.
4368          *----------
4369          */
4370         init_var(&x);
4371
4372         set_var_from_var(arg, &x);
4373
4374         if (x.sign == NUMERIC_NEG)
4375         {
4376                 xneg = TRUE;
4377                 x.sign = NUMERIC_POS;
4378         }
4379
4380         /* Extract the integer part, remove it from x */
4381         xintval = 0;
4382         while (x.weight >= 0)
4383         {
4384                 xintval *= NBASE;
4385                 if (x.ndigits > 0)
4386                 {
4387                         xintval += x.digits[0];
4388                         x.digits++;
4389                         x.ndigits--;
4390                 }
4391                 x.weight--;
4392                 /* Guard against overflow */
4393                 if (xintval >= NUMERIC_MAX_RESULT_SCALE * 3)
4394                         ereport(ERROR,
4395                                         (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
4396                                          errmsg("argument for function \"exp\" too big")));
4397         }
4398
4399         /* Select an appropriate scale for internal calculation */
4400         local_rscale = rscale + MUL_GUARD_DIGITS * 2;
4401
4402         /* Compute e^xfrac */
4403         exp_var_internal(&x, result, local_rscale);
4404
4405         /* If there's an integer part, multiply by e^xint */
4406         if (xintval > 0)
4407         {
4408                 NumericVar      e;
4409
4410                 init_var(&e);
4411                 exp_var_internal(&const_one, &e, local_rscale);
4412                 power_var_int(&e, xintval, &e, local_rscale);
4413                 mul_var(&e, result, result, local_rscale);
4414                 free_var(&e);
4415         }
4416
4417         /* Compensate for input sign, and round to requested rscale */
4418         if (xneg)
4419                 div_var(&const_one, result, result, rscale, true);
4420         else
4421                 round_var(result, rscale);
4422
4423         free_var(&x);
4424 }
4425
4426
4427 /*
4428  * exp_var_internal() -
4429  *
4430  *      Raise e to the power of x, where 0 <= x <= 1
4431  *
4432  * NB: the result should be good to at least rscale digits, but it has
4433  * *not* been rounded off; the caller must do that if wanted.
4434  */
4435 static void
4436 exp_var_internal(NumericVar *arg, NumericVar *result, int rscale)
4437 {
4438         NumericVar      x;
4439         NumericVar      xpow;
4440         NumericVar      ifac;
4441         NumericVar      elem;
4442         NumericVar      ni;
4443         int                     ndiv2 = 0;
4444         int                     local_rscale;
4445
4446         init_var(&x);
4447         init_var(&xpow);
4448         init_var(&ifac);
4449         init_var(&elem);
4450         init_var(&ni);
4451
4452         set_var_from_var(arg, &x);
4453
4454         Assert(x.sign == NUMERIC_POS);
4455
4456         local_rscale = rscale + 8;
4457
4458         /* Reduce input into range 0 <= x <= 0.01 */
4459         while (cmp_var(&x, &const_zero_point_01) > 0)
4460         {
4461                 ndiv2++;
4462                 local_rscale++;
4463                 mul_var(&x, &const_zero_point_five, &x, x.dscale + 1);
4464         }
4465
4466         /*
4467          * Use the Taylor series
4468          *
4469          * exp(x) = 1 + x + x^2/2! + x^3/3! + ...
4470          *
4471          * Given the limited range of x, this should converge reasonably quickly.
4472          * We run the series until the terms fall below the local_rscale limit.
4473          */
4474         add_var(&const_one, &x, result);
4475         set_var_from_var(&x, &xpow);
4476         set_var_from_var(&const_one, &ifac);
4477         set_var_from_var(&const_one, &ni);
4478
4479         for (;;)
4480         {
4481                 add_var(&ni, &const_one, &ni);
4482                 mul_var(&xpow, &x, &xpow, local_rscale);
4483                 mul_var(&ifac, &ni, &ifac, 0);
4484                 div_var(&xpow, &ifac, &elem, local_rscale, true);
4485
4486                 if (elem.ndigits == 0)
4487                         break;
4488
4489                 add_var(result, &elem, result);
4490         }
4491
4492         /* Compensate for argument range reduction */
4493         while (ndiv2-- > 0)
4494                 mul_var(result, result, result, local_rscale);
4495
4496         free_var(&x);
4497         free_var(&xpow);
4498         free_var(&ifac);
4499         free_var(&elem);
4500         free_var(&ni);
4501 }
4502
4503
4504 /*
4505  * ln_var() -
4506  *
4507  *      Compute the natural log of x
4508  */
4509 static void
4510 ln_var(NumericVar *arg, NumericVar *result, int rscale)
4511 {
4512         NumericVar      x;
4513         NumericVar      xx;
4514         NumericVar      ni;
4515         NumericVar      elem;
4516         NumericVar      fact;
4517         int                     local_rscale;
4518         int                     cmp;
4519
4520         cmp = cmp_var(arg, &const_zero);
4521         if (cmp == 0)
4522                 ereport(ERROR,
4523                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
4524                                  errmsg("cannot take logarithm of zero")));
4525         else if (cmp < 0)
4526                 ereport(ERROR,
4527                                 (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
4528                                  errmsg("cannot take logarithm of a negative number")));
4529
4530         local_rscale = rscale + 8;
4531
4532         init_var(&x);
4533         init_var(&xx);
4534         init_var(&ni);
4535         init_var(&elem);
4536         init_var(&fact);
4537
4538         set_var_from_var(arg, &x);
4539         set_var_from_var(&const_two, &fact);
4540
4541         /* Reduce input into range 0.9 < x < 1.1 */
4542         while (cmp_var(&x, &const_zero_point_nine) <= 0)
4543         {
4544                 local_rscale++;
4545                 sqrt_var(&x, &x, local_rscale);
4546                 mul_var(&fact, &const_two, &fact, 0);
4547         }
4548         while (cmp_var(&x, &const_one_point_one) >= 0)
4549         {
4550                 local_rscale++;
4551                 sqrt_var(&x, &x, local_rscale);
4552                 mul_var(&fact, &const_two, &fact, 0);
4553         }
4554
4555         /*
4556          * We use the Taylor series for 0.5 * ln((1+z)/(1-z)),
4557          *
4558          * z + z^3/3 + z^5/5 + ...
4559          *
4560          * where z = (x-1)/(x+1) is in the range (approximately) -0.053 .. 0.048
4561          * due to the above range-reduction of x.
4562          *
4563          * The convergence of this is not as fast as one would like, but is
4564          * tolerable given that z is small.
4565          */
4566         sub_var(&x, &const_one, result);
4567         add_var(&x, &const_one, &elem);
4568         div_var(result, &elem, result, local_rscale, true);
4569         set_var_from_var(result, &xx);
4570         mul_var(result, result, &x, local_rscale);
4571
4572         set_var_from_var(&const_one, &ni);
4573
4574         for (;;)
4575         {
4576                 add_var(&ni, &const_two, &ni);
4577                 mul_var(&xx, &x, &xx, local_rscale);
4578                 div_var(&xx, &ni, &elem, local_rscale, true);
4579
4580                 if (elem.ndigits == 0)
4581                         break;
4582
4583                 add_var(result, &elem, result);
4584
4585                 if (elem.weight < (result->weight - local_rscale * 2 / DEC_DIGITS))
4586                         break;
4587         }
4588
4589         /* Compensate for argument range reduction, round to requested rscale */
4590         mul_var(result, &fact, result, rscale);
4591
4592         free_var(&x);
4593         free_var(&xx);
4594         free_var(&ni);
4595         free_var(&elem);
4596         free_var(&fact);
4597 }
4598
4599
4600 /*
4601  * log_var() -
4602  *
4603  *      Compute the logarithm of num in a given base.
4604  *
4605  *      Note: this routine chooses dscale of the result.
4606  */
4607 static void
4608 log_var(NumericVar *base, NumericVar *num, NumericVar *result)
4609 {
4610         NumericVar      ln_base;
4611         NumericVar      ln_num;
4612         int                     dec_digits;
4613         int                     rscale;
4614         int                     local_rscale;
4615
4616         init_var(&ln_base);
4617         init_var(&ln_num);
4618
4619         /* Set scale for ln() calculations --- compare numeric_ln() */
4620
4621         /* Approx decimal digits before decimal point */
4622         dec_digits = (num->weight + 1) * DEC_DIGITS;
4623
4624         if (dec_digits > 1)
4625                 rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(dec_digits - 1);
4626         else if (dec_digits < 1)
4627                 rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(1 - dec_digits);
4628         else
4629                 rscale = NUMERIC_MIN_SIG_DIGITS;
4630
4631         rscale = Max(rscale, base->dscale);
4632         rscale = Max(rscale, num->dscale);
4633         rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
4634         rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
4635
4636         local_rscale = rscale + 8;
4637
4638         /* Form natural logarithms */
4639         ln_var(base, &ln_base, local_rscale);
4640         ln_var(num, &ln_num, local_rscale);
4641
4642         ln_base.dscale = rscale;
4643         ln_num.dscale = rscale;
4644
4645         /* Select scale for division result */
4646         rscale = select_div_scale(&ln_num, &ln_base);
4647
4648         div_var(&ln_num, &ln_base, result, rscale, true);
4649
4650         free_var(&ln_num);
4651         free_var(&ln_base);
4652 }
4653
4654
4655 /*
4656  * power_var() -
4657  *
4658  *      Raise base to the power of exp
4659  *
4660  *      Note: this routine chooses dscale of the result.
4661  */
4662 static void
4663 power_var(NumericVar *base, NumericVar *exp, NumericVar *result)
4664 {
4665         NumericVar      ln_base;
4666         NumericVar      ln_num;
4667         int                     dec_digits;
4668         int                     rscale;
4669         int                     local_rscale;
4670         double          val;
4671
4672         /* If exp can be represented as an integer, use power_var_int */
4673         if (exp->ndigits == 0 || exp->ndigits <= exp->weight + 1)
4674         {
4675                 /* exact integer, but does it fit in int? */
4676                 NumericVar      x;
4677                 int64           expval64;
4678
4679                 /* must copy because numericvar_to_int8() scribbles on input */
4680                 init_var(&x);
4681                 set_var_from_var(exp, &x);
4682                 if (numericvar_to_int8(&x, &expval64))
4683                 {
4684                         int                     expval = (int) expval64;
4685
4686                         /* Test for overflow by reverse-conversion. */
4687                         if ((int64) expval == expval64)
4688                         {
4689                                 /* Okay, select rscale */
4690                                 rscale = NUMERIC_MIN_SIG_DIGITS;
4691                                 rscale = Max(rscale, base->dscale);
4692                                 rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
4693                                 rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
4694
4695                                 power_var_int(base, expval, result, rscale);
4696
4697                                 free_var(&x);
4698                                 return;
4699                         }
4700                 }
4701                 free_var(&x);
4702         }
4703
4704         init_var(&ln_base);
4705         init_var(&ln_num);
4706
4707         /* Set scale for ln() calculation --- need extra accuracy here */
4708
4709         /* Approx decimal digits before decimal point */
4710         dec_digits = (base->weight + 1) * DEC_DIGITS;
4711
4712         if (dec_digits > 1)
4713                 rscale = NUMERIC_MIN_SIG_DIGITS * 2 - (int) log10(dec_digits - 1);
4714         else if (dec_digits < 1)
4715                 rscale = NUMERIC_MIN_SIG_DIGITS * 2 - (int) log10(1 - dec_digits);
4716         else
4717                 rscale = NUMERIC_MIN_SIG_DIGITS * 2;
4718
4719         rscale = Max(rscale, base->dscale * 2);
4720         rscale = Max(rscale, exp->dscale * 2);
4721         rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE * 2);
4722         rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE * 2);
4723
4724         local_rscale = rscale + 8;
4725
4726         ln_var(base, &ln_base, local_rscale);
4727
4728         mul_var(&ln_base, exp, &ln_num, local_rscale);
4729
4730         /* Set scale for exp() -- compare numeric_exp() */
4731
4732         /* convert input to float8, ignoring overflow */
4733         val = numericvar_to_double_no_overflow(&ln_num);
4734
4735         /*
4736          * log10(result) = num * log10(e), so this is approximately the weight:
4737          */
4738         val *= 0.434294481903252;
4739
4740         /* limit to something that won't cause integer overflow */
4741         val = Max(val, -NUMERIC_MAX_RESULT_SCALE);
4742         val = Min(val, NUMERIC_MAX_RESULT_SCALE);
4743
4744         rscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
4745         rscale = Max(rscale, base->dscale);
4746         rscale = Max(rscale, exp->dscale);
4747         rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
4748         rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
4749
4750         exp_var(&ln_num, result, rscale);
4751
4752         free_var(&ln_num);
4753         free_var(&ln_base);
4754 }
4755
4756 /*
4757  * power_var_int() -
4758  *
4759  *      Raise base to the power of exp, where exp is an integer.
4760  */
4761 static void
4762 power_var_int(NumericVar *base, int exp, NumericVar *result, int rscale)
4763 {
4764         bool            neg;
4765         NumericVar      base_prod;
4766         int                     local_rscale;
4767
4768         /* Detect some special cases, particularly 0^0. */
4769
4770         switch (exp)
4771         {
4772                 case 0:
4773                         if (base->ndigits == 0)
4774                                 ereport(ERROR,
4775                                                 (errcode(ERRCODE_FLOATING_POINT_EXCEPTION),
4776                                                  errmsg("zero raised to zero is undefined")));
4777                         set_var_from_var(&const_one, result);
4778                         result->dscale = rscale;        /* no need to round */
4779                         return;
4780                 case 1:
4781                         set_var_from_var(base, result);
4782                         round_var(result, rscale);
4783                         return;
4784                 case -1:
4785                         div_var(&const_one, base, result, rscale, true);
4786                         return;
4787                 case 2:
4788                         mul_var(base, base, result, rscale);
4789                         return;
4790                 default:
4791                         break;
4792         }
4793
4794         /*
4795          * The general case repeatedly multiplies base according to the bit
4796          * pattern of exp.      We do the multiplications with some extra precision.
4797          */
4798         neg = (exp < 0);
4799         exp = Abs(exp);
4800
4801         local_rscale = rscale + MUL_GUARD_DIGITS * 2;
4802
4803         init_var(&base_prod);
4804         set_var_from_var(base, &base_prod);
4805
4806         if (exp & 1)
4807                 set_var_from_var(base, result);
4808         else
4809                 set_var_from_var(&const_one, result);
4810
4811         while ((exp >>= 1) > 0)
4812         {
4813                 mul_var(&base_prod, &base_prod, &base_prod, local_rscale);
4814                 if (exp & 1)
4815                         mul_var(&base_prod, result, result, local_rscale);
4816         }
4817
4818         free_var(&base_prod);
4819
4820         /* Compensate for input sign, and round to requested rscale */
4821         if (neg)
4822                 div_var(&const_one, result, result, rscale, true);
4823         else
4824                 round_var(result, rscale);
4825 }
4826
4827
4828 /* ----------------------------------------------------------------------
4829  *
4830  * Following are the lowest level functions that operate unsigned
4831  * on the variable level
4832  *
4833  * ----------------------------------------------------------------------
4834  */
4835
4836
4837 /* ----------
4838  * cmp_abs() -
4839  *
4840  *      Compare the absolute values of var1 and var2
4841  *      Returns:        -1 for ABS(var1) < ABS(var2)
4842  *                              0  for ABS(var1) == ABS(var2)
4843  *                              1  for ABS(var1) > ABS(var2)
4844  * ----------
4845  */
4846 static int
4847 cmp_abs(NumericVar *var1, NumericVar *var2)
4848 {
4849         return cmp_abs_common(var1->digits, var1->ndigits, var1->weight,
4850                                                   var2->digits, var2->ndigits, var2->weight);
4851 }
4852
4853 /* ----------
4854  * cmp_abs_common() -
4855  *
4856  *  Main routine of cmp_abs(). This function can be used by both
4857  *  NumericVar and Numeric.
4858  * ----------
4859  */
4860 static int
4861 cmp_abs_common(const NumericDigit *var1digits, int var1ndigits, int var1weight,
4862                            const NumericDigit *var2digits, int var2ndigits, int var2weight)
4863 {
4864         int                     i1 = 0;
4865         int                     i2 = 0;
4866
4867         /* Check any digits before the first common digit */
4868
4869         while (var1weight > var2weight && i1 < var1ndigits)
4870         {
4871                 if (var1digits[i1++] != 0)
4872                         return 1;
4873                 var1weight--;
4874         }
4875         while (var2weight > var1weight && i2 < var2ndigits)
4876         {
4877                 if (var2digits[i2++] != 0)
4878                         return -1;
4879                 var2weight--;
4880         }
4881
4882         /* At this point, either w1 == w2 or we've run out of digits */
4883
4884         if (var1weight == var2weight)
4885         {
4886                 while (i1 < var1ndigits && i2 < var2ndigits)
4887                 {
4888                         int                     stat = var1digits[i1++] - var2digits[i2++];
4889
4890                         if (stat)
4891                         {
4892                                 if (stat > 0)
4893                                         return 1;
4894                                 return -1;
4895                         }
4896                 }
4897         }
4898
4899         /*
4900          * At this point, we've run out of digits on one side or the other; so any
4901          * remaining nonzero digits imply that side is larger
4902          */
4903         while (i1 < var1ndigits)
4904         {
4905                 if (var1digits[i1++] != 0)
4906                         return 1;
4907         }
4908         while (i2 < var2ndigits)
4909         {
4910                 if (var2digits[i2++] != 0)
4911                         return -1;
4912         }
4913
4914         return 0;
4915 }
4916
4917
4918 /*
4919  * add_abs() -
4920  *
4921  *      Add the absolute values of two variables into result.
4922  *      result might point to one of the operands without danger.
4923  */
4924 static void
4925 add_abs(NumericVar *var1, NumericVar *var2, NumericVar *result)
4926 {
4927         NumericDigit *res_buf;
4928         NumericDigit *res_digits;
4929         int                     res_ndigits;
4930         int                     res_weight;
4931         int                     res_rscale,
4932                                 rscale1,
4933                                 rscale2;
4934         int                     res_dscale;
4935         int                     i,
4936                                 i1,
4937                                 i2;
4938         int                     carry = 0;
4939
4940         /* copy these values into local vars for speed in inner loop */
4941         int                     var1ndigits = var1->ndigits;
4942         int                     var2ndigits = var2->ndigits;
4943         NumericDigit *var1digits = var1->digits;
4944         NumericDigit *var2digits = var2->digits;
4945
4946         res_weight = Max(var1->weight, var2->weight) + 1;
4947
4948         res_dscale = Max(var1->dscale, var2->dscale);
4949
4950         /* Note: here we are figuring rscale in base-NBASE digits */
4951         rscale1 = var1->ndigits - var1->weight - 1;
4952         rscale2 = var2->ndigits - var2->weight - 1;
4953         res_rscale = Max(rscale1, rscale2);
4954
4955         res_ndigits = res_rscale + res_weight + 1;
4956         if (res_ndigits <= 0)
4957                 res_ndigits = 1;
4958
4959         res_buf = digitbuf_alloc(res_ndigits + 1);
4960         res_buf[0] = 0;                         /* spare digit for later rounding */
4961         res_digits = res_buf + 1;
4962
4963         i1 = res_rscale + var1->weight + 1;
4964         i2 = res_rscale + var2->weight + 1;
4965         for (i = res_ndigits - 1; i >= 0; i--)
4966         {
4967                 i1--;
4968                 i2--;
4969                 if (i1 >= 0 && i1 < var1ndigits)
4970                         carry += var1digits[i1];
4971                 if (i2 >= 0 && i2 < var2ndigits)
4972                         carry += var2digits[i2];
4973
4974                 if (carry >= NBASE)
4975                 {
4976                         res_digits[i] = carry - NBASE;
4977                         carry = 1;
4978                 }
4979                 else
4980                 {
4981                         res_digits[i] = carry;
4982                         carry = 0;
4983                 }
4984         }
4985
4986         Assert(carry == 0);                     /* else we failed to allow for carry out */
4987
4988         digitbuf_free(result->buf);
4989         result->ndigits = res_ndigits;
4990         result->buf = res_buf;
4991         result->digits = res_digits;
4992         result->weight = res_weight;
4993         result->dscale = res_dscale;
4994
4995         /* Remove leading/trailing zeroes */
4996         strip_var(result);
4997 }
4998
4999
5000 /*
5001  * sub_abs()
5002  *
5003  *      Subtract the absolute value of var2 from the absolute value of var1
5004  *      and store in result. result might point to one of the operands
5005  *      without danger.
5006  *
5007  *      ABS(var1) MUST BE GREATER OR EQUAL ABS(var2) !!!
5008  */
5009 static void
5010 sub_abs(NumericVar *var1, NumericVar *var2, NumericVar *result)
5011 {
5012         NumericDigit *res_buf;
5013         NumericDigit *res_digits;
5014         int                     res_ndigits;
5015         int                     res_weight;
5016         int                     res_rscale,
5017                                 rscale1,
5018                                 rscale2;
5019         int                     res_dscale;
5020         int                     i,
5021                                 i1,
5022                                 i2;
5023         int                     borrow = 0;
5024
5025         /* copy these values into local vars for speed in inner loop */
5026         int                     var1ndigits = var1->ndigits;
5027         int                     var2ndigits = var2->ndigits;
5028         NumericDigit *var1digits = var1->digits;
5029         NumericDigit *var2digits = var2->digits;
5030
5031         res_weight = var1->weight;
5032
5033         res_dscale = Max(var1->dscale, var2->dscale);
5034
5035         /* Note: here we are figuring rscale in base-NBASE digits */
5036         rscale1 = var1->ndigits - var1->weight - 1;
5037         rscale2 = var2->ndigits - var2->weight - 1;
5038         res_rscale = Max(rscale1, rscale2);
5039
5040         res_ndigits = res_rscale + res_weight + 1;
5041         if (res_ndigits <= 0)
5042                 res_ndigits = 1;
5043
5044         res_buf = digitbuf_alloc(res_ndigits + 1);
5045         res_buf[0] = 0;                         /* spare digit for later rounding */
5046         res_digits = res_buf + 1;
5047
5048         i1 = res_rscale + var1->weight + 1;
5049         i2 = res_rscale + var2->weight + 1;
5050         for (i = res_ndigits - 1; i >= 0; i--)
5051         {
5052                 i1--;
5053                 i2--;
5054                 if (i1 >= 0 && i1 < var1ndigits)
5055                         borrow += var1digits[i1];
5056                 if (i2 >= 0 && i2 < var2ndigits)
5057                         borrow -= var2digits[i2];
5058
5059                 if (borrow < 0)
5060                 {
5061                         res_digits[i] = borrow + NBASE;
5062                         borrow = -1;
5063                 }
5064                 else
5065                 {
5066                         res_digits[i] = borrow;
5067                         borrow = 0;
5068                 }
5069         }
5070
5071         Assert(borrow == 0);            /* else caller gave us var1 < var2 */
5072
5073         digitbuf_free(result->buf);
5074         result->ndigits = res_ndigits;
5075         result->buf = res_buf;
5076         result->digits = res_digits;
5077         result->weight = res_weight;
5078         result->dscale = res_dscale;
5079
5080         /* Remove leading/trailing zeroes */
5081         strip_var(result);
5082 }
5083
5084 /*
5085  * round_var
5086  *
5087  * Round the value of a variable to no more than rscale decimal digits
5088  * after the decimal point.  NOTE: we allow rscale < 0 here, implying
5089  * rounding before the decimal point.
5090  */
5091 static void
5092 round_var(NumericVar *var, int rscale)
5093 {
5094         NumericDigit *digits = var->digits;
5095         int                     di;
5096         int                     ndigits;
5097         int                     carry;
5098
5099         var->dscale = rscale;
5100
5101         /* decimal digits wanted */
5102         di = (var->weight + 1) * DEC_DIGITS + rscale;
5103
5104         /*
5105          * If di = 0, the value loses all digits, but could round up to 1 if its
5106          * first extra digit is >= 5.  If di < 0 the result must be 0.
5107          */
5108         if (di < 0)
5109         {
5110                 var->ndigits = 0;
5111                 var->weight = 0;
5112                 var->sign = NUMERIC_POS;
5113         }
5114         else
5115         {
5116                 /* NBASE digits wanted */
5117                 ndigits = (di + DEC_DIGITS - 1) / DEC_DIGITS;
5118
5119                 /* 0, or number of decimal digits to keep in last NBASE digit */
5120                 di %= DEC_DIGITS;
5121
5122                 if (ndigits < var->ndigits ||
5123                         (ndigits == var->ndigits && di > 0))
5124                 {
5125                         var->ndigits = ndigits;
5126
5127 #if DEC_DIGITS == 1
5128                         /* di must be zero */
5129                         carry = (digits[ndigits] >= HALF_NBASE) ? 1 : 0;
5130 #else
5131                         if (di == 0)
5132                                 carry = (digits[ndigits] >= HALF_NBASE) ? 1 : 0;
5133                         else
5134                         {
5135                                 /* Must round within last NBASE digit */
5136                                 int                     extra,
5137                                                         pow10;
5138
5139 #if DEC_DIGITS == 4
5140                                 pow10 = round_powers[di];
5141 #elif DEC_DIGITS == 2
5142                                 pow10 = 10;
5143 #else
5144 #error unsupported NBASE
5145 #endif
5146                                 extra = digits[--ndigits] % pow10;
5147                                 digits[ndigits] -= extra;
5148                                 carry = 0;
5149                                 if (extra >= pow10 / 2)
5150                                 {
5151                                         pow10 += digits[ndigits];
5152                                         if (pow10 >= NBASE)
5153                                         {
5154                                                 pow10 -= NBASE;
5155                                                 carry = 1;
5156                                         }
5157                                         digits[ndigits] = pow10;
5158                                 }
5159                         }
5160 #endif
5161
5162                         /* Propagate carry if needed */
5163                         while (carry)
5164                         {
5165                                 carry += digits[--ndigits];
5166                                 if (carry >= NBASE)
5167                                 {
5168                                         digits[ndigits] = carry - NBASE;
5169                                         carry = 1;
5170                                 }
5171                                 else
5172                                 {
5173                                         digits[ndigits] = carry;
5174                                         carry = 0;
5175                                 }
5176                         }
5177
5178                         if (ndigits < 0)
5179                         {
5180                                 Assert(ndigits == -1);  /* better not have added > 1 digit */
5181                                 Assert(var->digits > var->buf);
5182                                 var->digits--;
5183                                 var->ndigits++;
5184                                 var->weight++;
5185                         }
5186                 }
5187         }
5188 }
5189
5190 /*
5191  * trunc_var
5192  *
5193  * Truncate the value of a variable at rscale decimal digits after the
5194  * decimal point.  NOTE: we allow rscale < 0 here, implying
5195  * truncation before the decimal point.
5196  */
5197 static void
5198 trunc_var(NumericVar *var, int rscale)
5199 {
5200         int                     di;
5201         int                     ndigits;
5202
5203         var->dscale = rscale;
5204
5205         /* decimal digits wanted */
5206         di = (var->weight + 1) * DEC_DIGITS + rscale;
5207
5208         /*
5209          * If di <= 0, the value loses all digits.
5210          */
5211         if (di <= 0)
5212         {
5213                 var->ndigits = 0;
5214                 var->weight = 0;
5215                 var->sign = NUMERIC_POS;
5216         }
5217         else
5218         {
5219                 /* NBASE digits wanted */
5220                 ndigits = (di + DEC_DIGITS - 1) / DEC_DIGITS;
5221
5222                 if (ndigits <= var->ndigits)
5223                 {
5224                         var->ndigits = ndigits;
5225
5226 #if DEC_DIGITS == 1
5227                         /* no within-digit stuff to worry about */
5228 #else
5229                         /* 0, or number of decimal digits to keep in last NBASE digit */
5230                         di %= DEC_DIGITS;
5231
5232                         if (di > 0)
5233                         {
5234                                 /* Must truncate within last NBASE digit */
5235                                 NumericDigit *digits = var->digits;
5236                                 int                     extra,
5237                                                         pow10;
5238
5239 #if DEC_DIGITS == 4
5240                                 pow10 = round_powers[di];
5241 #elif DEC_DIGITS == 2
5242                                 pow10 = 10;
5243 #else
5244 #error unsupported NBASE
5245 #endif
5246                                 extra = digits[--ndigits] % pow10;
5247                                 digits[ndigits] -= extra;
5248                         }
5249 #endif
5250                 }
5251         }
5252 }
5253
5254 /*
5255  * strip_var
5256  *
5257  * Strip any leading and trailing zeroes from a numeric variable
5258  */
5259 static void
5260 strip_var(NumericVar *var)
5261 {
5262         NumericDigit *digits = var->digits;
5263         int                     ndigits = var->ndigits;
5264
5265         /* Strip leading zeroes */
5266         while (ndigits > 0 && *digits == 0)
5267         {
5268                 digits++;
5269                 var->weight--;
5270                 ndigits--;
5271         }
5272
5273         /* Strip trailing zeroes */
5274         while (ndigits > 0 && digits[ndigits - 1] == 0)
5275                 ndigits--;
5276
5277         /* If it's zero, normalize the sign and weight */
5278         if (ndigits == 0)
5279         {
5280                 var->sign = NUMERIC_POS;
5281                 var->weight = 0;
5282         }
5283
5284         var->digits = digits;
5285         var->ndigits = ndigits;
5286 }