Re: integer arithmetics

[Jaroslav Hajek]

On Wed, Oct 1, 2008 at 10:36 PM, John W. Eaton <address@hidden> wrote:
> On  1-Oct-2008, Jaroslav Hajek wrote:
>
> | On Wed, Oct 1, 2008 at 5:14 PM, John W. Eaton <address@hidden> wrote:
> | > On 23-Sep-2008, Jaroslav Hajek wrote:
> | >
> | > | Sorry. I didn't realize that std::abs (int) is declared in cstdlib,
> | > | not cmath. Please try again, same location. This new one should also
> | > | get rid of the unused vars warnings you see.
> | >
> | > OK.
> | >
> | > I applied this changeset.
> | >
> |
> | Hmm, initially I didn't mean this to go into 3.2.x, but if you're fine
> | with that, then why not.
>
> Does it work?  If so, then I don't see why not.  I also built it with
> HAVE_FAST_INT_OPS undefined and it seemed to build and function OK, at
> least for the minimal tests I tried.
>
> | Please consider also the attached patch, that
> | fully implements the 64-bit arithmetics.
> | If at least 80-bit wide long double type is discovered, it is used to
> | implement 64-bit multiplication and mixed operations, otherwise, more
> | complicated non-inline code is employed. In the latter case, some
> | operations may suffer from precision loss (converting int64->double)
> | or early truncation - see comments in oct-inttypes.cc. I think that
> | the important cases (int64 +- double) are OK, and that one can live
> | with losing precision when doing, for instance, 4611686018427387901 *
> | 0.75. Maybe these emulations will be improved over time (though at a
> | performance cost).
>
> Do you mean this comment?
>
> | +  // XXX: This trick is to not lose precision for things like 3*2**63 + 
> intmin('int64')
>
> Can you elaborate?

OK, here's a detailed problem description: The problem is that
convertnig 64-bit integers to doubles is a lossy operation. If 64-bit
ints were handled exactly like 32-bit and smaller, then, for instance,
4611686018427387901 (int64) + (double) 1 would produce
4611686018427387904 = 2**62
Though this may be easily considered "a feature", it would cause
expressions like
a = a + 2 to produce unexpected results.

If we have at least 80 bits wide long double (gcc on x86 archs), the
problem is easily solved, as all 64-bit integers fit into long double
without losing precision.
If not, we're in trouble.
Addition and subtraction are still relatively easy to handle. When
doing x + y (x is octave_int64, y is double), I first check whether y
is within the range of int64. If yes, then x + octave_int64(y)
produces exactly the same results as the long double version would.
If not, the situation is more complicated, because, for instance,
expression like 3*2**62 + (1+intmin('int64')) should produce 2**62 +
1.
In other words, if the double operand is within twice the signed
integer range, the result may still not overflow. Here I use the fact
that if y is outside the signed int64 range, it must be divisible by 2
without going fractional. So I split it into two halfs (each of which
must be within int64 range or the result would overflow anyway) and do
two converts and two additions. The situation is slightly different
for unsigned integers (no /2 trick) but the code should also work
always.

I have not yet, however, figured out a similar bullet-proof trick for
int64 * double.

All this mess is caused by the wicked way how Matlab treats mixed
integer-real arithmetics. I guess they only realised too late that it
won't work for 64-bit integers. Perhaps that's why int64 is so
horribly crippled in Matlab 2007 - it's virtually useless.

>
> Also, please use FIXME instead of XXX since some people have
> complained in the past about spam filters flagging mails that contain
> XXX.

OK, no problem.

>
> | OK to push?
>
> Probably not yet.  What happened here:
>
> | diff -r 1651939a394b -r 836df84b7997 doc/interpreter/contributors.in
> | --- a/doc/interpreter/contributors.in Wed Oct 01 11:34:45 2008 -0400
> | +++ b/doc/interpreter/contributors.in Wed Oct 01 20:16:39 2008 +0200
> | @@ -193,7 +193,6 @@
> |  Doug Stewart
> |  Jonathan Stickel
> |  Thomas Stuart
> | -Ivan Sutoris
> |  John Swensen
> |  Ariel Tankus
> |  Georg Thimm
>
> ?  I just added Ivan to this list.
>
> jwe
>

Sorry, mea culpa. The problem was that I didn't quite expect you to
apply the first patch, and I have already finished the 64-bit integers
and refreshed it. So what I did was to pull from main repo, then a "hg
diff tip" and tried to strip the unrelated diffs (backward changes),
but I overlooked this one.

I'll resend an update.

regards

-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz