Re: numeric precision when raising one numeric to another.

It appears from checking the output of exponentiation of one numeric to
another, the output is actually in floating point. Is this normal and /
or expected?

Now, given that

create table test2 (i1 numeric(20,0), i2 numeric(20,0));
insert into test values (123456789012345,123456789012345);
select i1*i2 from test2;
gives:
?column?
-------------------------------
15241578753238669120562399025

it seems odd that

create table test (i1 numeric(20,0), i2 numeric(20,0));
insert into test values (2,55);
select i1^i2 from test;
gives:
?column?
---------------------
3.6028797018964e+16

Now, I can get an exact answer if I'm willing to twiddle with breaking
the exponent down:

select (2^60)::numeric;

Gives:
numeric
---------------------
1152921504606850000

While, select (2^30)::numeric*(2^30)::numeric;

Gives:
?column?
---------------------
1152921504606846976

So, numeric can hold the value, but it looks like the exponent math is
converting it to float.

I'm not bothered too much by it, as I don't really work with numbers
that big. I was mainly wondering if this is kosher is all.

Scott Marlowe <smarlowe(at)g2switchworks(dot)com> writes:
> It appears from checking the output of exponentiation of one numeric to
> another, the output is actually in floating point. Is this normal and /
> or expected?

Yes, seeing that the only ^ operator we have is float8.

regression=# \do ^
List of operators
Schema | Name | Left arg type | Right arg type | Result type |
Description
------------+------+------------------+------------------+------------------+----------------------
pg_catalog | ^ | double precision | double precision | double precision | exponentiation (x^y)
(1 row)

regards, tom lane

On Wed, 2005-05-18 at 16:42, Tom Lane wrote:
> Scott Marlowe <smarlowe(at)g2switchworks(dot)com> writes:
> > It appears from checking the output of exponentiation of one numeric to
> > another, the output is actually in floating point. Is this normal and /
> > or expected?
>
> Yes, seeing that the only ^ operator we have is float8.
>
> regression=# \do ^
> List of operators
> Schema | Name | Left arg type | Right arg type | Result type |
> Description
> ------------+------+------------------+------------------+------------------+----------------------
> pg_catalog | ^ | double precision | double precision | double precision | exponentiation (x^y)
> (1 row)

But is this proper behaviour?

Considering that the SQL spec says the result of multiplication of exact
numeric types is exact numeric types of precision S1+S2, and
exponentiation is nothing more than repeated multiplication, should
postgresql have a numeric capable exponentiation operator? Since I've
finally got a job where I can actually hack on the clock a bit, it might
be a nice trial balloon. It'll take a week or two to knock the rust off
my C skills though. :)

Scott Marlowe <smarlowe(at)g2switchworks(dot)com> writes:
> Considering that the SQL spec says the result of multiplication of exact
> numeric types is exact numeric types of precision S1+S2, and
> exponentiation is nothing more than repeated multiplication,

... not when the exponent is non-integral.

regards, tom lane

>> Considering that the SQL spec says the result of multiplication of
>> exact
>> numeric types is exact numeric types of precision S1+S2, and
>> exponentiation is nothing more than repeated multiplication,
>
> ... not when the exponent is non-integral.

For one thing. For another, I believe the standard C library only has
floating point exponentiation functions, not that there aren't plenty
of numeric libraries with integral ones. Finally, exponentiated
numbers get real big, real fast, and the floating point types can hold
much larger magnitudes than the integer types, albeit inexactly. For
example, on the Mac I'm using now, long long ints max out at about
10^19, while long doubles can represent 10^308.

- John Burger
MITRE

On Wed, May 18, 2005 at 10:46:50PM -0400, John Burger wrote:
> >>Considering that the SQL spec says the result of multiplication of
> >>exact
> >>numeric types is exact numeric types of precision S1+S2, and
> >>exponentiation is nothing more than repeated multiplication,
> >
> >... not when the exponent is non-integral.
>
> For one thing. For another, I believe the standard C library only has
> floating point exponentiation functions, not that there aren't plenty
> of numeric libraries with integral ones. Finally, exponentiated
> numbers get real big, real fast, and the floating point types can hold
> much larger magnitudes than the integer types, albeit inexactly. For
> example, on the Mac I'm using now, long long ints max out at about
> 10^19, while long doubles can represent 10^308.

Well, we already have an interesting library of mathematical functions
for NUMERIC (which is an arbitrary precision type, so it wouldn't matter
how big the result would get). I think the only reason we don't have a
NUMERIC exponentiation function is that nobody has implemented it.

--
Alvaro Herrera (<alvherre[a]surnet.cl>)
"People get annoyed when you try to debug them." (Larry Wall)

On Wed, May 18, 2005 at 11:32:40PM -0400, Alvaro Herrera wrote:
> On Wed, May 18, 2005 at 10:46:50PM -0400, John Burger wrote:
> > For one thing. For another, I believe the standard C library only has
> > floating point exponentiation functions, not that there aren't plenty
> > of numeric libraries with integral ones. Finally, exponentiated
> > numbers get real big, real fast, and the floating point types can hold
> > much larger magnitudes than the integer types, albeit inexactly. For
> > example, on the Mac I'm using now, long long ints max out at about
> > 10^19, while long doubles can represent 10^308.
>
> Well, we already have an interesting library of mathematical functions
> for NUMERIC (which is an arbitrary precision type, so it wouldn't matter
> how big the result would get). I think the only reason we don't have a
> NUMERIC exponentiation function is that nobody has implemented it.

The prerequisites for such a function would be a log() and exp()
function for numeric. And the real question there would be, what's a
sufficient accuracy? Numbers people actually use rarely have even
rational logarithms, so there is no way to store them 100% accurate.

As long as you're using integral exponents you can get away with
multiplication. BTW, the commandline utility "bc" has arbitrary number
arithmatic, maybe we can see how they do it? It defaults to 20 digits
precision, which is obviously not enough for large exponents.

Hmm, it looks like even they don't support raising to fractional
powers. When calculating 2^100, you need a precision of at least 35
decimal places to get in the ballpark of the correct figure using
log/exp, 30 isn't enough. Maybe do exact for integer exponents and
approx for non-integer?

kleptog(at)vali:~$ bc -l
bc 1.06
Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
> 2^100
1267650600228229401496703205376
> 2^100.1
Runtime warning (func=(main), adr=11): non-zero scale in exponent
1267650600228229401496703205376
> e(l(2)*100)
1267650600228229400579922894637.90158245154400629512
> scale=30
> e(l(2)*100)
1267650600228229401496703205353.617337311111135194699059124092
> scale=35
> e(l(2)*100)
1267650600228229401496703205375.99897630874075350752485091801369515

Hope this helps,
--
Martijn van Oosterhout <kleptog(at)svana(dot)org> http://svana.org/kleptog/
> Patent. n. Genius is 5% inspiration and 95% perspiration. A patent is a
> tool for doing 5% of the work and then sitting around waiting for someone
> else to do the other 95% so you can sue them.

Why are we allowing implicit casts from numeric to floating point?
Doesn't that violate the principle of not doing any implicit casts that
would potentially drop precision? It seems that about half the arguments
here are related to getting unexpected or inaccurate results, presumably
from the implicit cast; if it was removed at least then people would
know their results might be imprecise, where as now there's no
indication of that at all.

On Wed, May 18, 2005 at 05:42:28PM -0400, Tom Lane wrote:
> Scott Marlowe <smarlowe(at)g2switchworks(dot)com> writes:
> > It appears from checking the output of exponentiation of one numeric to
> > another, the output is actually in floating point. Is this normal and /
> > or expected?
>
> Yes, seeing that the only ^ operator we have is float8.
>
> regression=# \do ^
> List of operators
> Schema | Name | Left arg type | Right arg type | Result type |
> Description
> ------------+------+------------------+------------------+------------------+----------------------
> pg_catalog | ^ | double precision | double precision | double precision | exponentiation (x^y)
> (1 row)
>
> regards, tom lane
>
> ---------------------------(end of broadcast)---------------------------
> TIP 3: if posting/reading through Usenet, please send an appropriate
> subscribe-nomail command to majordomo(at)postgresql(dot)org so that your
> message can get through to the mailing list cleanly
>

--
Jim C. Nasby, Database Consultant decibel(at)decibel(dot)org
Give your computer some brain candy! www.distributed.net Team #1828

Windows: "Where do you want to go today?"
Linux: "Where do you want to go tomorrow?"
FreeBSD: "Are you guys coming, or what?"

"Jim C. Nasby" <decibel(at)decibel(dot)org> writes:
> Why are we allowing implicit casts from numeric to floating point?

Because the SQL spec requires it.

2) If the data type of either operand of a dyadic arithmetic op-
erator is approximate numeric, then the data type of the re-
sult is approximate numeric.

It doesn't say to throw an error for mixed-type arithmetic.

Now it also says

1) If the data type of both operands of a dyadic arithmetic opera-
tor is exact numeric, then the data type of the result is exact
numeric, ...

which you could take as requiring us to provide numeric equivalents of
every floating-point operator, but I don't find that argument very
convincing for operations that are inherently not going to give exact
results. The spec demands exact results from addition, subtraction,
and multiplication, but as soon as you get to division they punt; let
alone transcendental functions.

But having said that, I don't have a problem with putting in a
pg_operator entry for numeric_power. And if someone wants to improve
the scale factor calculations therein, go for it. But so far there's
been an extremely low signal-to-noise ratio in this thread ...

regards, tom lane

On Fri, May 20, 2005 at 01:03:08PM -0400, Tom Lane wrote:
> "Jim C. Nasby" <decibel(at)decibel(dot)org> writes:
> > Why are we allowing implicit casts from numeric to floating point?
>
> Because the SQL spec requires it.
>
> 2) If the data type of either operand of a dyadic arithmetic op-
> erator is approximate numeric, then the data type of the re-
> sult is approximate numeric.
>
> It doesn't say to throw an error for mixed-type arithmetic.
>
> Now it also says
>
> 1) If the data type of both operands of a dyadic arithmetic opera-
> tor is exact numeric, then the data type of the result is exact
> numeric, ...

But isn't NUMERIC exact numeric and not approximate?

> which you could take as requiring us to provide numeric equivalents of
> every floating-point operator, but I don't find that argument very
> convincing for operations that are inherently not going to give exact
> results. The spec demands exact results from addition, subtraction,
> and multiplication, but as soon as you get to division they punt; let
> alone transcendental functions.

ISTM what's more important than be exact is respecting precision. If I'm
remembering this correctly from high school, multiplying two numbers
each having 10 significant digits means you then have 20 significant
digits, so we should at least respect that. Which means
numeric(500)^numeric(500) should give an exact numeric(1000), which I
don't think is a given when casting to a double.

I'm not sure how this changes if you're using a fractional exponent. But
it seems like a pretty serious issue if you're doing financial
calculations and those are sometimes done in floating point under the
covers.
--
Jim C. Nasby, Database Consultant decibel(at)decibel(dot)org
Give your computer some brain candy! www.distributed.net Team #1828

Windows: "Where do you want to go today?"
Linux: "Where do you want to go tomorrow?"
FreeBSD: "Are you guys coming, or what?"

On Fri, 2005-05-20 at 12:03, Tom Lane wrote:
> "Jim C. Nasby" <decibel(at)decibel(dot)org> writes:
> > Why are we allowing implicit casts from numeric to floating point?
>
> Because the SQL spec requires it.
>
> 2) If the data type of either operand of a dyadic arithmetic op-
> erator is approximate numeric, then the data type of the re-
> sult is approximate numeric.
>
> It doesn't say to throw an error for mixed-type arithmetic.
>
> Now it also says
>
> 1) If the data type of both operands of a dyadic arithmetic opera-
> tor is exact numeric, then the data type of the result is exact
> numeric, ...
>
> which you could take as requiring us to provide numeric equivalents of
> every floating-point operator, but I don't find that argument very
> convincing for operations that are inherently not going to give exact
> results.

Are you saying that the exponent operator will return inexact results?
OR talking about other operators

> The spec demands exact results from addition, subtraction,
> and multiplication, but as soon as you get to division they punt; let
> alone transcendental functions.

If you're quoting the 92 spec, it seems to say that multiplication
precision is also implementation specific.

> But having said that, I don't have a problem with putting in a
> pg_operator entry for numeric_power. And if someone wants to improve
> the scale factor calculations therein, go for it.

OK, I'm gonna look at it this weekend. I might have some questions
before I really get anything working, this being my first real adventure
hacking pgsql.

> But so far there's
> been an extremely low signal-to-noise ratio in this thread ...

Really, I've found it quite informative. I see no reason to insult the
people who've contributed to it.

Scott Marlowe <smarlowe(at)g2switchworks(dot)com> writes:
> Are you saying that the exponent operator will return inexact results?

For a fractional exponent, it generally has to, because there is no
finite exact result.

> If you're quoting the 92 spec, it seems to say that multiplication
> precision is also implementation specific.

You're misreading it: the scale part is what's important. Precision
means the implementation gets to set a limit on the total number of
digits it will store. The scale rules effectively say that for add/sub/mul
you must either deliver an exact result or report overflow. But that
requirement is not placed on division, and (by implication) not on other
operations that are incapable of delivering exact results every time.

numeric_power can in theory deliver an exact answer when the exponent is
a positive integer. Division can deliver an exact answer in some cases
too --- but the spec doesn't say it must do so when possible. So I
would say that there is no spec requirement for special behavior for
integral exponents.

We could try to deliver an exact answer for an integral exponent by
selecting output scale = input scale times exponent. But that doesn't
work for any but very small exponents --- as the exponent gets bigger
you really have to drop fractional precision, or you're going to hit
overflow, which is not an improvement. (We do set a limit on total
number of digits...) So it's a question of tradeoffs, not black and
white.

regards, tom lane

On Fri, May 20, 2005 at 12:22:33PM -0500, Jim C. Nasby wrote:
> > which you could take as requiring us to provide numeric equivalents of
> > every floating-point operator, but I don't find that argument very
> > convincing for operations that are inherently not going to give exact
> > results. The spec demands exact results from addition, subtraction,
> > and multiplication, but as soon as you get to division they punt; let
> > alone transcendental functions.
>
> ISTM what's more important than be exact is respecting precision. If I'm
> remembering this correctly from high school, multiplying two numbers
> each having 10 significant digits means you then have 20 significant
> digits, so we should at least respect that. Which means
> numeric(500)^numeric(500) should give an exact numeric(1000), which I
> don't think is a given when casting to a double.

Wrong.

numeric(500) * numeric(500) = numeric(1000)
numeric(500) ^ numeric(500) = numeric(10 ^ 503) >> googleplex

You do not have enough memory to store the exact result. There are not
enough atoms in the universe to store this result. That's one reason
why you can't guarentee an exact result. Even numeric(20) ^ numeric(20)
= numeric( 10 ^ 22 )

> I'm not sure how this changes if you're using a fractional exponent. But
> it seems like a pretty serious issue if you're doing financial
> calculations and those are sometimes done in floating point under the
> covers.

Financial calculations are a red herring. They don't deal with less
than hundredths of a cent or more than trillions of dollars so 20
significant digits is easily enough. I would say to place an upper
limit at say 100 digits. It you want better, go get a real math
package.

Have a nice day,
--
Martijn van Oosterhout <kleptog(at)svana(dot)org> http://svana.org/kleptog/
> Patent. n. Genius is 5% inspiration and 95% perspiration. A patent is a
> tool for doing 5% of the work and then sitting around waiting for someone
> else to do the other 95% so you can sue them.

On Fri, May 20, 2005 at 01:03:08PM -0400, Tom Lane wrote:

> But having said that, I don't have a problem with putting in a
> pg_operator entry for numeric_power. And if someone wants to improve
> the scale factor calculations therein, go for it.

Oh, and while at it, it would be nice to solve the modulo bug that still
lurks there:

alvherre=# select 12345678901234567890 % 123;
?column?
----------
-45
(1 fila)

alvherre=# select 12345678901234567890 % 123::numeric(4,1);
?column?
----------
78.0
(1 fila)

alvherre=# select 12345678901234567890 % 123::numeric(3,0);
?column?
----------
-45
(1 fila)

alvherre=# select version();
version
----------------------------------------------------------------------------------------------
PostgreSQL 8.1devel on i686-pc-linux-gnu, compiled by GCC gcc (GCC) 3.3.6 (Debian 1:3.3.6-4)
(1 fila)

--
Alvaro Herrera (<alvherre[a]surnet.cl>)
"I call it GNU/Linux. Except the GNU/ is silent." (Ben Reiter)

Alvaro Herrera wrote:
> On Fri, May 20, 2005 at 01:03:08PM -0400, Tom Lane wrote:
>
> Oh, and while at it, it would be nice to solve the modulo bug that still
> lurks there:
>
> alvherre=# select 12345678901234567890 % 123;
> ?column?
> ----------
> -45
> (1 fila)
>
> alvherre=# select 12345678901234567890 % 123::numeric(4,1);
> ?column?
> ----------
> 78.0
> (1 fila)
>
> alvherre=# select 12345678901234567890 % 123::numeric(3,0);
> ?column?
> ----------
> -45
> (1 fila)
>
> alvherre=# select version();
> version
> ----------------------------------------------------------------------------------------------
> PostgreSQL 8.1devel on i686-pc-linux-gnu, compiled by GCC gcc (GCC) 3.3.6 (Debian 1:3.3.6-4)
> (1 fila)

I poked around on this one and found this in the comments in
numeric::mod_var():

/* ---------
* We do this using the equation
* mod(x,y) = x - trunc(x/y)*y
* We set rscale the same way numeric_div and numeric_mul do
* to get the right answer from the equation. The final result,
* however, need not be displayed to more precision than the inputs.
* ----------
*/

so I tried it:

test=> select 12345678901234567890 % 123;
?column?
----------
-45
(1 row)

test=> select 12345678901234567890 / 123;
?column?
--------------------
100371373180768845
(1 row)

test=> select 100371373180768845::numeric * 123::numeric;
?column?
----------------------
12345678901234567935
(1 row)

test=> select 12345678901234567890 - 12345678901234567935;
?column?
----------
-45
(1 row)

and I was quite surprised at the result. Basically, it looks like the
division is rounding _up_ the next integer on the /123 division, and
that is causing the modulo error. In fact, should the /123 round up
with numeric? I think there is an assumption in our code that div_var()
will not round up, but in fact it does in this case.

Here is 'calc' showing the same calculation:

> 12345678901234567890 % 123
78
> 12345678901234567890 / 123
~100371373180768844.63414634146341463414
> 100371373180768845 * 123

^^^^^^^^^^^^^^^^^^ rounded up by me

12345678901234567935
> 12345678901234567890 - 12345678901234567935
-45

and here is 'bc' doing integer division:

12345678901234567890 / 123
100371373180768844
100371373180768844 * 123
12345678901234567812
12345678901234567890 - 12345678901234567812
78

This is why 123::numeric(4,1) fixes it because the division returns on
digit that is truncated, rather than rounding up to the next whole
number.

I am not sure how to fix this. Adding extra scale to the division would
help, but if the division returned .999 and we had a scale of 2, it
would still round up and the truncate would not see it.

--
Bruce Momjian | http://candle.pha.pa.us
pgman(at)candle(dot)pha(dot)pa(dot)us | (610) 359-1001
+ If your life is a hard drive, | 13 Roberts Road
+ Christ can be your backup. | Newtown Square, Pennsylvania 19073

On Mon, May 30, 2005 at 11:29:48PM -0400, Bruce Momjian wrote:

> test=> select 12345678901234567890 / 123;
> ?column?
> --------------------
> 100371373180768845
> (1 row)

Well, that's a bug, right?

--
Alvaro Herrera (<alvherre[a]surnet.cl>)
"Industry suffers from the managerial dogma that for the sake of stability
and continuity, the company should be independent of the competence of
individual employees." (E. Dijkstra)

Alvaro Herrera wrote:
> On Mon, May 30, 2005 at 11:29:48PM -0400, Bruce Momjian wrote:
>
> > test=> select 12345678901234567890 / 123;
> > ?column?
> > --------------------
> > 100371373180768845
> > (1 row)
>
> Well, that's a bug, right?

I don't think so. The fuller answer is
100371373180768844.63414634146341463414, and that rounded to the nearest
integer is 100371373180768845. I think people expect % do to that,
except for integers. You could argue that numerics with zero scale are
integers, but NUMERIC % NUMERIC doesn't behave like an integer operator
--- it rounds to the proper precision.

--
Bruce Momjian | http://candle.pha.pa.us
pgman(at)candle(dot)pha(dot)pa(dot)us | (610) 359-1001
+ If your life is a hard drive, | 13 Roberts Road
+ Christ can be your backup. | Newtown Square, Pennsylvania 19073

Bruce Momjian wrote:
> Alvaro Herrera wrote:
> > On Mon, May 30, 2005 at 11:29:48PM -0400, Bruce Momjian wrote:
> >
> > > test=> select 12345678901234567890 / 123;
> > > ?column?
> > > --------------------
> > > 100371373180768845
> > > (1 row)
> >
> > Well, that's a bug, right?
>
> I don't think so. The fuller answer is
> 100371373180768844.63414634146341463414, and that rounded to the nearest
> integer is 100371373180768845. I think people expect % do to that,
> except for integers. You could argue that numerics with zero scale are
> integers, but NUMERIC % NUMERIC doesn't behave like an integer operator
> --- it rounds to the proper precision.

Attached is a fix for the problem:

test=> select 12345678901234567890 % 123;
?column?
----------
78
(1 row)

I could have just forced X digits after the decimal point, but there was
really no _right_ number of digits, so it would have been a hack. The
proper solution is to change div_var to accept 'round' as true/false,
and that's what I did in the patch.

--
Bruce Momjian | http://candle.pha.pa.us
pgman(at)candle(dot)pha(dot)pa(dot)us | (610) 359-1001
+ If your life is a hard drive, | 13 Roberts Road
+ Christ can be your backup. | Newtown Square, Pennsylvania 19073

On 5/20/2005 2:26 PM, Tom Lane wrote:

> numeric_power can in theory deliver an exact answer when the exponent is
> a positive integer. Division can deliver an exact answer in some cases
> too --- but the spec doesn't say it must do so when possible. So I
> would say that there is no spec requirement for special behavior for
> integral exponents.

There are cases where a numeric_power could in theory deliver an exact
answer for a fractional exponent. That is when the exponent is a natural
fraction because the result is the m'th root of x^n (for n/m). As an
example 4^1.5 = 8. Of course does the m'th root need to produce a finite
result, which I think is not guaranteed for arbitrary numbers.

I'm not advocating to do that, just saying it is theoretically possible
for a subset of possible inputs.

Jan

--
#======================================================================#
# It's easier to get forgiveness for being wrong than for being right. #
# Let's break this rule - forgive me. #
#================================================== JanWieck(at)Yahoo(dot)com #