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