Re: Numeric x^y for negative x

On Thu, 22 Jul 2021 at 16:19, Dean Rasheed <dean(dot)a(dot)rasheed(at)gmail(dot)com> wrote:
>
> On Thu, 22 Jul 2021 at 06:13, Yugo NAGATA <nagata(at)sraoss(dot)co(dot)jp> wrote:
> >
> > Thank you for updating the patch. I am fine with the additional comments.
> > I don't think there is any other problem left, so I marked it Ready-for-Committers.
>
> Thanks for looking. Barring any further comments, I'll push this in a few days.
>

So I have been testing this a lot over the last couple of days, and I
have concluded that the patch works well as far as it goes, but I did
manage to construct another case where numeric_power() loses
precision. I think, though, that it would be better to tackle that as
a separate patch.

In writing the commit message for this patch, I realised that it was
possible to tidy up the local_rscale calculation part of it a little,
to make it more obvious what was going wrong, so attached is a
slightly tweaked version. I'll hold off on committing it for a few
more days in case anyone else wants to have a look. Tom?

The other issue I found is related to the first part of power_var(),
where it does a low-precision calculation to get an estimate for the
weight of the result. It occurred to me that for certain input bases,
that calculation could be made to be quite inaccurate, and therefore
lead to choosing the wrong rscale later. This is the test case I
constructed:

(1-1.5123456789012345678e-1000) ^ 1.15e1003

Here, the base is a sliver under 1, and so ln(base) is approximately
-1.5e-1000, and ln_dweight is -1000 (the decimal weight of ln(base)).
The problem is that the local_rscale used for the first low-precision
calculation is limited to NUMERIC_MAX_DISPLAY_SCALE (which is 1000),
so we only compute ln_base to a scale of 1000 at that stage, and the
result is rounded to exactly 2e-1000, which is off by a factor of
around 1.333333.

That makes it think the result weight will be -998, when actually it's
-755, so it then chooses a local_rscale for the full calculation
that's far too small, and the result is very inaccurate.

To fix this, I think it's necessary to remove the line that limits the
initial local_rscale. I tried that in a debugger, and managed to get a
result that agreed with the result from "bc -l" with a scale of 2000.

The reason I think it will be OK to remove that line is that it only
ever comes into play when ln_dweight is a large negative number (and
the smallest it can be is -16383). But that can only happen in
instances like this, where the base is very very close to 1. In such
cases, the ln(base) calculation is very fast, because it basically
only has to do a couple of Taylor series terms, and it's done. This
will still only be a low-precision estimate of the result (about 8
significant digits, shifted down a long way).

It might also be necessary to re-think the choice of local_rscale for
the mul_var() that follows. If the weight of exp is much larger than
the weight of ln_base (or vice versa), it doesn't really make sense to
compute the product to the same local_rscale. That might be a source
of other inaccuracies. I'll try to investigate some more.

Anyway, I don't think any of that should get in the way of committing
the current patch.

Regards,
Dean