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