How to handle reductions on interval arrays

[Joel Dahne]

Hi Oliver,

A new week and some new questions!

I have almost finished the work on vectorization. The only thing that is
left is how to handle ctc_union.m, ctc_intersection.m and
fsolve.m. These functions are a bit more complicated however so I will
wait a bit with them.

Instead I have started to work on folding/reduction functions (Octave
seems to use the word reduction in the source code at least). I have
implemented support for N-dimensional arrays in sum.m (and also in
mpfr_vector_sum_d.cc). In that change there are also some other changes,
I have moved all of the vectorization inside the oct-file (yielding
about a 2 times speed up) and have slightly modified how it handles empty
arrays to be consistent with how the standard sum-function does it.

When I were to work on dot.m I realized that Octaves standard
implementation is a bit odd. It is inconsistent with both the
sum-function and Matlab. See bug #51333,
https://savannah.gnu.org/bugs/?51333, on Savannah for more detailed
info. I think that the interval package should follow Octave, but in
this case Octave is weird. Probably we should wait and see what happens
with my bug report.

I have also worked on prod.m and while doing that I encountered
something which I think is a bug. But I'm not sure so I thought I would
ask first. A minimal example is the code

> x = infsup (1:4);
> i = logical ([0 1 0 0]);
> x(i)
ans = [2]

where everything is as it should. However when I try to put this code
inside an interval function, for example testme.m, it does not
work. With

function testme (arg)
  x = infsup (1:4);
  i = logical ([0 1 0 0]);
  x(i);
endfunction

we get

> testme (infsup ([]))
error: x(2): out of bound 1
error: called from
    testme at line 7 column 3

It is important that testme.m is inside the @infsup folder, otherwise it
seems to work. Do you have any idea of what this is?

Best,
Joel Dahne