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  *        $Header: /cvsroot/pgsql/src/backend/utils/adt/numeric.c,v 1.63 2003/07/27 04:53:07 tgl Exp $
  18  *
  19  *-------------------------------------------------------------------------
  20  */
  21
  22 #include "postgres.h"
  23
  24 #include <ctype.h>
  25 #include <float.h>
  26 #include <limits.h>
  27 #include <math.h>
  28 #include <errno.h>
  29
  30 #include "catalog/pg_type.h"
  31 #include "libpq/pqformat.h"
  32 #include "utils/array.h"
  33 #include "utils/builtins.h"
  34 #include "utils/int8.h"
  35 #include "utils/numeric.h"
  36
  37 /* ----------
  38  * Uncomment the following to enable compilation of dump_numeric()
  39  * and dump_var() and to get a dump of any result produced by make_result().
  40  * ----------
  41 #define NUMERIC_DEBUG
  42  */
  43
  44
  45 /* ----------
  46  * Local definitions
  47  * ----------
  48  */
  49 #ifndef NAN
  50 #define NAN             (0.0/0.0)
  51 #endif
  52
  53
  54 /* ----------
  55  * Local data types
  56  *
  57  * Numeric values are represented in a base-NBASE floating point format.
  58  * Each "digit" ranges from 0 to NBASE-1.  The type NumericDigit is signed
  59  * and wide enough to store a digit.  We assume that NBASE*NBASE can fit in
  60  * an int.  Although the purely calculational routines could handle any even
  61  * NBASE that's less than sqrt(INT_MAX), in practice we are only interested
  62  * in NBASE a power of ten, so that I/O conversions and decimal rounding
  63  * are easy.  Also, it's actually more efficient if NBASE is rather less than
  64  * sqrt(INT_MAX), so that there is "headroom" for mul_var and div_var to
  65  * postpone processing carries.
  66  * ----------
  67  */
  68
  69 #if 0
  70 #define NBASE           10
  71 #define HALF_NBASE      5
  72 #define DEC_DIGITS      1                       /* decimal digits per NBASE digit */
  73 #define MUL_GUARD_DIGITS        4       /* these are measured in NBASE digits */
  74 #define DIV_GUARD_DIGITS        8
  75
  76 typedef signed char NumericDigit;
  77 #endif
  78
  79 #if 0
  80 #define NBASE           100
  81 #define HALF_NBASE      50
  82 #define DEC_DIGITS      2                       /* decimal digits per NBASE digit */
  83 #define MUL_GUARD_DIGITS        3       /* these are measured in NBASE digits */
  84 #define DIV_GUARD_DIGITS        6
  85
  86 typedef signed char NumericDigit;
  87 #endif
  88
  89 #if 1
  90 #define NBASE           10000
  91 #define HALF_NBASE      5000
  92 #define DEC_DIGITS      4                       /* decimal digits per NBASE digit */
  93 #define MUL_GUARD_DIGITS        2       /* these are measured in NBASE digits */
  94 #define DIV_GUARD_DIGITS        4
  95
  96 typedef int16 NumericDigit;
  97 #endif
  98
  99
 100 /* ----------
 101  * The value represented by a NumericVar is determined by the sign, weight,
 102  * ndigits, and digits[] array.
 103  * Note: the first digit of a NumericVar's value is assumed to be multiplied
 104  * by NBASE ** weight.  Another way to say it is that there are weight+1
 105  * digits before the decimal point.  It is possible to have weight < 0.
 106  *
 107  * buf points at the physical start of the palloc'd digit buffer for the
 108  * NumericVar.  digits points at the first digit in actual use (the one
 109  * with the specified weight).  We normally leave an unused digit or two
 110  * (preset to zeroes) between buf and digits, so that there is room to store
 111  * a carry out of the top digit without special pushups.  We just need to
 112  * decrement digits (and increment weight) to make room for the carry digit.
 113  * (There is no such extra space in a numeric value stored in the database,
 114  * only in a NumericVar in memory.)
 115  *
 116  * If buf is NULL then the digit buffer isn't actually palloc'd and should
 117  * not be freed --- see the constants below for an example.
 118  *
 119  * dscale, or display scale, is the nominal precision expressed as number
 120  * of digits after the decimal point (it must always be >= 0 at present).
 121  * dscale may be more than the number of physically stored fractional digits,
 122  * implying that we have suppressed storage of significant trailing zeroes.
 123  * It should never be less than the number of stored digits, since that would
 124  * imply hiding digits that are present.  NOTE that dscale is always expressed
 125  * in *decimal* digits, and so it may correspond to a fractional number of
 126  * base-NBASE digits --- divide by DEC_DIGITS to convert to NBASE digits.
 127  *
 128  * rscale, or result scale, is the target precision for a computation.
 129  * Like dscale it is expressed as number of *decimal* digits after the decimal
 130  * point, and is always >= 0 at present.
 131  * Note that rscale is not stored in variables --- it's figured on-the-fly
 132  * from the dscales of the inputs.
 133  *
 134  * NB: All the variable-level functions are written in a style that makes it
 135  * possible to give one and the same variable as argument and destination.
 136  * This is feasible because the digit buffer is separate from the variable.
 137  * ----------
 138  */
 139 typedef struct NumericVar
 140 {
 141         int                     ndigits;                /* # of digits in digits[] - can be 0! */
 142         int                     weight;                 /* weight of first digit */
 143         int                     sign;                   /* NUMERIC_POS, NUMERIC_NEG, or
 144                                                                  * NUMERIC_NAN */
 145         int                     dscale;                 /* display scale */
 146         NumericDigit *buf;                      /* start of palloc'd space for digits[] */
 147         NumericDigit *digits;           /* base-NBASE digits */
 148 } NumericVar;
 149
 150
 151 /* ----------
 152  * Some preinitialized constants
 153  * ----------
 154  */
 155 static NumericDigit const_zero_data[1] = {0};
 156 static NumericVar const_zero =
 157 {0, 0, NUMERIC_POS, 0, NULL, const_zero_data};
 158
 159 static NumericDigit const_one_data[1] = {1};
 160 static NumericVar const_one =
 161 {1, 0, NUMERIC_POS, 0, NULL, const_one_data};
 162
 163 static NumericDigit const_two_data[1] = {2};
 164 static NumericVar const_two =
 165 {1, 0, NUMERIC_POS, 0, NULL, const_two_data};
 166
 167 #if DEC_DIGITS == 4
 168 static NumericDigit const_zero_point_five_data[1] = {5000};
 169 #elif DEC_DIGITS == 2
 170 static NumericDigit const_zero_point_five_data[1] = {50};
 171 #elif DEC_DIGITS == 1
 172 static NumericDigit const_zero_point_five_data[1] = {5};
 173 #endif
 174 static NumericVar const_zero_point_five =
 175 {1, -1, NUMERIC_POS, 1, NULL, const_zero_point_five_data};
 176
 177 #if DEC_DIGITS == 4
 178 static NumericDigit const_zero_point_nine_data[1] = {9000};
 179 #elif DEC_DIGITS == 2
 180 static NumericDigit const_zero_point_nine_data[1] = {90};
 181 #elif DEC_DIGITS == 1
 182 static NumericDigit const_zero_point_nine_data[1] = {9};
 183 #endif
 184 static NumericVar const_zero_point_nine =
 185 {1, -1, NUMERIC_POS, 1, NULL, const_zero_point_nine_data};
 186
 187 #if DEC_DIGITS == 4
 188 static NumericDigit const_zero_point_01_data[1] = {100};
 189 static NumericVar const_zero_point_01 =
 190 {1, -1, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
 191 #elif DEC_DIGITS == 2
 192 static NumericDigit const_zero_point_01_data[1] = {1};
 193 static NumericVar const_zero_point_01 =
 194 {1, -1, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
 195 #elif DEC_DIGITS == 1
 196 static NumericDigit const_zero_point_01_data[1] = {1};
 197 static NumericVar const_zero_point_01 =
 198 {1, -2, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
 199 #endif
 200
 201 #if DEC_DIGITS == 4
 202 static NumericDigit const_one_point_one_data[2] = {1, 1000};
 203 #elif DEC_DIGITS == 2
 204 static NumericDigit const_one_point_one_data[2] = {1, 10};
 205 #elif DEC_DIGITS == 1
 206 static NumericDigit const_one_point_one_data[2] = {1, 1};
 207 #endif
 208 static NumericVar const_one_point_one =
 209 {2, 0, NUMERIC_POS, 1, NULL, const_one_point_one_data};
 210
 211 static NumericVar const_nan =
 212 {0, 0, NUMERIC_NAN, 0, NULL, NULL};
 213
 214 #if DEC_DIGITS == 4
 215 static const int        round_powers[4] = { 0, 1000, 100, 10 };
 216 #endif
 217
 218
 219 /* ----------
 220  * Local functions
 221  * ----------
 222  */
 223
 224 #ifdef NUMERIC_DEBUG
 225 static void dump_numeric(const char *str, Numeric num);
 226 static void dump_var(const char *str, NumericVar *var);
 227
 228 #else
 229 #define dump_numeric(s,n)
 230 #define dump_var(s,v)
 231 #endif
 232
 233 #define digitbuf_alloc(ndigits)  \
 234         ((NumericDigit *) palloc((ndigits) * sizeof(NumericDigit)))
 235 #define digitbuf_free(buf)      \
 236         do { \
 237                  if ((buf) != NULL) \
 238                          pfree(buf); \
 239         } while (0)
 240
 241 #define init_var(v)             MemSetAligned(v, 0, sizeof(NumericVar))
 242
 243 static void alloc_var(NumericVar *var, int ndigits);
 244 static void free_var(NumericVar *var);
 245 static void zero_var(NumericVar *var);
 246
 247 static void set_var_from_str(const char *str, NumericVar *dest);
 248 static void set_var_from_num(Numeric value, NumericVar *dest);
 249 static void set_var_from_var(NumericVar *value, NumericVar *dest);
 250 static char *get_str_from_var(NumericVar *var, int dscale);
 251
 252 static Numeric make_result(NumericVar *var);
 253
 254 static void apply_typmod(NumericVar *var, int32 typmod);
 255
 256 static bool numericvar_to_int8(NumericVar *var, int64 *result);
 257 static void int8_to_numericvar(int64 val, NumericVar *var);
 258 static double numeric_to_double_no_overflow(Numeric num);
 259 static double numericvar_to_double_no_overflow(NumericVar *var);
 260
 261 static int      cmp_numerics(Numeric num1, Numeric num2);
 262 static int      cmp_var(NumericVar *var1, NumericVar *var2);
 263 static void add_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
 264 static void sub_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
 265 static void mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 266                                         int rscale);
 267 static void div_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 268                                         int rscale);
 269 static int      select_div_scale(NumericVar *var1, NumericVar *var2);
 270 static void mod_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
 271 static void ceil_var(NumericVar *var, NumericVar *result);
 272 static void floor_var(NumericVar *var, NumericVar *result);
 273
 274 static void sqrt_var(NumericVar *arg, NumericVar *result, int rscale);
 275 static void exp_var(NumericVar *arg, NumericVar *result, int rscale);
 276 static void exp_var_internal(NumericVar *arg, NumericVar *result, int rscale);
 277 static void ln_var(NumericVar *arg, NumericVar *result, int rscale);
 278 static void log_var(NumericVar *base, NumericVar *num, NumericVar *result);
 279 static void power_var(NumericVar *base, NumericVar *exp, NumericVar *result);
 280 static void power_var_int(NumericVar *base, int exp, NumericVar *result,
 281                                                   int rscale);
 282
 283 static int      cmp_abs(NumericVar *var1, NumericVar *var2);
 284 static void add_abs(NumericVar *var1, NumericVar *var2, NumericVar *result);
 285 static void sub_abs(NumericVar *var1, NumericVar *var2, NumericVar *result);
 286 static void round_var(NumericVar *var, int rscale);
 287 static void trunc_var(NumericVar *var, int rscale);
 288 static void strip_var(NumericVar *var);
 289
 290
 291 /* ----------------------------------------------------------------------
 292  *
 293  * Input-, output- and rounding-functions
 294  *
 295  * ----------------------------------------------------------------------
 296  */
 297
 298
 299 /*
 300  * numeric_in() -
 301  *
 302  *      Input function for numeric data type
 303  */
 304 Datum
 305 numeric_in(PG_FUNCTION_ARGS)
 306 {
 307         char       *str = PG_GETARG_CSTRING(0);
 308
 309 #ifdef NOT_USED
 310         Oid                     typelem = PG_GETARG_OID(1);
 311 #endif
 312         int32           typmod = PG_GETARG_INT32(2);
 313         NumericVar      value;
 314         Numeric         res;
 315
 316         /*
 317          * Check for NaN
 318          */
 319         if (strcmp(str, "NaN") == 0)
 320                 PG_RETURN_NUMERIC(make_result(&const_nan));
 321
 322         /*
 323          * Use set_var_from_str() to parse the input string and return it in
 324          * the packed DB storage format
 325          */
 326         init_var(&value);
 327         set_var_from_str(str, &value);
 328
 329         apply_typmod(&value, typmod);
 330
 331         res = make_result(&value);
 332         free_var(&value);
 333
 334         PG_RETURN_NUMERIC(res);
 335 }
 336
 337
 338 /*
 339  * numeric_out() -
 340  *
 341  *      Output function for numeric data type
 342  */
 343 Datum
 344 numeric_out(PG_FUNCTION_ARGS)
 345 {
 346         Numeric         num = PG_GETARG_NUMERIC(0);
 347         NumericVar      x;
 348         char       *str;
 349
 350         /*
 351          * Handle NaN
 352          */
 353         if (NUMERIC_IS_NAN(num))
 354                 PG_RETURN_CSTRING(pstrdup("NaN"));
 355
 356         /*
 357          * Get the number in the variable format.
 358          *
 359          * Even if we didn't need to change format, we'd still need to copy the
 360          * value to have a modifiable copy for rounding.  set_var_from_num()
 361          * also guarantees there is extra digit space in case we produce a
 362          * carry out from rounding.
 363          */
 364         init_var(&x);
 365         set_var_from_num(num, &x);
 366
 367         str = get_str_from_var(&x, x.dscale);
 368
 369         free_var(&x);
 370
 371         PG_RETURN_CSTRING(str);
 372 }
 373
 374 /*
 375  *              numeric_recv                    - converts external binary format to numeric
 376  *
 377  * External format is a sequence of int16's:
 378  * ndigits, weight, sign, dscale, NumericDigits.
 379  */
 380 Datum
 381 numeric_recv(PG_FUNCTION_ARGS)
 382 {
 383         StringInfo      buf = (StringInfo) PG_GETARG_POINTER(0);
 384         NumericVar      value;
 385         Numeric         res;
 386         int                     len,
 387                                 i;
 388
 389         init_var(&value);
 390
 391         len = (uint16) pq_getmsgint(buf, sizeof(uint16));
 392         if (len < 0 || len > NUMERIC_MAX_PRECISION + NUMERIC_MAX_RESULT_SCALE)
 393                 ereport(ERROR,
 394                                 (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
 395                                  errmsg("invalid length in external numeric")));
 396
 397         alloc_var(&value, len);
 398
 399         value.weight = (int16) pq_getmsgint(buf, sizeof(int16));
 400         value.sign = (uint16) pq_getmsgint(buf, sizeof(uint16));
 401         if (!(value.sign == NUMERIC_POS ||
 402                   value.sign == NUMERIC_NEG ||
 403                   value.sign == NUMERIC_NAN))
 404                 ereport(ERROR,
 405                                 (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
 406                                  errmsg("invalid sign in external numeric")));
 407
 408         value.dscale = (uint16) pq_getmsgint(buf, sizeof(uint16));
 409         for (i = 0; i < len; i++)
 410         {
 411                 NumericDigit    d = pq_getmsgint(buf, sizeof(NumericDigit));
 412
 413                 if (d < 0 || d >= NBASE)
 414                         ereport(ERROR,
 415                                         (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
 416                                          errmsg("invalid digit in external numeric")));
 417                 value.digits[i] = d;
 418         }
 419
 420         res = make_result(&value);
 421         free_var(&value);
 422
 423         PG_RETURN_NUMERIC(res);
 424 }
 425
 426 /*
 427  *              numeric_send                    - converts numeric to binary format
 428  */
 429 Datum
 430 numeric_send(PG_FUNCTION_ARGS)
 431 {
 432         Numeric         num = PG_GETARG_NUMERIC(0);
 433         NumericVar      x;
 434         StringInfoData buf;
 435         int                     i;
 436
 437         init_var(&x);
 438         set_var_from_num(num, &x);
 439
 440         pq_begintypsend(&buf);
 441
 442         pq_sendint(&buf, x.ndigits, sizeof(int16));
 443         pq_sendint(&buf, x.weight, sizeof(int16));
 444         pq_sendint(&buf, x.sign, sizeof(int16));
 445         pq_sendint(&buf, x.dscale, sizeof(int16));
 446         for (i = 0; i < x.ndigits; i++)
 447                 pq_sendint(&buf, x.digits[i], sizeof(NumericDigit));
 448
 449         free_var(&x);
 450
 451         PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
 452 }
 453
 454
 455 /*
 456  * numeric() -
 457  *
 458  *      This is a special function called by the Postgres database system
 459  *      before a value is stored in a tuple's attribute. The precision and
 460  *      scale of the attribute have to be applied on the value.
 461  */
 462 Datum
 463 numeric(PG_FUNCTION_ARGS)
 464 {
 465         Numeric         num = PG_GETARG_NUMERIC(0);
 466         int32           typmod = PG_GETARG_INT32(1);
 467         Numeric         new;
 468         int32           tmp_typmod;
 469         int                     precision;
 470         int                     scale;
 471         int                     ddigits;
 472         int                     maxdigits;
 473         NumericVar      var;
 474
 475         /*
 476          * Handle NaN
 477          */
 478         if (NUMERIC_IS_NAN(num))
 479                 PG_RETURN_NUMERIC(make_result(&const_nan));
 480
 481         /*
 482          * If the value isn't a valid type modifier, simply return a copy of
 483          * the input value
 484          */
 485         if (typmod < (int32) (VARHDRSZ))
 486         {
 487                 new = (Numeric) palloc(num->varlen);
 488                 memcpy(new, num, num->varlen);
 489                 PG_RETURN_NUMERIC(new);
 490         }
 491
 492         /*
 493          * Get the precision and scale out of the typmod value
 494          */
 495         tmp_typmod = typmod - VARHDRSZ;
 496         precision = (tmp_typmod >> 16) & 0xffff;
 497         scale = tmp_typmod & 0xffff;
 498         maxdigits = precision - scale;
 499
 500         /*
 501          * If the number is certainly in bounds and due to the target scale no
 502          * rounding could be necessary, just make a copy of the input and
 503          * modify its scale fields.  (Note we assume the existing dscale is
 504          * honest...)
 505          */
 506         ddigits = (num->n_weight + 1) * DEC_DIGITS;
 507         if (ddigits <= maxdigits && scale >= NUMERIC_DSCALE(num))
 508         {
 509                 new = (Numeric) palloc(num->varlen);
 510                 memcpy(new, num, num->varlen);
 511                 new->n_sign_dscale = NUMERIC_SIGN(new) |
 512                         ((uint16) scale & NUMERIC_DSCALE_MASK);
 513                 PG_RETURN_NUMERIC(new);
 514         }
 515
 516         /*
 517          * We really need to fiddle with things - unpack the number into a
 518          * variable and let apply_typmod() do it.
 519          */
 520         init_var(&var);
 521
 522         set_var_from_num(num, &var);
 523         apply_typmod(&var, typmod);
 524         new = make_result(&var);
 525
 526         free_var(&var);
 527
 528         PG_RETURN_NUMERIC(new);
 529 }
 530
 531
 532 /* ----------------------------------------------------------------------
 533  *
 534  * Sign manipulation, rounding and the like
 535  *
 536  * ----------------------------------------------------------------------
 537  */
 538
 539 Datum
 540 numeric_abs(PG_FUNCTION_ARGS)
 541 {
 542         Numeric         num = PG_GETARG_NUMERIC(0);
 543         Numeric         res;
 544
 545         /*
 546          * Handle NaN
 547          */
 548         if (NUMERIC_IS_NAN(num))
 549                 PG_RETURN_NUMERIC(make_result(&const_nan));
 550
 551         /*
 552          * Do it the easy way directly on the packed format
 553          */
 554         res = (Numeric) palloc(num->varlen);
 555         memcpy(res, num, num->varlen);
 556
 557         res->n_sign_dscale = NUMERIC_POS | NUMERIC_DSCALE(num);
 558
 559         PG_RETURN_NUMERIC(res);
 560 }
 561
 562
 563 Datum
 564 numeric_uminus(PG_FUNCTION_ARGS)
 565 {
 566         Numeric         num = PG_GETARG_NUMERIC(0);
 567         Numeric         res;
 568
 569         /*
 570          * Handle NaN
 571          */
 572         if (NUMERIC_IS_NAN(num))
 573                 PG_RETURN_NUMERIC(make_result(&const_nan));
 574
 575         /*
 576          * Do it the easy way directly on the packed format
 577          */
 578         res = (Numeric) palloc(num->varlen);
 579         memcpy(res, num, num->varlen);
 580
 581         /*
 582          * The packed format is known to be totally zero digit trimmed always.
 583          * So we can identify a ZERO by the fact that there are no digits at
 584          * all.  Do nothing to a zero.
 585          */
 586         if (num->varlen != NUMERIC_HDRSZ)
 587         {
 588                 /* Else, flip the sign */
 589                 if (NUMERIC_SIGN(num) == NUMERIC_POS)
 590                         res->n_sign_dscale = NUMERIC_NEG | NUMERIC_DSCALE(num);
 591                 else
 592                         res->n_sign_dscale = NUMERIC_POS | NUMERIC_DSCALE(num);
 593         }
 594
 595         PG_RETURN_NUMERIC(res);
 596 }
 597
 598
 599 Datum
 600 numeric_uplus(PG_FUNCTION_ARGS)
 601 {
 602         Numeric         num = PG_GETARG_NUMERIC(0);
 603         Numeric         res;
 604
 605         res = (Numeric) palloc(num->varlen);
 606         memcpy(res, num, num->varlen);
 607
 608         PG_RETURN_NUMERIC(res);
 609 }
 610
 611 /*
 612  * numeric_sign() -
 613  *
 614  * returns -1 if the argument is less than 0, 0 if the argument is equal
 615  * to 0, and 1 if the argument is greater than zero.
 616  */
 617 Datum
 618 numeric_sign(PG_FUNCTION_ARGS)
 619 {
 620         Numeric         num = PG_GETARG_NUMERIC(0);
 621         Numeric         res;
 622         NumericVar      result;
 623
 624         /*
 625          * Handle NaN
 626          */
 627         if (NUMERIC_IS_NAN(num))
 628                 PG_RETURN_NUMERIC(make_result(&const_nan));
 629
 630         init_var(&result);
 631
 632         /*
 633          * The packed format is known to be totally zero digit trimmed always.
 634          * So we can identify a ZERO by the fact that there are no digits at
 635          * all.
 636          */
 637         if (num->varlen == NUMERIC_HDRSZ)
 638                 set_var_from_var(&const_zero, &result);
 639         else
 640         {
 641                 /*
 642                  * And if there are some, we return a copy of ONE with the sign of
 643                  * our argument
 644                  */
 645                 set_var_from_var(&const_one, &result);
 646                 result.sign = NUMERIC_SIGN(num);
 647         }
 648
 649         res = make_result(&result);
 650         free_var(&result);
 651
 652         PG_RETURN_NUMERIC(res);
 653 }
 654
 655
 656 /*
 657  * numeric_round() -
 658  *
 659  *      Round a value to have 'scale' digits after the decimal point.
 660  *      We allow negative 'scale', implying rounding before the decimal
 661  *      point --- Oracle interprets rounding that way.
 662  */
 663 Datum
 664 numeric_round(PG_FUNCTION_ARGS)
 665 {
 666         Numeric         num = PG_GETARG_NUMERIC(0);
 667         int32           scale = PG_GETARG_INT32(1);
 668         Numeric         res;
 669         NumericVar      arg;
 670
 671         /*
 672          * Handle NaN
 673          */
 674         if (NUMERIC_IS_NAN(num))
 675                 PG_RETURN_NUMERIC(make_result(&const_nan));
 676
 677         /*
 678          * Limit the scale value to avoid possible overflow in calculations
 679          */
 680         scale = Max(scale, -NUMERIC_MAX_RESULT_SCALE);
 681         scale = Min(scale, NUMERIC_MAX_RESULT_SCALE);
 682
 683         /*
 684          * Unpack the argument and round it at the proper digit position
 685          */
 686         init_var(&arg);
 687         set_var_from_num(num, &arg);
 688
 689         round_var(&arg, scale);
 690
 691         /* We don't allow negative output dscale */
 692         if (scale < 0)
 693                 arg.dscale = 0;
 694
 695         /*
 696          * Return the rounded result
 697          */
 698         res = make_result(&arg);
 699
 700         free_var(&arg);
 701         PG_RETURN_NUMERIC(res);
 702 }
 703
 704
 705 /*
 706  * numeric_trunc() -
 707  *
 708  *      Truncate a value to have 'scale' digits after the decimal point.
 709  *      We allow negative 'scale', implying a truncation before the decimal
 710  *      point --- Oracle interprets truncation that way.
 711  */
 712 Datum
 713 numeric_trunc(PG_FUNCTION_ARGS)
 714 {
 715         Numeric         num = PG_GETARG_NUMERIC(0);
 716         int32           scale = PG_GETARG_INT32(1);
 717         Numeric         res;
 718         NumericVar      arg;
 719
 720         /*
 721          * Handle NaN
 722          */
 723         if (NUMERIC_IS_NAN(num))
 724                 PG_RETURN_NUMERIC(make_result(&const_nan));
 725
 726         /*
 727          * Limit the scale value to avoid possible overflow in calculations
 728          */
 729         scale = Max(scale, -NUMERIC_MAX_RESULT_SCALE);
 730         scale = Min(scale, NUMERIC_MAX_RESULT_SCALE);
 731
 732         /*
 733          * Unpack the argument and truncate it at the proper digit position
 734          */
 735         init_var(&arg);
 736         set_var_from_num(num, &arg);
 737
 738         trunc_var(&arg, scale);
 739
 740         /* We don't allow negative output dscale */
 741         if (scale < 0)
 742                 arg.dscale = 0;
 743
 744         /*
 745          * Return the truncated result
 746          */
 747         res = make_result(&arg);
 748
 749         free_var(&arg);
 750         PG_RETURN_NUMERIC(res);
 751 }
 752
 753
 754 /*
 755  * numeric_ceil() -
 756  *
 757  *      Return the smallest integer greater than or equal to the argument
 758  */
 759 Datum
 760 numeric_ceil(PG_FUNCTION_ARGS)
 761 {
 762         Numeric         num = PG_GETARG_NUMERIC(0);
 763         Numeric         res;
 764         NumericVar      result;
 765
 766         if (NUMERIC_IS_NAN(num))
 767                 PG_RETURN_NUMERIC(make_result(&const_nan));
 768
 769         init_var(&result);
 770
 771         set_var_from_num(num, &result);
 772         ceil_var(&result, &result);
 773
 774         res = make_result(&result);
 775         free_var(&result);
 776
 777         PG_RETURN_NUMERIC(res);
 778 }
 779
 780
 781 /*
 782  * numeric_floor() -
 783  *
 784  *      Return the largest integer equal to or less than the argument
 785  */
 786 Datum
 787 numeric_floor(PG_FUNCTION_ARGS)
 788 {
 789         Numeric         num = PG_GETARG_NUMERIC(0);
 790         Numeric         res;
 791         NumericVar      result;
 792
 793         if (NUMERIC_IS_NAN(num))
 794                 PG_RETURN_NUMERIC(make_result(&const_nan));
 795
 796         init_var(&result);
 797
 798         set_var_from_num(num, &result);
 799         floor_var(&result, &result);
 800
 801         res = make_result(&result);
 802         free_var(&result);
 803
 804         PG_RETURN_NUMERIC(res);
 805 }
 806
 807
 808 /* ----------------------------------------------------------------------
 809  *
 810  * Comparison functions
 811  *
 812  * Note: btree indexes need these routines not to leak memory; therefore,
 813  * be careful to free working copies of toasted datums.  Most places don't
 814  * need to be so careful.
 815  * ----------------------------------------------------------------------
 816  */
 817
 818
 819 Datum
 820 numeric_cmp(PG_FUNCTION_ARGS)
 821 {
 822         Numeric         num1 = PG_GETARG_NUMERIC(0);
 823         Numeric         num2 = PG_GETARG_NUMERIC(1);
 824         int                     result;
 825
 826         result = cmp_numerics(num1, num2);
 827
 828         PG_FREE_IF_COPY(num1, 0);
 829         PG_FREE_IF_COPY(num2, 1);
 830
 831         PG_RETURN_INT32(result);
 832 }
 833
 834
 835 Datum
 836 numeric_eq(PG_FUNCTION_ARGS)
 837 {
 838         Numeric         num1 = PG_GETARG_NUMERIC(0);
 839         Numeric         num2 = PG_GETARG_NUMERIC(1);
 840         bool            result;
 841
 842         result = cmp_numerics(num1, num2) == 0;
 843
 844         PG_FREE_IF_COPY(num1, 0);
 845         PG_FREE_IF_COPY(num2, 1);
 846
 847         PG_RETURN_BOOL(result);
 848 }
 849
 850 Datum
 851 numeric_ne(PG_FUNCTION_ARGS)
 852 {
 853         Numeric         num1 = PG_GETARG_NUMERIC(0);
 854         Numeric         num2 = PG_GETARG_NUMERIC(1);
 855         bool            result;
 856
 857         result = cmp_numerics(num1, num2) != 0;
 858
 859         PG_FREE_IF_COPY(num1, 0);
 860         PG_FREE_IF_COPY(num2, 1);
 861
 862         PG_RETURN_BOOL(result);
 863 }
 864
 865 Datum
 866 numeric_gt(PG_FUNCTION_ARGS)
 867 {
 868         Numeric         num1 = PG_GETARG_NUMERIC(0);
 869         Numeric         num2 = PG_GETARG_NUMERIC(1);
 870         bool            result;
 871
 872         result = cmp_numerics(num1, num2) > 0;
 873
 874         PG_FREE_IF_COPY(num1, 0);
 875         PG_FREE_IF_COPY(num2, 1);
 876
 877         PG_RETURN_BOOL(result);
 878 }
 879
 880 Datum
 881 numeric_ge(PG_FUNCTION_ARGS)
 882 {
 883         Numeric         num1 = PG_GETARG_NUMERIC(0);
 884         Numeric         num2 = PG_GETARG_NUMERIC(1);
 885         bool            result;
 886
 887         result = cmp_numerics(num1, num2) >= 0;
 888
 889         PG_FREE_IF_COPY(num1, 0);
 890         PG_FREE_IF_COPY(num2, 1);
 891
 892         PG_RETURN_BOOL(result);
 893 }
 894
 895 Datum
 896 numeric_lt(PG_FUNCTION_ARGS)
 897 {
 898         Numeric         num1 = PG_GETARG_NUMERIC(0);
 899         Numeric         num2 = PG_GETARG_NUMERIC(1);
 900         bool            result;
 901
 902         result = cmp_numerics(num1, num2) < 0;
 903
 904         PG_FREE_IF_COPY(num1, 0);
 905         PG_FREE_IF_COPY(num2, 1);
 906
 907         PG_RETURN_BOOL(result);
 908 }
 909
 910 Datum
 911 numeric_le(PG_FUNCTION_ARGS)
 912 {
 913         Numeric         num1 = PG_GETARG_NUMERIC(0);
 914         Numeric         num2 = PG_GETARG_NUMERIC(1);
 915         bool            result;
 916
 917         result = cmp_numerics(num1, num2) <= 0;
 918
 919         PG_FREE_IF_COPY(num1, 0);
 920         PG_FREE_IF_COPY(num2, 1);
 921
 922         PG_RETURN_BOOL(result);
 923 }
 924
 925 static int
 926 cmp_numerics(Numeric num1, Numeric num2)
 927 {
 928         int                     result;
 929
 930         /*
 931          * We consider all NANs to be equal and larger than any non-NAN. This
 932          * is somewhat arbitrary; the important thing is to have a consistent
 933          * sort order.
 934          */
 935         if (NUMERIC_IS_NAN(num1))
 936         {
 937                 if (NUMERIC_IS_NAN(num2))
 938                         result = 0;                     /* NAN = NAN */
 939                 else
 940                         result = 1;                     /* NAN > non-NAN */
 941         }
 942         else if (NUMERIC_IS_NAN(num2))
 943         {
 944                 result = -1;                    /* non-NAN < NAN */
 945         }
 946         else
 947         {
 948                 NumericVar      arg1;
 949                 NumericVar      arg2;
 950
 951                 init_var(&arg1);
 952                 init_var(&arg2);
 953
 954                 set_var_from_num(num1, &arg1);
 955                 set_var_from_num(num2, &arg2);
 956
 957                 result = cmp_var(&arg1, &arg2);
 958
 959                 free_var(&arg1);
 960                 free_var(&arg2);
 961         }
 962
 963         return result;
 964 }
 965
 966
 967 /* ----------------------------------------------------------------------
 968  *
 969  * Basic arithmetic functions
 970  *
 971  * ----------------------------------------------------------------------
 972  */
 973
 974
 975 /*
 976  * numeric_add() -
 977  *
 978  *      Add two numerics
 979  */
 980 Datum
 981 numeric_add(PG_FUNCTION_ARGS)
 982 {
 983         Numeric         num1 = PG_GETARG_NUMERIC(0);
 984         Numeric         num2 = PG_GETARG_NUMERIC(1);
 985         NumericVar      arg1;
 986         NumericVar      arg2;
 987         NumericVar      result;
 988         Numeric         res;
 989
 990         /*
 991          * Handle NaN
 992          */
 993         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
 994                 PG_RETURN_NUMERIC(make_result(&const_nan));
 995
 996         /*
 997          * Unpack the values, let add_var() compute the result and return it.
 998          */
 999         init_var(&arg1);
1000         init_var(&arg2);
1001         init_var(&result);
1002
1003         set_var_from_num(num1, &arg1);
1004         set_var_from_num(num2, &arg2);
1005
1006         add_var(&arg1, &arg2, &result);
1007
1008         res = make_result(&result);
1009
1010         free_var(&arg1);
1011         free_var(&arg2);
1012         free_var(&result);
1013
1014         PG_RETURN_NUMERIC(res);
1015 }
1016
1017
1018 /*
1019  * numeric_sub() -
1020  *
1021  *      Subtract one numeric from another
1022  */
1023 Datum
1024 numeric_sub(PG_FUNCTION_ARGS)
1025 {
1026         Numeric         num1 = PG_GETARG_NUMERIC(0);
1027         Numeric         num2 = PG_GETARG_NUMERIC(1);
1028         NumericVar      arg1;
1029         NumericVar      arg2;
1030         NumericVar      result;
1031         Numeric         res;
1032
1033         /*
1034          * Handle NaN
1035          */
1036         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1037                 PG_RETURN_NUMERIC(make_result(&const_nan));
1038
1039         /*
1040          * Unpack the values, let sub_var() compute the result and return it.
1041          */
1042         init_var(&arg1);
1043         init_var(&arg2);
1044         init_var(&result);
1045
1046         set_var_from_num(num1, &arg1);
1047         set_var_from_num(num2, &arg2);
1048
1049         sub_var(&arg1, &arg2, &result);
1050
1051         res = make_result(&result);
1052
1053         free_var(&arg1);
1054         free_var(&arg2);
1055         free_var(&result);
1056
1057         PG_RETURN_NUMERIC(res);
1058 }
1059
1060
1061 /*
1062  * numeric_mul() -
1063  *
1064  *      Calculate the product of two numerics
1065  */
1066 Datum
1067 numeric_mul(PG_FUNCTION_ARGS)
1068 {
1069         Numeric         num1 = PG_GETARG_NUMERIC(0);
1070         Numeric         num2 = PG_GETARG_NUMERIC(1);
1071         NumericVar      arg1;
1072         NumericVar      arg2;
1073         NumericVar      result;
1074         Numeric         res;
1075
1076         /*
1077          * Handle NaN
1078          */
1079         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1080                 PG_RETURN_NUMERIC(make_result(&const_nan));
1081
1082         /*
1083          * Unpack the values, let mul_var() compute the result and return it.
1084          * Unlike add_var() and sub_var(), mul_var() will round its result.
1085          * In the case of numeric_mul(), which is invoked for the * operator on
1086          * numerics, we request exact representation for the product (rscale =
1087          * sum(dscale of arg1, dscale of arg2)).
1088          */
1089         init_var(&arg1);
1090         init_var(&arg2);
1091         init_var(&result);
1092
1093         set_var_from_num(num1, &arg1);
1094         set_var_from_num(num2, &arg2);
1095
1096         mul_var(&arg1, &arg2, &result, arg1.dscale + arg2.dscale);
1097
1098         res = make_result(&result);
1099
1100         free_var(&arg1);
1101         free_var(&arg2);
1102         free_var(&result);
1103
1104         PG_RETURN_NUMERIC(res);
1105 }
1106
1107
1108 /*
1109  * numeric_div() -
1110  *
1111  *      Divide one numeric into another
1112  */
1113 Datum
1114 numeric_div(PG_FUNCTION_ARGS)
1115 {
1116         Numeric         num1 = PG_GETARG_NUMERIC(0);
1117         Numeric         num2 = PG_GETARG_NUMERIC(1);
1118         NumericVar      arg1;
1119         NumericVar      arg2;
1120         NumericVar      result;
1121         Numeric         res;
1122         int                     rscale;
1123
1124         /*
1125          * Handle NaN
1126          */
1127         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1128                 PG_RETURN_NUMERIC(make_result(&const_nan));
1129
1130         /*
1131          * Unpack the arguments
1132          */
1133         init_var(&arg1);
1134         init_var(&arg2);
1135         init_var(&result);
1136
1137         set_var_from_num(num1, &arg1);
1138         set_var_from_num(num2, &arg2);
1139
1140         /*
1141          * Select scale for division result
1142          */
1143         rscale = select_div_scale(&arg1, &arg2);
1144
1145         /*
1146          * Do the divide and return the result
1147          */
1148         div_var(&arg1, &arg2, &result, rscale);
1149
1150         res = make_result(&result);
1151
1152         free_var(&arg1);
1153         free_var(&arg2);
1154         free_var(&result);
1155
1156         PG_RETURN_NUMERIC(res);
1157 }
1158
1159
1160 /*
1161  * numeric_mod() -
1162  *
1163  *      Calculate the modulo of two numerics
1164  */
1165 Datum
1166 numeric_mod(PG_FUNCTION_ARGS)
1167 {
1168         Numeric         num1 = PG_GETARG_NUMERIC(0);
1169         Numeric         num2 = PG_GETARG_NUMERIC(1);
1170         Numeric         res;
1171         NumericVar      arg1;
1172         NumericVar      arg2;
1173         NumericVar      result;
1174
1175         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1176                 PG_RETURN_NUMERIC(make_result(&const_nan));
1177
1178         init_var(&arg1);
1179         init_var(&arg2);
1180         init_var(&result);
1181
1182         set_var_from_num(num1, &arg1);
1183         set_var_from_num(num2, &arg2);
1184
1185         mod_var(&arg1, &arg2, &result);
1186
1187         res = make_result(&result);
1188
1189         free_var(&result);
1190         free_var(&arg2);
1191         free_var(&arg1);
1192
1193         PG_RETURN_NUMERIC(res);
1194 }
1195
1196
1197 /*
1198  * numeric_inc() -
1199  *
1200  *      Increment a number by one
1201  */
1202 Datum
1203 numeric_inc(PG_FUNCTION_ARGS)
1204 {
1205         Numeric         num = PG_GETARG_NUMERIC(0);
1206         NumericVar      arg;
1207         Numeric         res;
1208
1209         /*
1210          * Handle NaN
1211          */
1212         if (NUMERIC_IS_NAN(num))
1213                 PG_RETURN_NUMERIC(make_result(&const_nan));
1214
1215         /*
1216          * Compute the result and return it
1217          */
1218         init_var(&arg);
1219
1220         set_var_from_num(num, &arg);
1221
1222         add_var(&arg, &const_one, &arg);
1223
1224         res = make_result(&arg);
1225
1226         free_var(&arg);
1227
1228         PG_RETURN_NUMERIC(res);
1229 }
1230
1231
1232 /*
1233  * numeric_smaller() -
1234  *
1235  *      Return the smaller of two numbers
1236  */
1237 Datum
1238 numeric_smaller(PG_FUNCTION_ARGS)
1239 {
1240         Numeric         num1 = PG_GETARG_NUMERIC(0);
1241         Numeric         num2 = PG_GETARG_NUMERIC(1);
1242         NumericVar      arg1;
1243         NumericVar      arg2;
1244         Numeric         res;
1245
1246         /*
1247          * Handle NaN
1248          */
1249         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1250                 PG_RETURN_NUMERIC(make_result(&const_nan));
1251
1252         /*
1253          * Unpack the values, and decide which is the smaller one
1254          */
1255         init_var(&arg1);
1256         init_var(&arg2);
1257
1258         set_var_from_num(num1, &arg1);
1259         set_var_from_num(num2, &arg2);
1260
1261         if (cmp_var(&arg1, &arg2) <= 0)
1262                 res = make_result(&arg1);
1263         else
1264                 res = make_result(&arg2);
1265
1266         free_var(&arg1);
1267         free_var(&arg2);
1268
1269         PG_RETURN_NUMERIC(res);
1270 }
1271
1272
1273 /*
1274  * numeric_larger() -
1275  *
1276  *      Return the larger of two numbers
1277  */
1278 Datum
1279 numeric_larger(PG_FUNCTION_ARGS)
1280 {
1281         Numeric         num1 = PG_GETARG_NUMERIC(0);
1282         Numeric         num2 = PG_GETARG_NUMERIC(1);
1283         NumericVar      arg1;
1284         NumericVar      arg2;
1285         Numeric         res;
1286
1287         /*
1288          * Handle NaN
1289          */
1290         if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
1291                 PG_RETURN_NUMERIC(make_result(&const_nan));
1292
1293         /*
1294          * Unpack the values, and decide which is the larger one
1295          */
1296         init_var(&arg1);
1297         init_var(&arg2);
1298
1299         set_var_from_num(num1, &arg1);
1300         set_var_from_num(num2, &arg2);
1301
1302         if (cmp_var(&arg1, &arg2) >= 0)
1303                 res = make_result(&arg1);
1304         else
1305                 res = make_result(&arg2);
1306
1307         free_var(&arg1);
1308         free_var(&arg2);
1309
1310         PG_RETURN_NUMERIC(res);
1311 }
1312
1313
1314 /* ----------------------------------------------------------------------
1315  *
1316  * Advanced math functions
1317  *
1318  * ----------------------------------------------------------------------
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 decimal
1411          * 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(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(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((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
2451 #endif   /* NUMERIC_DEBUG */
2452
2453
2454 /* ----------------------------------------------------------------------
2455  *
2456  * Local functions follow
2457  *
2458  * In general, these do not support NaNs --- callers must eliminate
2459  * the possibility of NaN first.  (make_result() is an exception.)
2460  *
2461  * ----------------------------------------------------------------------
2462  */
2463
2464
2465 /*
2466  * alloc_var() -
2467  *
2468  *      Allocate a digit buffer of ndigits digits (plus a spare digit for rounding)
2469  */
2470 static void
2471 alloc_var(NumericVar *var, int ndigits)
2472 {
2473         digitbuf_free(var->buf);
2474         var->buf = digitbuf_alloc(ndigits + 1);
2475         var->buf[0] = 0;                                /* spare digit for rounding */
2476         var->digits = var->buf + 1;
2477         var->ndigits = ndigits;
2478 }
2479
2480
2481 /*
2482  * free_var() -
2483  *
2484  *      Return the digit buffer of a variable to the free pool
2485  */
2486 static void
2487 free_var(NumericVar *var)
2488 {
2489         digitbuf_free(var->buf);
2490         var->buf = NULL;
2491         var->digits = NULL;
2492         var->sign = NUMERIC_NAN;
2493 }
2494
2495
2496 /*
2497  * zero_var() -
2498  *
2499  *      Set a variable to ZERO.
2500  *      Note: its dscale is not touched.
2501  */
2502 static void
2503 zero_var(NumericVar *var)
2504 {
2505         digitbuf_free(var->buf);
2506         var->buf = NULL;
2507         var->digits = NULL;
2508         var->ndigits = 0;
2509         var->weight = 0;                        /* by convention; doesn't really matter */
2510         var->sign = NUMERIC_POS;        /* anything but NAN... */
2511 }
2512
2513
2514 /*
2515  * set_var_from_str()
2516  *
2517  *      Parse a string and put the number into a variable
2518  */
2519 static void
2520 set_var_from_str(const char *str, NumericVar *dest)
2521 {
2522         const char *cp = str;
2523         bool            have_dp = FALSE;
2524         int                     i;
2525         unsigned char *decdigits;
2526         int                     sign = NUMERIC_POS;
2527         int                     dweight = -1;
2528         int                     ddigits;
2529         int                     dscale = 0;
2530         int                     weight;
2531         int                     ndigits;
2532         int                     offset;
2533         NumericDigit *digits;
2534
2535         /*
2536          * We first parse the string to extract decimal digits and determine the
2537          * correct decimal weight.  Then convert to NBASE representation.
2538          */
2539
2540         /* skip leading spaces */
2541         while (*cp)
2542         {
2543                 if (!isspace((unsigned char) *cp))
2544                         break;
2545                 cp++;
2546         }
2547
2548         switch (*cp)
2549         {
2550                 case '+':
2551                         sign = NUMERIC_POS;
2552                         cp++;
2553                         break;
2554
2555                 case '-':
2556                         sign = NUMERIC_NEG;
2557                         cp++;
2558                         break;
2559         }
2560
2561         if (*cp == '.')
2562         {
2563                 have_dp = TRUE;
2564                 cp++;
2565         }
2566
2567         if (!isdigit((unsigned char) *cp))
2568                 ereport(ERROR,
2569                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2570                                  errmsg("invalid input syntax for numeric: \"%s\"", str)));
2571
2572         decdigits = (unsigned char *) palloc(strlen(cp) + DEC_DIGITS*2);
2573
2574         /* leading padding for digit alignment later */
2575         memset(decdigits, 0, DEC_DIGITS);
2576         i = DEC_DIGITS;
2577
2578         while (*cp)
2579         {
2580                 if (isdigit((unsigned char) *cp))
2581                 {
2582                         decdigits[i++] = *cp++ - '0';
2583                         if (!have_dp)
2584                                 dweight++;
2585                         else
2586                                 dscale++;
2587                 }
2588                 else if (*cp == '.')
2589                 {
2590                         if (have_dp)
2591                                 ereport(ERROR,
2592                                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2593                                                  errmsg("invalid input syntax for numeric: \"%s\"",
2594                                                                 str)));
2595                         have_dp = TRUE;
2596                         cp++;
2597                 }
2598                 else
2599                         break;
2600         }
2601
2602         ddigits = i - DEC_DIGITS;
2603         /* trailing padding for digit alignment later */
2604         memset(decdigits + i, 0, DEC_DIGITS-1);
2605
2606         /* Handle exponent, if any */
2607         if (*cp == 'e' || *cp == 'E')
2608         {
2609                 long            exponent;
2610                 char       *endptr;
2611
2612                 cp++;
2613                 exponent = strtol(cp, &endptr, 10);
2614                 if (endptr == cp)
2615                         ereport(ERROR,
2616                                         (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2617                                          errmsg("invalid input syntax for numeric: \"%s\"",
2618                                                         str)));
2619                 cp = endptr;
2620                 if (exponent > NUMERIC_MAX_PRECISION ||
2621                         exponent < -NUMERIC_MAX_PRECISION)
2622                         ereport(ERROR,
2623                                         (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2624                                          errmsg("invalid input syntax for numeric: \"%s\"",
2625                                                         str)));
2626                 dweight += (int) exponent;
2627                 dscale -= (int) exponent;
2628                 if (dscale < 0)
2629                         dscale = 0;
2630         }
2631
2632         /* Should be nothing left but spaces */
2633         while (*cp)
2634         {
2635                 if (!isspace((unsigned char) *cp))
2636                         ereport(ERROR,
2637                                         (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
2638                                          errmsg("invalid input syntax for numeric: \"%s\"",
2639                                                         str)));
2640                 cp++;
2641         }
2642
2643         /*
2644          * Okay, convert pure-decimal representation to base NBASE.  First we
2645          * need to determine the converted weight and ndigits.  offset is the
2646          * number of decimal zeroes to insert before the first given digit to
2647          * have a correctly aligned first NBASE digit.
2648          */
2649         if (dweight >= 0)
2650                 weight = (dweight + 1 + DEC_DIGITS-1) / DEC_DIGITS - 1;
2651         else
2652                 weight = - ((-dweight - 1) / DEC_DIGITS + 1);
2653         offset = (weight + 1) * DEC_DIGITS - (dweight + 1);
2654         ndigits = (ddigits + offset + DEC_DIGITS-1) / DEC_DIGITS;
2655
2656         alloc_var(dest, ndigits);
2657         dest->sign = sign;
2658         dest->weight = weight;
2659         dest->dscale = dscale;
2660
2661         i = DEC_DIGITS - offset;
2662         digits = dest->digits;
2663
2664         while (ndigits-- > 0)
2665         {
2666 #if DEC_DIGITS == 4
2667                 *digits++ = ((decdigits[i] * 10 + decdigits[i+1]) * 10 +
2668                                          decdigits[i+2]) * 10 + decdigits[i+3];
2669 #elif DEC_DIGITS == 2
2670                 *digits++ = decdigits[i] * 10 + decdigits[i+1];
2671 #elif DEC_DIGITS == 1
2672                 *digits++ = decdigits[i];
2673 #else
2674 #error unsupported NBASE
2675 #endif
2676                 i += DEC_DIGITS;
2677         }
2678
2679         pfree(decdigits);
2680
2681         /* Strip any leading/trailing zeroes, and normalize weight if zero */
2682         strip_var(dest);
2683 }
2684
2685
2686 /*
2687  * set_var_from_num() -
2688  *
2689  *      Convert the packed db format into a variable
2690  */
2691 static void
2692 set_var_from_num(Numeric num, NumericVar *dest)
2693 {
2694         int                     ndigits;
2695
2696         ndigits = (num->varlen - NUMERIC_HDRSZ) / sizeof(NumericDigit);
2697
2698         alloc_var(dest, ndigits);
2699
2700         dest->weight = num->n_weight;
2701         dest->sign = NUMERIC_SIGN(num);
2702         dest->dscale = NUMERIC_DSCALE(num);
2703
2704         memcpy(dest->digits, num->n_data, ndigits * sizeof(NumericDigit));
2705 }
2706
2707
2708 /*
2709  * set_var_from_var() -
2710  *
2711  *      Copy one variable into another
2712  */
2713 static void
2714 set_var_from_var(NumericVar *value, NumericVar *dest)
2715 {
2716         NumericDigit *newbuf;
2717
2718         newbuf = digitbuf_alloc(value->ndigits + 1);
2719         newbuf[0] = 0;                          /* spare digit for rounding */
2720         memcpy(newbuf + 1, value->digits, value->ndigits * sizeof(NumericDigit));
2721
2722         digitbuf_free(dest->buf);
2723
2724         memcpy(dest, value, sizeof(NumericVar));
2725         dest->buf = newbuf;
2726         dest->digits = newbuf + 1;
2727 }
2728
2729
2730 /*
2731  * get_str_from_var() -
2732  *
2733  *      Convert a var to text representation (guts of numeric_out).
2734  *      CAUTION: var's contents may be modified by rounding!
2735  *      Returns a palloc'd string.
2736  */
2737 static char *
2738 get_str_from_var(NumericVar *var, int dscale)
2739 {
2740         char       *str;
2741         char       *cp;
2742         char       *endcp;
2743         int                     i;
2744         int                     d;
2745         NumericDigit    dig;
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.
2762          * dscale is the # of decimal digits we will print after decimal point.
2763          * We may generate as many as DEC_DIGITS-1 excess digits at the end,
2764          * and in 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("ABS(value) >= 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 point.
3048          */
3049         weight = var->weight;
3050         Assert(weight >= 0 && ndigits <= weight+1);
3051
3052         /* Construct the result */
3053         digits = var->digits;
3054         neg = (var->sign == NUMERIC_NEG);
3055         val = digits[0];
3056         for (i = 1; i <= weight; i++)
3057         {
3058                 oldval = val;
3059                 val *= NBASE;
3060                 if (i < ndigits)
3061                         val += digits[i];
3062                 /*
3063                  * The overflow check is a bit tricky because we want to accept
3064                  * INT64_MIN, which will overflow the positive accumulator.  We
3065                  * can detect this case easily though because INT64_MIN is the
3066                  * only nonzero value for which -val == val (on a two's complement
3067                  * machine, anyway).
3068                  */
3069                 if ((val / NBASE) != oldval)    /* possible overflow? */
3070                 {
3071                         if (!neg || (-val) != val || val == 0 || oldval < 0)
3072                                 return false;
3073                 }
3074         }
3075
3076         *result = neg ? -val : val;
3077         return true;
3078 }
3079
3080 /*
3081  * Convert int8 value to numeric.
3082  */
3083 static void
3084 int8_to_numericvar(int64 val, NumericVar *var)
3085 {
3086         uint64          uval,
3087                                 newuval;
3088         NumericDigit *ptr;
3089         int                     ndigits;
3090
3091         /* int8 can require at most 19 decimal digits; add one for safety */
3092         alloc_var(var, 20/DEC_DIGITS);
3093         if (val < 0)
3094         {
3095                 var->sign = NUMERIC_NEG;
3096                 uval = -val;
3097         }
3098         else
3099         {
3100                 var->sign = NUMERIC_POS;
3101                 uval = val;
3102         }
3103         var->dscale = 0;
3104         if (val == 0)
3105         {
3106                 var->ndigits = 0;
3107                 var->weight = 0;
3108                 return;
3109         }
3110         ptr = var->digits + var->ndigits;
3111         ndigits = 0;
3112         do {
3113                 ptr--;
3114                 ndigits++;
3115                 newuval = uval / NBASE;
3116                 *ptr = uval - newuval * NBASE;
3117                 uval = newuval;
3118         } while (uval);
3119         var->digits = ptr;
3120         var->ndigits = ndigits;
3121         var->weight = ndigits - 1;
3122 }
3123
3124 /*
3125  * Convert numeric to float8; if out of range, return +/- HUGE_VAL
3126  */
3127 static double
3128 numeric_to_double_no_overflow(Numeric num)
3129 {
3130         char       *tmp;
3131         double          val;
3132         char       *endptr;
3133
3134         tmp = DatumGetCString(DirectFunctionCall1(numeric_out,
3135                                                                                           NumericGetDatum(num)));
3136
3137         /* unlike float8in, we ignore ERANGE from strtod */
3138         val = strtod(tmp, &endptr);
3139         if (*endptr != '\0')
3140         {
3141                 /* shouldn't happen ... */
3142                 ereport(ERROR,
3143                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
3144                                  errmsg("invalid input syntax for float8: \"%s\"",
3145                                                 tmp)));
3146         }
3147
3148         pfree(tmp);
3149
3150         return val;
3151 }
3152
3153 /* As above, but work from a NumericVar */
3154 static double
3155 numericvar_to_double_no_overflow(NumericVar *var)
3156 {
3157         char       *tmp;
3158         double          val;
3159         char       *endptr;
3160
3161         tmp = get_str_from_var(var, var->dscale);
3162
3163         /* unlike float8in, we ignore ERANGE from strtod */
3164         val = strtod(tmp, &endptr);
3165         if (*endptr != '\0')
3166         {
3167                 /* shouldn't happen ... */
3168                 ereport(ERROR,
3169                                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
3170                                  errmsg("invalid input syntax for float8: \"%s\"",
3171                                                 tmp)));
3172         }
3173
3174         pfree(tmp);
3175
3176         return val;
3177 }
3178
3179
3180 /*
3181  * cmp_var() -
3182  *
3183  *      Compare two values on variable level.  We assume zeroes have been
3184  *      truncated to no digits.
3185  */
3186 static int
3187 cmp_var(NumericVar *var1, NumericVar *var2)
3188 {
3189         if (var1->ndigits == 0)
3190         {
3191                 if (var2->ndigits == 0)
3192                         return 0;
3193                 if (var2->sign == NUMERIC_NEG)
3194                         return 1;
3195                 return -1;
3196         }
3197         if (var2->ndigits == 0)
3198         {
3199                 if (var1->sign == NUMERIC_POS)
3200                         return 1;
3201                 return -1;
3202         }
3203
3204         if (var1->sign == NUMERIC_POS)
3205         {
3206                 if (var2->sign == NUMERIC_NEG)
3207                         return 1;
3208                 return cmp_abs(var1, var2);
3209         }
3210
3211         if (var2->sign == NUMERIC_POS)
3212                 return -1;
3213
3214         return cmp_abs(var2, var1);
3215 }
3216
3217
3218 /*
3219  * add_var() -
3220  *
3221  *      Full version of add functionality on variable level (handling signs).
3222  *      result might point to one of the operands too without danger.
3223  */
3224 static void
3225 add_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
3226 {
3227         /*
3228          * Decide on the signs of the two variables what to do
3229          */
3230         if (var1->sign == NUMERIC_POS)
3231         {
3232                 if (var2->sign == NUMERIC_POS)
3233                 {
3234                         /*
3235                          * Both are positive result = +(ABS(var1) + ABS(var2))
3236                          */
3237                         add_abs(var1, var2, result);
3238                         result->sign = NUMERIC_POS;
3239                 }
3240                 else
3241                 {
3242                         /*
3243                          * var1 is positive, var2 is negative Must compare absolute
3244                          * values
3245                          */
3246                         switch (cmp_abs(var1, var2))
3247                         {
3248                                 case 0:
3249                                         /* ----------
3250                                          * ABS(var1) == ABS(var2)
3251                                          * result = ZERO
3252                                          * ----------
3253                                          */
3254                                         zero_var(result);
3255                                         result->dscale = Max(var1->dscale, var2->dscale);
3256                                         break;
3257
3258                                 case 1:
3259                                         /* ----------
3260                                          * ABS(var1) > ABS(var2)
3261                                          * result = +(ABS(var1) - ABS(var2))
3262                                          * ----------
3263                                          */
3264                                         sub_abs(var1, var2, result);
3265                                         result->sign = NUMERIC_POS;
3266                                         break;
3267
3268                                 case -1:
3269                                         /* ----------
3270                                          * ABS(var1) < ABS(var2)
3271                                          * result = -(ABS(var2) - ABS(var1))
3272                                          * ----------
3273                                          */
3274                                         sub_abs(var2, var1, result);
3275                                         result->sign = NUMERIC_NEG;
3276                                         break;
3277                         }
3278                 }
3279         }
3280         else
3281         {
3282                 if (var2->sign == NUMERIC_POS)
3283                 {
3284                         /* ----------
3285                          * var1 is negative, var2 is positive
3286                          * Must compare absolute values
3287                          * ----------
3288                          */
3289                         switch (cmp_abs(var1, var2))
3290                         {
3291                                 case 0:
3292                                         /* ----------
3293                                          * ABS(var1) == ABS(var2)
3294                                          * result = ZERO
3295                                          * ----------
3296                                          */
3297                                         zero_var(result);
3298                                         result->dscale = Max(var1->dscale, var2->dscale);
3299                                         break;
3300
3301                                 case 1:
3302                                         /* ----------
3303                                          * ABS(var1) > ABS(var2)
3304                                          * result = -(ABS(var1) - ABS(var2))
3305                                          * ----------
3306                                          */
3307                                         sub_abs(var1, var2, result);
3308                                         result->sign = NUMERIC_NEG;
3309                                         break;
3310
3311                                 case -1:
3312                                         /* ----------
3313                                          * ABS(var1) < ABS(var2)
3314                                          * result = +(ABS(var2) - ABS(var1))
3315                                          * ----------
3316                                          */
3317                                         sub_abs(var2, var1, result);
3318                                         result->sign = NUMERIC_POS;
3319                                         break;
3320                         }
3321                 }
3322                 else
3323                 {
3324                         /* ----------
3325                          * Both are negative
3326                          * result = -(ABS(var1) + ABS(var2))
3327                          * ----------
3328                          */
3329                         add_abs(var1, var2, result);
3330                         result->sign = NUMERIC_NEG;
3331                 }
3332         }
3333 }
3334
3335
3336 /*
3337  * sub_var() -
3338  *
3339  *      Full version of sub functionality on variable level (handling signs).
3340  *      result might point to one of the operands too without danger.
3341  */
3342 static void
3343 sub_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
3344 {
3345         /*
3346          * Decide on the signs of the two variables what to do
3347          */
3348         if (var1->sign == NUMERIC_POS)
3349         {
3350                 if (var2->sign == NUMERIC_NEG)
3351                 {
3352                         /* ----------
3353                          * var1 is positive, var2 is negative
3354                          * result = +(ABS(var1) + ABS(var2))
3355                          * ----------
3356                          */
3357                         add_abs(var1, var2, result);
3358                         result->sign = NUMERIC_POS;
3359                 }
3360                 else
3361                 {
3362                         /* ----------
3363                          * Both are positive
3364                          * Must compare absolute values
3365                          * ----------
3366                          */
3367                         switch (cmp_abs(var1, var2))
3368                         {
3369                                 case 0:
3370                                         /* ----------
3371                                          * ABS(var1) == ABS(var2)
3372                                          * result = ZERO
3373                                          * ----------
3374                                          */
3375                                         zero_var(result);
3376                                         result->dscale = Max(var1->dscale, var2->dscale);
3377                                         break;
3378
3379                                 case 1:
3380                                         /* ----------
3381                                          * ABS(var1) > ABS(var2)
3382                                          * result = +(ABS(var1) - ABS(var2))
3383                                          * ----------
3384                                          */
3385                                         sub_abs(var1, var2, result);
3386                                         result->sign = NUMERIC_POS;
3387                                         break;
3388
3389                                 case -1:
3390                                         /* ----------
3391                                          * ABS(var1) < ABS(var2)
3392                                          * result = -(ABS(var2) - ABS(var1))
3393                                          * ----------
3394                                          */
3395                                         sub_abs(var2, var1, result);
3396                                         result->sign = NUMERIC_NEG;
3397                                         break;
3398                         }
3399                 }
3400         }
3401         else
3402         {
3403                 if (var2->sign == NUMERIC_NEG)
3404                 {
3405                         /* ----------
3406                          * Both are negative
3407                          * Must compare absolute values
3408                          * ----------
3409                          */
3410                         switch (cmp_abs(var1, var2))
3411                         {
3412                                 case 0:
3413                                         /* ----------
3414                                          * ABS(var1) == ABS(var2)
3415                                          * result = ZERO
3416                                          * ----------
3417                                          */
3418                                         zero_var(result);
3419                                         result->dscale = Max(var1->dscale, var2->dscale);
3420                                         break;
3421
3422                                 case 1:
3423                                         /* ----------
3424                                          * ABS(var1) > ABS(var2)
3425                                          * result = -(ABS(var1) - ABS(var2))
3426                                          * ----------
3427                                          */
3428                                         sub_abs(var1, var2, result);
3429                                         result->sign = NUMERIC_NEG;
3430                                         break;
3431
3432                                 case -1:
3433                                         /* ----------
3434                                          * ABS(var1) < ABS(var2)
3435                                          * result = +(ABS(var2) - ABS(var1))
3436                                          * ----------
3437                                          */
3438                                         sub_abs(var2, var1, result);
3439                                         result->sign = NUMERIC_POS;
3440                                         break;
3441                         }
3442                 }
3443                 else
3444                 {
3445                         /* ----------
3446                          * var1 is negative, var2 is positive
3447                          * result = -(ABS(var1) + ABS(var2))
3448                          * ----------
3449                          */
3450                         add_abs(var1, var2, result);
3451                         result->sign = NUMERIC_NEG;
3452                 }
3453         }
3454 }
3455
3456
3457 /*
3458  * mul_var() -
3459  *
3460  *      Multiplication on variable level. Product of var1 * var2 is stored
3461  *      in result.  Result is rounded to no more than rscale fractional digits.
3462  */
3463 static void
3464 mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
3465                 int rscale)
3466 {
3467         int                     res_ndigits;
3468         int                     res_sign;
3469         int                     res_weight;
3470         int                     maxdigits;
3471         int                *dig;
3472         int                     carry;
3473         int                     maxdig;
3474         int                     newdig;
3475         NumericDigit *res_digits;
3476         int                     i,
3477                                 ri,
3478                                 i1,
3479                                 i2;
3480         /* copy these values into local vars for speed in inner loop */
3481         int                     var1ndigits = var1->ndigits;
3482         int                     var2ndigits = var2->ndigits;
3483         NumericDigit *var1digits = var1->digits;
3484         NumericDigit *var2digits = var2->digits;
3485
3486         if (var1ndigits == 0 || var2ndigits == 0)
3487         {
3488                 /* one or both inputs is zero; so is result */
3489                 zero_var(result);
3490                 result->dscale = rscale;
3491                 return;
3492         }
3493
3494         /* Determine result sign and (maximum possible) weight */
3495         if (var1->sign == var2->sign)
3496                 res_sign = NUMERIC_POS;
3497         else
3498                 res_sign = NUMERIC_NEG;
3499         res_weight = var1->weight + var2->weight + 2;
3500
3501         /*
3502          * Determine number of result digits to compute.  If the exact result
3503          * would have more than rscale fractional digits, truncate the computation
3504          * with MUL_GUARD_DIGITS guard digits.  We do that by pretending that
3505          * one or both inputs have fewer digits than they really do.
3506          */
3507         res_ndigits = var1ndigits + var2ndigits + 1;
3508         maxdigits = res_weight + 1 + (rscale * DEC_DIGITS) + MUL_GUARD_DIGITS;
3509         if (res_ndigits > maxdigits)
3510         {
3511                 if (maxdigits < 3)
3512                 {
3513                         /* no useful precision at all in the result... */
3514                         zero_var(result);
3515                         result->dscale = rscale;
3516                         return;
3517                 }
3518                 /* force maxdigits odd so that input ndigits can be equal */
3519                 if ((maxdigits & 1) == 0)
3520                         maxdigits++;
3521                 if (var1ndigits > var2ndigits)
3522                 {
3523                         var1ndigits -= res_ndigits - maxdigits;
3524                         if (var1ndigits < var2ndigits)
3525                                 var1ndigits = var2ndigits = (var1ndigits + var2ndigits) / 2;
3526                 }
3527                 else
3528                 {
3529                         var2ndigits -= res_ndigits - maxdigits;
3530                         if (var2ndigits < var1ndigits)
3531                                 var1ndigits = var2ndigits = (var1ndigits + var2ndigits) / 2;
3532                 }
3533                 res_ndigits = maxdigits;
3534                 Assert(res_ndigits == var1ndigits + var2ndigits + 1);
3535         }
3536
3537         /*
3538          * We do the arithmetic in an array "dig[]" of signed int's.  Since
3539          * INT_MAX is noticeably larger than NBASE*NBASE, this gives us headroom
3540          * to avoid normalizing carries immediately.
3541          *
3542          * maxdig tracks the maximum possible value of any dig[] entry;
3543          * when this threatens to exceed INT_MAX, we take the time to propagate
3544          * carries.  To avoid overflow in maxdig itself, it actually represents
3545          * the max possible value divided by NBASE-1.
3546          */
3547         dig = (int *) palloc0(res_ndigits * sizeof(int));
3548         maxdig = 0;
3549
3550         ri = res_ndigits - 1;
3551         for (i1 = var1ndigits - 1; i1 >= 0; ri--, i1--)
3552         {
3553                 int             var1digit = var1digits[i1];
3554
3555                 if (var1digit == 0)
3556                         continue;
3557
3558                 /* Time to normalize? */
3559                 maxdig += var1digit;
3560                 if (maxdig > INT_MAX/(NBASE-1))
3561                 {
3562                         /* Yes, do it */
3563                         carry = 0;
3564                         for (i = res_ndigits-1; i >= 0; i--)
3565                         {
3566                                 newdig = dig[i] + carry;
3567                                 if (newdig >= NBASE)
3568                                 {
3569                                         carry = newdig/NBASE;
3570                                         newdig -= carry*NBASE;
3571                                 }
3572                                 else
3573                                         carry = 0;
3574                                 dig[i] = newdig;
3575                         }
3576                         Assert(carry == 0);
3577                         /* Reset maxdig to indicate new worst-case */
3578                         maxdig = 1 + var1digit;
3579                 }
3580
3581                 /* Add appropriate multiple of var2 into the accumulator */
3582                 i = ri;
3583                 for (i2 = var2ndigits - 1; i2 >= 0; i2--)
3584                 {
3585                         dig[i--] += var1digit * var2digits[i2];
3586                 }
3587         }
3588
3589         /*
3590          * Now we do a final carry propagation pass to normalize the result,
3591          * which we combine with storing the result digits into the output.
3592          * Note that this is still done at full precision w/guard digits.
3593          */
3594         alloc_var(result, res_ndigits);
3595         res_digits = result->digits;
3596         carry = 0;
3597         for (i = res_ndigits-1; i >= 0; i--)
3598         {
3599                 newdig = dig[i] + carry;
3600                 if (newdig >= NBASE)
3601                 {
3602                         carry = newdig/NBASE;
3603                         newdig -= carry*NBASE;
3604                 }
3605                 else
3606                         carry = 0;
3607                 res_digits[i] = newdig;
3608         }
3609         Assert(carry == 0);
3610
3611         pfree(dig);
3612
3613         /*
3614          * Finally, round the result to the requested precision.
3615          */
3616         result->weight = res_weight;
3617         result->sign = res_sign;
3618
3619         /* Round to target rscale (and set result->dscale) */
3620         round_var(result, rscale);
3621
3622         /* Strip leading and trailing zeroes */
3623         strip_var(result);
3624 }
3625
3626
3627 /*
3628  * div_var() -
3629  *
3630  *      Division on variable level. Quotient of var1 / var2 is stored
3631  *      in result.  Result is rounded to no more than rscale fractional digits.
3632  */
3633 static void
3634 div_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
3635                 int rscale)
3636 {
3637         int                     div_ndigits;
3638         int                     res_sign;
3639         int                     res_weight;
3640         int                *div;
3641         int                     qdigit;
3642         int                     carry;
3643         int                     maxdiv;
3644         int                     newdig;
3645         NumericDigit *res_digits;
3646         double          fdividend,
3647                                 fdivisor,
3648                                 fdivisorinverse,
3649                                 fquotient;
3650         int                     qi;
3651         int                     i;
3652         /* copy these values into local vars for speed in inner loop */
3653         int                     var1ndigits = var1->ndigits;
3654         int                     var2ndigits = var2->ndigits;
3655         NumericDigit *var1digits = var1->digits;
3656         NumericDigit *var2digits = var2->digits;
3657
3658         /*
3659          * First of all division by zero check; we must not be handed an
3660          * unnormalized divisor.
3661          */
3662         if (var2ndigits == 0 || var2digits[0] == 0)
3663                 ereport(ERROR,
3664                                 (errcode(ERRCODE_DIVISION_BY_ZERO),
3665                                  errmsg("division by zero")));
3666
3667         /*
3668          * Now result zero check
3669          */
3670         if (var1ndigits == 0)
3671         {
3672                 zero_var(result);
3673                 result->dscale = rscale;
3674                 return;
3675         }
3676
3677         /*
3678          * Determine the result sign, weight and number of digits to calculate
3679          */
3680         if (var1->sign == var2->sign)
3681                 res_sign = NUMERIC_POS;
3682         else
3683                 res_sign = NUMERIC_NEG;
3684         res_weight = var1->weight - var2->weight + 1;
3685         /* The number of accurate result digits we need to produce: */
3686         div_ndigits = res_weight + 1 + (rscale + DEC_DIGITS-1)/DEC_DIGITS;
3687         /* Add guard digits for roundoff error */
3688         div_ndigits += DIV_GUARD_DIGITS;
3689         if (div_ndigits < DIV_GUARD_DIGITS)
3690                 div_ndigits = DIV_GUARD_DIGITS;
3691         /* Must be at least var1ndigits, too, to simplify data-loading loop */
3692         if (div_ndigits < var1ndigits)
3693                 div_ndigits = var1ndigits;
3694
3695         /*
3696          * We do the arithmetic in an array "div[]" of signed int's.  Since
3697          * INT_MAX is noticeably larger than NBASE*NBASE, this gives us headroom
3698          * to avoid normalizing carries immediately.
3699          *
3700          * We start with div[] containing one zero digit followed by the
3701          * dividend's digits (plus appended zeroes to reach the desired
3702          * precision including guard digits).  Each step of the main loop
3703          * computes an (approximate) quotient digit and stores it into div[],
3704          * removing one position of dividend space.  A final pass of carry
3705          * propagation takes care of any mistaken quotient digits.
3706          */
3707         div = (int *) palloc0((div_ndigits + 1) * sizeof(int));
3708         for (i = 0; i < var1ndigits; i++)
3709                 div[i+1] = var1digits[i];
3710
3711         /*
3712          * We estimate each quotient digit using floating-point arithmetic,
3713          * taking the first four digits of the (current) dividend and divisor.
3714          * This must be float to avoid overflow.
3715          */
3716         fdivisor = (double) var2digits[0];
3717         for (i = 1; i < 4; i++)
3718         {
3719                 fdivisor *= NBASE;
3720                 if (i < var2ndigits)
3721                         fdivisor += (double) var2digits[i];
3722         }
3723         fdivisorinverse = 1.0 / fdivisor;
3724
3725         /*
3726          * maxdiv tracks the maximum possible absolute value of any div[] entry;
3727          * when this threatens to exceed INT_MAX, we take the time to propagate
3728          * carries.  To avoid overflow in maxdiv itself, it actually represents
3729          * the max possible abs. value divided by NBASE-1.
3730          */
3731         maxdiv = 1;
3732
3733         /*
3734          * Outer loop computes next quotient digit, which will go into div[qi]
3735          */
3736         for (qi = 0; qi < div_ndigits; qi++)
3737         {
3738                 /* Approximate the current dividend value */
3739                 fdividend = (double) div[qi];
3740                 for (i = 1; i < 4; i++)
3741                 {
3742                         fdividend *= NBASE;
3743                         if (qi+i <= div_ndigits)
3744                                 fdividend += (double) div[qi+i];
3745                 }
3746                 /* Compute the (approximate) quotient digit */
3747                 fquotient = fdividend * fdivisorinverse;
3748                 qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
3749                         (((int) fquotient) - 1); /* truncate towards -infinity */
3750
3751                 if (qdigit != 0)
3752                 {
3753                         /* Do we need to normalize now? */
3754                         maxdiv += Abs(qdigit);
3755                         if (maxdiv > INT_MAX/(NBASE-1))
3756                         {
3757                                 /* Yes, do it */
3758                                 carry = 0;
3759                                 for (i = div_ndigits; i > qi; i--)
3760                                 {
3761                                         newdig = div[i] + carry;
3762                                         if (newdig < 0)
3763                                         {
3764                                                 carry = -((-newdig-1)/NBASE) - 1;
3765                                                 newdig -= carry*NBASE;
3766                                         }
3767                                         else if (newdig >= NBASE)
3768                                         {
3769                                                 carry = newdig/NBASE;
3770                                                 newdig -= carry*NBASE;
3771                                         }
3772                                         else
3773                                                 carry = 0;
3774                                         div[i] = newdig;
3775                                 }
3776                                 newdig = div[qi] + carry;
3777                                 div[qi] = newdig;
3778                                 /*
3779                                  * All the div[] digits except possibly div[qi] are now
3780                                  * in the range 0..NBASE-1.
3781                                  */
3782                                 maxdiv = Abs(newdig) / (NBASE-1);
3783                                 maxdiv = Max(maxdiv, 1);
3784                                 /*
3785                                  * Recompute the quotient digit since new info may have
3786                                  * propagated into the top four dividend digits
3787                                  */
3788                                 fdividend = (double) div[qi];
3789                                 for (i = 1; i < 4; i++)
3790                                 {
3791                                         fdividend *= NBASE;
3792                                         if (qi+i <= div_ndigits)
3793                                                 fdividend += (double) div[qi+i];
3794                                 }
3795                                 /* Compute the (approximate) quotient digit */
3796                                 fquotient = fdividend * fdivisorinverse;
3797                                 qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
3798                                         (((int) fquotient) - 1); /* truncate towards -infinity */
3799                                 maxdiv += Abs(qdigit);
3800                         }
3801
3802                         /* Subtract off the appropriate multiple of the divisor */
3803                         if (qdigit != 0)
3804                         {
3805                                 int             istop = Min(var2ndigits, div_ndigits-qi+1);
3806
3807                                 for (i = 0; i < istop; i++)
3808                                         div[qi+i] -= qdigit * var2digits[i];
3809                         }
3810                 }
3811
3812                 /*
3813                  * The dividend digit we are about to replace might still be nonzero.
3814                  * Fold it into the next digit position.  We don't need to worry about
3815                  * overflow here since this should nearly cancel with the subtraction
3816                  * of the divisor.
3817                  */
3818                 div[qi+1] += div[qi] * NBASE;
3819
3820                 div[qi] = qdigit;
3821         }
3822
3823         /*
3824          * Approximate and store the last quotient digit (div[div_ndigits])
3825          */
3826         fdividend = (double) div[qi];
3827         for (i = 1; i < 4; i++)
3828         {
3829                 fdividend *= NBASE;
3830         }
3831         fquotient = fdividend * fdivisorinverse;
3832         qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
3833                 (((int) fquotient) - 1); /* truncate towards -infinity */
3834         div[qi] = qdigit;
3835
3836         /*
3837          * Now we do a final carry propagation pass to normalize the result,
3838          * which we combine with storing the result digits into the output.
3839          * Note that this is still done at full precision w/guard digits.
3840          */
3841         alloc_var(result, div_ndigits+1);
3842         res_digits = result->digits;
3843         carry = 0;
3844         for (i = div_ndigits; i >= 0; i--)
3845         {
3846                 newdig = div[i] + carry;
3847                 if (newdig < 0)
3848                 {
3849                         carry = -((-newdig-1)/NBASE) - 1;
3850                         newdig -= carry*NBASE;
3851                 }
3852                 else if (newdig >= NBASE)
3853                 {
3854                         carry = newdig/NBASE;
3855                         newdig -= carry*NBASE;
3856                 }
3857                 else
3858                         carry = 0;
3859                 res_digits[i] = newdig;
3860         }
3861         Assert(carry == 0);
3862
3863         pfree(div);
3864
3865         /*
3866          * Finally, round the result to the requested precision.
3867          */
3868         result->weight = res_weight;
3869         result->sign = res_sign;
3870
3871         /* Round to target rscale (and set result->dscale) */
3872         round_var(result, rscale);
3873
3874         /* Strip leading and trailing zeroes */
3875         strip_var(result);
3876 }
3877
3878
3879 /*
3880  * Default scale selection for division
3881  *
3882  * Returns the appropriate result scale for the division result.
3883  */
3884 static int
3885 select_div_scale(NumericVar *var1, NumericVar *var2)
3886 {
3887         int                     weight1,
3888                                 weight2,
3889                                 qweight,
3890                                 i;
3891         NumericDigit firstdigit1,
3892                                 firstdigit2;
3893         int                     rscale;
3894
3895         /*
3896          * The result scale of a division isn't specified in any SQL standard.
3897          * For PostgreSQL we select a result scale that will give at least
3898          * NUMERIC_MIN_SIG_DIGITS significant digits, so that numeric gives a
3899          * result no less accurate than float8; but use a scale not less than
3900          * either input's display scale.
3901          */
3902
3903         /* Get the actual (normalized) weight and first digit of each input */
3904
3905         weight1 = 0;                            /* values to use if var1 is zero */
3906         firstdigit1 = 0;
3907         for (i = 0; i < var1->ndigits; i++)
3908         {
3909                 firstdigit1 = var1->digits[i];
3910                 if (firstdigit1 != 0)
3911                 {
3912                         weight1 = var1->weight - i;
3913                         break;
3914                 }
3915         }
3916
3917         weight2 = 0;                            /* values to use if var2 is zero */
3918         firstdigit2 = 0;
3919         for (i = 0; i < var2->ndigits; i++)
3920         {
3921                 firstdigit2 = var2->digits[i];
3922                 if (firstdigit2 != 0)
3923                 {
3924                         weight2 = var2->weight - i;
3925                         break;
3926                 }
3927         }
3928
3929         /*
3930          * Estimate weight of quotient.  If the two first digits are equal,
3931          * we can't be sure, but assume that var1 is less than var2.
3932          */
3933         qweight = weight1 - weight2;
3934         if (firstdigit1 <= firstdigit2)
3935                 qweight--;
3936
3937         /* Select result scale */
3938         rscale = NUMERIC_MIN_SIG_DIGITS - qweight * DEC_DIGITS;
3939         rscale = Max(rscale, var1->dscale);
3940         rscale = Max(rscale, var2->dscale);
3941         rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
3942         rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
3943
3944         return rscale;
3945 }
3946
3947
3948 /*
3949  * mod_var() -
3950  *
3951  *      Calculate the modulo of two numerics at variable level
3952  */
3953 static void
3954 mod_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
3955 {
3956         NumericVar      tmp;
3957         int                     rscale;
3958
3959         init_var(&tmp);
3960
3961         /* ---------
3962          * We do this using the equation
3963          *              mod(x,y) = x - trunc(x/y)*y
3964          * We set rscale the same way numeric_div and numeric_mul do
3965          * to get the right answer from the equation.  The final result,
3966          * however, need not be displayed to more precision than the inputs.
3967          * ----------
3968          */
3969         rscale = select_div_scale(var1, var2);
3970
3971         div_var(var1, var2, &tmp, rscale);
3972
3973         trunc_var(&tmp, 0);
3974
3975         mul_var(var2, &tmp, &tmp, var2->dscale + tmp.dscale);
3976
3977         sub_var(var1, &tmp, result);
3978
3979         round_var(result, Max(var1->dscale, var2->dscale));
3980
3981         free_var(&tmp);
3982 }
3983
3984
3985 /*
3986  * ceil_var() -
3987  *
3988  *      Return the smallest integer greater than or equal to the argument
3989  *      on variable level
3990  */
3991 static void
3992 ceil_var(NumericVar *var, NumericVar *result)
3993 {
3994         NumericVar      tmp;
3995
3996         init_var(&tmp);
3997         set_var_from_var(var, &tmp);
3998
3999         trunc_var(&tmp, 0);
4000
4001         if (var->sign == NUMERIC_POS && cmp_var(var, &tmp) != 0)
4002                 add_var(&tmp, &const_one, &tmp);
4003
4004         set_var_from_var(&tmp, result);
4005         free_var(&tmp);
4006 }
4007
4008
4009 /*
4010  * floor_var() -
4011  *
4012  *      Return the largest integer equal to or less than the argument
4013  *      on variable level
4014  */
4015 static void
4016 floor_var(NumericVar *var, NumericVar *result)
4017 {
4018         NumericVar      tmp;
4019
4020         init_var(&tmp);
4021         set_var_from_var(var, &tmp);
4022
4023         trunc_var(&tmp, 0);
4024
4025         if (var->sign == NUMERIC_NEG && cmp_var(var, &tmp) != 0)
4026                 sub_var(&tmp, &const_one, &tmp);
4027
4028         set_var_from_var(&tmp, result);
4029         free_var(&tmp);
4030 }
4031
4032
4033 /*
4034  * sqrt_var() -
4035  *
4036  *      Compute the square root of x using Newton's algorithm
4037  */
4038 static void
4039 sqrt_var(NumericVar *arg, NumericVar *result, int rscale)
4040 {
4041         NumericVar      tmp_arg;
4042         NumericVar      tmp_val;
4043         NumericVar      last_val;
4044         int                     local_rscale;
4045         int                     stat;
4046
4047         local_rscale = rscale + 8;
4048
4049         stat = cmp_var(arg, &const_zero);
4050         if (stat == 0)
4051         {
4052                 zero_var(result);
4053                 result->dscale = rscale;
4054                 return;
4055         }
4056
4057         if (stat < 0)
4058                 ereport(ERROR,
4059                                 (errcode(ERRCODE_FLOATING_POINT_EXCEPTION),
4060                                  errmsg("cannot take square root of a negative number")));
4061
4062         init_var(&tmp_arg);
4063         init_var(&tmp_val);
4064         init_var(&last_val);
4065
4066         /* Copy arg in case it is the same var as result */
4067         set_var_from_var(arg, &tmp_arg);
4068
4069         /*
4070          * Initialize the result to the first guess
4071          */
4072         alloc_var(result, 1);
4073         result->digits[0] = tmp_arg.digits[0] / 2;
4074         if (result->digits[0] == 0)
4075                 result->digits[0] = 1;
4076         result->weight = tmp_arg.weight / 2;
4077         result->sign = NUMERIC_POS;
4078
4079         set_var_from_var(result, &last_val);
4080
4081         for (;;)
4082         {
4083                 div_var(&tmp_arg, result, &tmp_val, local_rscale);
4084
4085                 add_var(result, &tmp_val, result);
4086                 mul_var(result, &const_zero_point_five, result, local_rscale);
4087
4088                 if (cmp_var(&last_val, result) == 0)
4089                         break;
4090                 set_var_from_var(result, &last_val);
4091         }
4092
4093         free_var(&last_val);
4094         free_var(&tmp_val);
4095         free_var(&tmp_arg);
4096
4097         /* Round to requested precision */
4098         round_var(result, rscale);
4099 }
4100
4101
4102 /*
4103  * exp_var() -
4104  *
4105  *      Raise e to the power of x
4106  */
4107 static void
4108 exp_var(NumericVar *arg, NumericVar *result, int rscale)
4109 {
4110         NumericVar      x;
4111         int                     xintval;
4112         bool            xneg = FALSE;
4113         int                     local_rscale;
4114
4115         /*----------
4116          * We separate the integral and fraction parts of x, then compute
4117          *              e^x = e^xint * e^xfrac
4118          * where e = exp(1) and e^xfrac = exp(xfrac) are computed by
4119          * exp_var_internal; the limited range of inputs allows that routine
4120          * to do a good job with a simple Taylor series.  Raising e^xint is
4121          * done by repeated multiplications in power_var_int.
4122          *----------
4123          */
4124         init_var(&x);
4125
4126         set_var_from_var(arg, &x);
4127
4128         if (x.sign == NUMERIC_NEG)
4129         {
4130                 xneg = TRUE;
4131                 x.sign = NUMERIC_POS;
4132         }
4133
4134         /* Extract the integer part, remove it from x */
4135         xintval = 0;
4136         while (x.weight >= 0)
4137         {
4138                 xintval *= NBASE;
4139                 if (x.ndigits > 0)
4140                 {
4141                         xintval += x.digits[0];
4142                         x.digits++;
4143                         x.ndigits--;
4144                 }
4145                 x.weight--;
4146                 /* Guard against overflow */
4147                 if (xintval >= NUMERIC_MAX_RESULT_SCALE * 3)
4148                         ereport(ERROR,
4149                                         (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
4150                                          errmsg("argument for EXP() too big")));
4151         }
4152
4153         /* Select an appropriate scale for internal calculation */
4154         local_rscale = rscale + MUL_GUARD_DIGITS * 2;
4155
4156         /* Compute e^xfrac */
4157         exp_var_internal(&x, result, local_rscale);
4158
4159         /* If there's an integer part, multiply by e^xint */
4160         if (xintval > 0)
4161         {
4162                 NumericVar      e;
4163
4164                 init_var(&e);
4165                 exp_var_internal(&const_one, &e, local_rscale);
4166                 power_var_int(&e, xintval, &e, local_rscale);
4167                 mul_var(&e, result, result, local_rscale);
4168                 free_var(&e);
4169         }
4170
4171         /* Compensate for input sign, and round to requested rscale */
4172         if (xneg)
4173                 div_var(&const_one, result, result, rscale);
4174         else
4175                 round_var(result, rscale);
4176
4177         free_var(&x);
4178 }
4179
4180
4181 /*
4182  * exp_var_internal() -
4183  *
4184  *      Raise e to the power of x, where 0 <= x <= 1
4185  *
4186  * NB: the result should be good to at least rscale digits, but it has
4187  * *not* been rounded off; the caller must do that if wanted.
4188  */
4189 static void
4190 exp_var_internal(NumericVar *arg, NumericVar *result, int rscale)
4191 {
4192         NumericVar      x;
4193         NumericVar      xpow;
4194         NumericVar      ifac;
4195         NumericVar      elem;
4196         NumericVar      ni;
4197         int                     ndiv2 = 0;
4198         int                     local_rscale;
4199
4200         init_var(&x);
4201         init_var(&xpow);
4202         init_var(&ifac);
4203         init_var(&elem);
4204         init_var(&ni);
4205
4206         set_var_from_var(arg, &x);
4207
4208         Assert(x.sign == NUMERIC_POS);
4209
4210         local_rscale = rscale + 8;
4211
4212         /* Reduce input into range 0 <= x <= 0.01 */
4213         while (cmp_var(&x, &const_zero_point_01) > 0)
4214         {
4215                 ndiv2++;
4216                 local_rscale++;
4217                 mul_var(&x, &const_zero_point_five, &x, x.dscale+1);
4218         }
4219
4220         /*
4221          * Use the Taylor series
4222          *
4223          *              exp(x) = 1 + x + x^2/2! + x^3/3! + ...
4224          *
4225          * Given the limited range of x, this should converge reasonably quickly.
4226          * We run the series until the terms fall below the local_rscale limit.
4227          */
4228         add_var(&const_one, &x, result);
4229         set_var_from_var(&x, &xpow);
4230         set_var_from_var(&const_one, &ifac);
4231         set_var_from_var(&const_one, &ni);
4232
4233         for (;;)
4234         {
4235                 add_var(&ni, &const_one, &ni);
4236                 mul_var(&xpow, &x, &xpow, local_rscale);
4237                 mul_var(&ifac, &ni, &ifac, 0);
4238                 div_var(&xpow, &ifac, &elem, local_rscale);
4239
4240                 if (elem.ndigits == 0)
4241                         break;
4242
4243                 add_var(result, &elem, result);
4244         }
4245
4246         /* Compensate for argument range reduction */
4247         while (ndiv2-- > 0)
4248                 mul_var(result, result, result, local_rscale);
4249
4250         free_var(&x);
4251         free_var(&xpow);
4252         free_var(&ifac);
4253         free_var(&elem);
4254         free_var(&ni);
4255 }
4256
4257
4258 /*
4259  * ln_var() -
4260  *
4261  *      Compute the natural log of x
4262  */
4263 static void
4264 ln_var(NumericVar *arg, NumericVar *result, int rscale)
4265 {
4266         NumericVar      x;
4267         NumericVar      xx;
4268         NumericVar      ni;
4269         NumericVar      elem;
4270         NumericVar      fact;
4271         int                     local_rscale;
4272
4273         if (cmp_var(arg, &const_zero) <= 0)
4274                 ereport(ERROR,
4275                                 (errcode(ERRCODE_FLOATING_POINT_EXCEPTION),
4276                                  errmsg("cannot take log of a negative number")));
4277
4278         local_rscale = rscale + 8;
4279
4280         init_var(&x);
4281         init_var(&xx);
4282         init_var(&ni);
4283         init_var(&elem);
4284         init_var(&fact);
4285
4286         set_var_from_var(arg, &x);
4287         set_var_from_var(&const_two, &fact);
4288
4289         /* Reduce input into range 0.9 < x < 1.1 */
4290         while (cmp_var(&x, &const_zero_point_nine) <= 0)
4291         {
4292                 local_rscale++;
4293                 sqrt_var(&x, &x, local_rscale);
4294                 mul_var(&fact, &const_two, &fact, 0);
4295         }
4296         while (cmp_var(&x, &const_one_point_one) >= 0)
4297         {
4298                 local_rscale++;
4299                 sqrt_var(&x, &x, local_rscale);
4300                 mul_var(&fact, &const_two, &fact, 0);
4301         }
4302
4303         /*
4304          * We use the Taylor series for 0.5 * ln((1+z)/(1-z)),
4305          *
4306          *              z + z^3/3 + z^5/5 + ...
4307          *
4308          * where z = (x-1)/(x+1) is in the range (approximately) -0.053 .. 0.048
4309          * due to the above range-reduction of x.
4310          *
4311          * The convergence of this is not as fast as one would like, but is
4312          * tolerable given that z is small.
4313          */
4314         sub_var(&x, &const_one, result);
4315         add_var(&x, &const_one, &elem);
4316         div_var(result, &elem, result, local_rscale);
4317         set_var_from_var(result, &xx);
4318         mul_var(result, result, &x, local_rscale);
4319
4320         set_var_from_var(&const_one, &ni);
4321
4322         for (;;)
4323         {
4324                 add_var(&ni, &const_two, &ni);
4325                 mul_var(&xx, &x, &xx, local_rscale);
4326                 div_var(&xx, &ni, &elem, local_rscale);
4327
4328                 if (elem.ndigits == 0)
4329                         break;
4330
4331                 add_var(result, &elem, result);
4332
4333                 if (elem.weight < (result->weight - local_rscale * 2/DEC_DIGITS))
4334                         break;
4335         }
4336
4337         /* Compensate for argument range reduction, round to requested rscale */
4338         mul_var(result, &fact, result, rscale);
4339
4340         free_var(&x);
4341         free_var(&xx);
4342         free_var(&ni);
4343         free_var(&elem);
4344         free_var(&fact);
4345 }
4346
4347
4348 /*
4349  * log_var() -
4350  *
4351  *      Compute the logarithm of num in a given base.
4352  *
4353  *      Note: this routine chooses dscale of the result.
4354  */
4355 static void
4356 log_var(NumericVar *base, NumericVar *num, NumericVar *result)
4357 {
4358         NumericVar      ln_base;
4359         NumericVar      ln_num;
4360         int                     dec_digits;
4361         int                     rscale;
4362         int                     local_rscale;
4363
4364         init_var(&ln_base);
4365         init_var(&ln_num);
4366
4367         /* Set scale for ln() calculations --- compare numeric_ln() */
4368
4369         /* Approx decimal digits before decimal point */
4370         dec_digits = (num->weight + 1) * DEC_DIGITS;
4371
4372         if (dec_digits > 1)
4373                 rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(dec_digits - 1);
4374         else if (dec_digits < 1)
4375                 rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(1 - dec_digits);
4376         else
4377                 rscale = NUMERIC_MIN_SIG_DIGITS;
4378
4379         rscale = Max(rscale, base->dscale);
4380         rscale = Max(rscale, num->dscale);
4381         rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
4382         rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
4383
4384         local_rscale = rscale + 8;
4385
4386         /* Form natural logarithms */
4387         ln_var(base, &ln_base, local_rscale);
4388         ln_var(num, &ln_num, local_rscale);
4389
4390         ln_base.dscale = rscale;
4391         ln_num.dscale = rscale;
4392
4393         /* Select scale for division result */
4394         rscale = select_div_scale(&ln_num, &ln_base);
4395
4396         div_var(&ln_num, &ln_base, result, rscale);
4397
4398         free_var(&ln_num);
4399         free_var(&ln_base);
4400 }
4401
4402
4403 /*
4404  * power_var() -
4405  *
4406  *      Raise base to the power of exp
4407  *
4408  *      Note: this routine chooses dscale of the result.
4409  */
4410 static void
4411 power_var(NumericVar *base, NumericVar *exp, NumericVar *result)
4412 {
4413         NumericVar      ln_base;
4414         NumericVar      ln_num;
4415         int                     dec_digits;
4416         int                     rscale;
4417         int                     local_rscale;
4418         double          val;
4419
4420         /* If exp can be represented as an integer, use power_var_int */
4421         if (exp->ndigits == 0 || exp->ndigits <= exp->weight + 1)
4422         {
4423                 /* exact integer, but does it fit in int? */
4424                 NumericVar      x;
4425                 int64           expval64;
4426
4427                 /* must copy because numericvar_to_int8() scribbles on input */
4428                 init_var(&x);
4429                 set_var_from_var(exp, &x);
4430                 if (numericvar_to_int8(&x, &expval64))
4431                 {
4432                         int             expval = (int) expval64;
4433
4434                         /* Test for overflow by reverse-conversion. */
4435                         if ((int64) expval == expval64)
4436                         {
4437                                 /* Okay, select rscale */
4438                                 rscale = NUMERIC_MIN_SIG_DIGITS;
4439                                 rscale = Max(rscale, base->dscale);
4440                                 rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
4441                                 rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
4442
4443                                 power_var_int(base, expval, result, rscale);
4444
4445                                 free_var(&x);
4446                                 return;
4447                         }
4448                 }
4449                 free_var(&x);
4450         }
4451
4452         init_var(&ln_base);
4453         init_var(&ln_num);
4454
4455         /* Set scale for ln() calculation --- need extra accuracy here */
4456
4457         /* Approx decimal digits before decimal point */
4458         dec_digits = (base->weight + 1) * DEC_DIGITS;
4459
4460         if (dec_digits > 1)
4461                 rscale = NUMERIC_MIN_SIG_DIGITS*2 - (int) log10(dec_digits - 1);
4462         else if (dec_digits < 1)
4463                 rscale = NUMERIC_MIN_SIG_DIGITS*2 - (int) log10(1 - dec_digits);
4464         else
4465                 rscale = NUMERIC_MIN_SIG_DIGITS*2;
4466
4467         rscale = Max(rscale, base->dscale * 2);
4468         rscale = Max(rscale, exp->dscale * 2);
4469         rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE * 2);
4470         rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE * 2);
4471
4472         local_rscale = rscale + 8;
4473
4474         ln_var(base, &ln_base, local_rscale);
4475
4476         mul_var(&ln_base, exp, &ln_num, local_rscale);
4477
4478         /* Set scale for exp() -- compare numeric_exp() */
4479
4480         /* convert input to float8, ignoring overflow */
4481         val = numericvar_to_double_no_overflow(&ln_num);
4482
4483         /* log10(result) = num * log10(e), so this is approximately the weight: */
4484         val *= 0.434294481903252;
4485
4486         /* limit to something that won't cause integer overflow */
4487         val = Max(val, -NUMERIC_MAX_RESULT_SCALE);
4488         val = Min(val, NUMERIC_MAX_RESULT_SCALE);
4489
4490         rscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
4491         rscale = Max(rscale, base->dscale);
4492         rscale = Max(rscale, exp->dscale);
4493         rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
4494         rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
4495
4496         exp_var(&ln_num, result, rscale);
4497
4498         free_var(&ln_num);
4499         free_var(&ln_base);
4500 }
4501
4502 /*
4503  * power_var_int() -
4504  *
4505  *      Raise base to the power of exp, where exp is an integer.
4506  */
4507 static void
4508 power_var_int(NumericVar *base, int exp, NumericVar *result, int rscale)
4509 {
4510         bool            neg;
4511         NumericVar      base_prod;
4512         int                     local_rscale;
4513
4514         /* Detect some special cases, particularly 0^0. */
4515
4516         switch (exp)
4517         {
4518                 case 0:
4519                         if (base->ndigits == 0)
4520                                 ereport(ERROR,
4521                                                 (errcode(ERRCODE_FLOATING_POINT_EXCEPTION),
4522                                                  errmsg("zero raised to zero is undefined")));
4523                         set_var_from_var(&const_one, result);
4524                         result->dscale = rscale; /* no need to round */
4525                         return;
4526                 case 1:
4527                         set_var_from_var(base, result);
4528                         round_var(result, rscale);
4529                         return;
4530                 case -1:
4531                         div_var(&const_one, base, result, rscale);
4532                         return;
4533                 case 2:
4534                         mul_var(base, base, result, rscale);
4535                         return;
4536                 default:
4537                         break;
4538         }
4539
4540         /*
4541          * The general case repeatedly multiplies base according to the
4542          * bit pattern of exp.  We do the multiplications with some extra
4543          * precision.
4544          */
4545         neg = (exp < 0);
4546         exp = Abs(exp);
4547
4548         local_rscale = rscale + MUL_GUARD_DIGITS * 2;
4549
4550         init_var(&base_prod);
4551         set_var_from_var(base, &base_prod);
4552
4553         if (exp & 1)
4554                 set_var_from_var(base, result);
4555         else
4556                 set_var_from_var(&const_one, result);
4557
4558         while ((exp >>= 1) > 0)
4559         {
4560                 mul_var(&base_prod, &base_prod, &base_prod, local_rscale);
4561                 if (exp & 1)
4562                         mul_var(&base_prod, result, result, local_rscale);
4563         }
4564
4565         free_var(&base_prod);
4566
4567         /* Compensate for input sign, and round to requested rscale */
4568         if (neg)
4569                 div_var(&const_one, result, result, rscale);
4570         else
4571                 round_var(result, rscale);
4572 }
4573
4574
4575 /* ----------------------------------------------------------------------
4576  *
4577  * Following are the lowest level functions that operate unsigned
4578  * on the variable level
4579  *
4580  * ----------------------------------------------------------------------
4581  */
4582
4583
4584 /* ----------
4585  * cmp_abs() -
4586  *
4587  *      Compare the absolute values of var1 and var2
4588  *      Returns:        -1 for ABS(var1) < ABS(var2)
4589  *                              0  for ABS(var1) == ABS(var2)
4590  *                              1  for ABS(var1) > ABS(var2)
4591  * ----------
4592  */
4593 static int
4594 cmp_abs(NumericVar *var1, NumericVar *var2)
4595 {
4596         NumericDigit *var1digits = var1->digits;
4597         NumericDigit *var2digits = var2->digits;
4598         int                     i1 = 0;
4599         int                     i2 = 0;
4600         int                     w1 = var1->weight;
4601         int                     w2 = var2->weight;
4602
4603         /* Check any digits before the first common digit */
4604
4605         while (w1 > w2 && i1 < var1->ndigits)
4606         {
4607                 if (var1digits[i1++] != 0)
4608                         return 1;
4609                 w1--;
4610         }
4611         while (w2 > w1 && i2 < var2->ndigits)
4612         {
4613                 if (var2digits[i2++] != 0)
4614                         return -1;
4615                 w2--;
4616         }
4617
4618         /* At this point, either w1 == w2 or we've run out of digits */
4619
4620         if (w1 == w2)
4621         {
4622                 while (i1 < var1->ndigits && i2 < var2->ndigits)
4623                 {
4624                         int                     stat = var1digits[i1++] - var2digits[i2++];
4625
4626                         if (stat)
4627                         {
4628                                 if (stat > 0)
4629                                         return 1;
4630                                 return -1;
4631                         }
4632                 }
4633         }
4634
4635         /*
4636          * At this point, we've run out of digits on one side or the other;
4637          * so any remaining nonzero digits imply that side is larger
4638          */
4639         while (i1 < var1->ndigits)
4640         {
4641                 if (var1digits[i1++] != 0)
4642                         return 1;
4643         }
4644         while (i2 < var2->ndigits)
4645         {
4646                 if (var2digits[i2++] != 0)
4647                         return -1;
4648         }
4649
4650         return 0;
4651 }
4652
4653
4654 /*
4655  * add_abs() -
4656  *
4657  *      Add the absolute values of two variables into result.
4658  *      result might point to one of the operands without danger.
4659  */
4660 static void
4661 add_abs(NumericVar *var1, NumericVar *var2, NumericVar *result)
4662 {
4663         NumericDigit *res_buf;
4664         NumericDigit *res_digits;
4665         int                     res_ndigits;
4666         int                     res_weight;
4667         int                     res_rscale,
4668                                 rscale1,
4669                                 rscale2;
4670         int                     res_dscale;
4671         int                     i,
4672                                 i1,
4673                                 i2;
4674         int                     carry = 0;
4675
4676         /* copy these values into local vars for speed in inner loop */
4677         int                     var1ndigits = var1->ndigits;
4678         int                     var2ndigits = var2->ndigits;
4679         NumericDigit *var1digits = var1->digits;
4680         NumericDigit *var2digits = var2->digits;
4681
4682         res_weight = Max(var1->weight, var2->weight) + 1;
4683
4684         res_dscale = Max(var1->dscale, var2->dscale);
4685
4686         /* Note: here we are figuring rscale in base-NBASE digits */
4687         rscale1 = var1->ndigits - var1->weight - 1;
4688         rscale2 = var2->ndigits - var2->weight - 1;
4689         res_rscale = Max(rscale1, rscale2);
4690
4691         res_ndigits = res_rscale + res_weight + 1;
4692         if (res_ndigits <= 0)
4693                 res_ndigits = 1;
4694
4695         res_buf = digitbuf_alloc(res_ndigits + 1);
4696         res_buf[0] = 0;                         /* spare digit for later rounding */
4697         res_digits = res_buf + 1;
4698
4699         i1 = res_rscale + var1->weight + 1;
4700         i2 = res_rscale + var2->weight + 1;
4701         for (i = res_ndigits - 1; i >= 0; i--)
4702         {
4703                 i1--;
4704                 i2--;
4705                 if (i1 >= 0 && i1 < var1ndigits)
4706                         carry += var1digits[i1];
4707                 if (i2 >= 0 && i2 < var2ndigits)
4708                         carry += var2digits[i2];
4709
4710                 if (carry >= NBASE)
4711                 {
4712                         res_digits[i] = carry - NBASE;
4713                         carry = 1;
4714                 }
4715                 else
4716                 {
4717                         res_digits[i] = carry;
4718                         carry = 0;
4719                 }
4720         }
4721
4722         Assert(carry == 0);                     /* else we failed to allow for carry out */
4723
4724         digitbuf_free(result->buf);
4725         result->ndigits = res_ndigits;
4726         result->buf = res_buf;
4727         result->digits = res_digits;
4728         result->weight = res_weight;
4729         result->dscale = res_dscale;
4730
4731         /* Remove leading/trailing zeroes */
4732         strip_var(result);
4733 }
4734
4735
4736 /*
4737  * sub_abs()
4738  *
4739  *      Subtract the absolute value of var2 from the absolute value of var1
4740  *      and store in result. result might point to one of the operands
4741  *      without danger.
4742  *
4743  *      ABS(var1) MUST BE GREATER OR EQUAL ABS(var2) !!!
4744  */
4745 static void
4746 sub_abs(NumericVar *var1, NumericVar *var2, NumericVar *result)
4747 {
4748         NumericDigit *res_buf;
4749         NumericDigit *res_digits;
4750         int                     res_ndigits;
4751         int                     res_weight;
4752         int                     res_rscale,
4753                                 rscale1,
4754                                 rscale2;
4755         int                     res_dscale;
4756         int                     i,
4757                                 i1,
4758                                 i2;
4759         int                     borrow = 0;
4760
4761         /* copy these values into local vars for speed in inner loop */
4762         int                     var1ndigits = var1->ndigits;
4763         int                     var2ndigits = var2->ndigits;
4764         NumericDigit *var1digits = var1->digits;
4765         NumericDigit *var2digits = var2->digits;
4766
4767         res_weight = var1->weight;
4768
4769         res_dscale = Max(var1->dscale, var2->dscale);
4770
4771         /* Note: here we are figuring rscale in base-NBASE digits */
4772         rscale1 = var1->ndigits - var1->weight - 1;
4773         rscale2 = var2->ndigits - var2->weight - 1;
4774         res_rscale = Max(rscale1, rscale2);
4775
4776         res_ndigits = res_rscale + res_weight + 1;
4777         if (res_ndigits <= 0)
4778                 res_ndigits = 1;
4779
4780         res_buf = digitbuf_alloc(res_ndigits + 1);
4781         res_buf[0] = 0;                         /* spare digit for later rounding */
4782         res_digits = res_buf + 1;
4783
4784         i1 = res_rscale + var1->weight + 1;
4785         i2 = res_rscale + var2->weight + 1;
4786         for (i = res_ndigits - 1; i >= 0; i--)
4787         {
4788                 i1--;
4789                 i2--;
4790                 if (i1 >= 0 && i1 < var1ndigits)
4791                         borrow += var1digits[i1];
4792                 if (i2 >= 0 && i2 < var2ndigits)
4793                         borrow -= var2digits[i2];
4794
4795                 if (borrow < 0)
4796                 {
4797                         res_digits[i] = borrow + NBASE;
4798                         borrow = -1;
4799                 }
4800                 else
4801                 {
4802                         res_digits[i] = borrow;
4803                         borrow = 0;
4804                 }
4805         }
4806
4807         Assert(borrow == 0);            /* else caller gave us var1 < var2 */
4808
4809         digitbuf_free(result->buf);
4810         result->ndigits = res_ndigits;
4811         result->buf = res_buf;
4812         result->digits = res_digits;
4813         result->weight = res_weight;
4814         result->dscale = res_dscale;
4815
4816         /* Remove leading/trailing zeroes */
4817         strip_var(result);
4818 }
4819
4820 /*
4821  * round_var
4822  *
4823  * Round the value of a variable to no more than rscale decimal digits
4824  * after the decimal point.  NOTE: we allow rscale < 0 here, implying
4825  * rounding before the decimal point.
4826  */
4827 static void
4828 round_var(NumericVar *var, int rscale)
4829 {
4830         NumericDigit   *digits = var->digits;
4831         int                     di;
4832         int                     ndigits;
4833         int                     carry;
4834
4835         var->dscale = rscale;
4836
4837         /* decimal digits wanted */
4838         di = (var->weight + 1) * DEC_DIGITS + rscale;
4839
4840         /*
4841          * If di = 0, the value loses all digits, but could round up to 1
4842          * if its first extra digit is >= 5.  If di < 0 the result must be 0.
4843          */
4844         if (di < 0)
4845         {
4846                 var->ndigits = 0;
4847                 var->weight = 0;
4848                 var->sign = NUMERIC_POS;
4849         }
4850         else
4851         {
4852                 /* NBASE digits wanted */
4853                 ndigits = (di + DEC_DIGITS-1) / DEC_DIGITS;
4854
4855                 /* 0, or number of decimal digits to keep in last NBASE digit */
4856                 di %= DEC_DIGITS;
4857
4858                 if (ndigits < var->ndigits ||
4859                         (ndigits == var->ndigits && di > 0))
4860                 {
4861                         var->ndigits = ndigits;
4862
4863 #if DEC_DIGITS == 1
4864                         /* di must be zero */
4865                         carry = (digits[ndigits] >= HALF_NBASE) ? 1 : 0;
4866 #else
4867                         if (di == 0)
4868                         {
4869                                 carry = (digits[ndigits] >= HALF_NBASE) ? 1 : 0;
4870                         }
4871                         else
4872                         {
4873                                 /* Must round within last NBASE digit */
4874                                 int             extra,
4875                                                 pow10;
4876
4877 #if DEC_DIGITS == 4
4878                                 pow10 = round_powers[di];
4879 #elif DEC_DIGITS == 2
4880                                 pow10 = 10;
4881 #else
4882 #error unsupported NBASE
4883 #endif
4884                                 extra = digits[--ndigits] % pow10;
4885                                 digits[ndigits] -= extra;
4886                                 carry = 0;
4887                                 if (extra >= pow10/2)
4888                                 {
4889                                         pow10 += digits[ndigits];
4890                                         if (pow10 >= NBASE)
4891                                         {
4892                                                 pow10 -= NBASE;
4893                                                 carry = 1;
4894                                         }
4895                                         digits[ndigits] = pow10;
4896                                 }
4897                         }
4898 #endif
4899
4900                         /* Propagate carry if needed */
4901                         while (carry)
4902                         {
4903                                 carry += digits[--ndigits];
4904                                 if (carry >= NBASE)
4905                                 {
4906                                         digits[ndigits] = carry - NBASE;
4907                                         carry = 1;
4908                                 }
4909                                 else
4910                                 {
4911                                         digits[ndigits] = carry;
4912                                         carry = 0;
4913                                 }
4914                         }
4915
4916                         if (ndigits < 0)
4917                         {
4918                                 Assert(ndigits == -1);  /* better not have added > 1 digit */
4919                                 Assert(var->digits > var->buf);
4920                                 var->digits--;
4921                                 var->ndigits++;
4922                                 var->weight++;
4923                         }
4924                 }
4925         }
4926 }
4927
4928 /*
4929  * trunc_var
4930  *
4931  * Truncate the value of a variable at rscale decimal digits after the
4932  * decimal point.  NOTE: we allow rscale < 0 here, implying
4933  * truncation before the decimal point.
4934  */
4935 static void
4936 trunc_var(NumericVar *var, int rscale)
4937 {
4938         int                     di;
4939         int                     ndigits;
4940
4941         var->dscale = rscale;
4942
4943         /* decimal digits wanted */
4944         di = (var->weight + 1) * DEC_DIGITS + rscale;
4945
4946         /*
4947          * If di <= 0, the value loses all digits.
4948          */
4949         if (di <= 0)
4950         {
4951                 var->ndigits = 0;
4952                 var->weight = 0;
4953                 var->sign = NUMERIC_POS;
4954         }
4955         else
4956         {
4957                 /* NBASE digits wanted */
4958                 ndigits = (di + DEC_DIGITS-1) / DEC_DIGITS;
4959
4960                 if (ndigits <= var->ndigits)
4961                 {
4962                         var->ndigits = ndigits;
4963
4964 #if DEC_DIGITS == 1
4965                         /* no within-digit stuff to worry about */
4966 #else
4967                         /* 0, or number of decimal digits to keep in last NBASE digit */
4968                         di %= DEC_DIGITS;
4969
4970                         if (di > 0)
4971                         {
4972                                 /* Must truncate within last NBASE digit */
4973                                 NumericDigit   *digits = var->digits;
4974                                 int             extra,
4975                                                 pow10;
4976
4977 #if DEC_DIGITS == 4
4978                                 pow10 = round_powers[di];
4979 #elif DEC_DIGITS == 2
4980                                 pow10 = 10;
4981 #else
4982 #error unsupported NBASE
4983 #endif
4984                                 extra = digits[--ndigits] % pow10;
4985                                 digits[ndigits] -= extra;
4986                         }
4987 #endif
4988                 }
4989         }
4990 }
4991
4992 /*
4993  * strip_var
4994  *
4995  * Strip any leading and trailing zeroes from a numeric variable
4996  */
4997 static void
4998 strip_var(NumericVar *var)
4999 {
5000         NumericDigit   *digits = var->digits;
5001         int                     ndigits = var->ndigits;
5002
5003         /* Strip leading zeroes */
5004         while (ndigits > 0 && *digits == 0)
5005         {
5006                 digits++;
5007                 var->weight--;
5008                 ndigits--;
5009         }
5010
5011         /* Strip trailing zeroes */
5012         while (ndigits > 0 && digits[ndigits - 1] == 0)
5013                 ndigits--;
5014
5015         /* If it's zero, normalize the sign and weight */
5016         if (ndigits == 0)
5017         {
5018                 var->sign = NUMERIC_POS;
5019                 var->weight = 0;
5020         }
5021
5022         var->digits = digits;
5023         var->ndigits = ndigits;
5024 }