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