Re: Does Ada need elemental functions to make it suitable for scientific work?

[Brian Drummond <brian@shapes.demon.co.uk>]

On Tue, 10 Jul 2012 01:02:13 -0700, gautier.de.montmollin wrote:

> Le mardi 10 juillet 2012 07:22:58 UTC+2, Ada novice a écrit :
> 
>> Believe me, for someone who do numerical computations on a daily basis,
>> an operation like V**2 is a NECESSITY. It's a pity that Ada does
>> not offer such a facility. I do not blame those who stick to Matlab (or
>> Fortran?) for doing numerical computations.
> 
> Wait... There are *such* facilities and several built-in operations
> *like* that:
> http://www.ada-auth.org/standards/12rm/html/RM-G-3-1.html
> 
> Apparently it is not enough. But it is possible to make a proposal for
> the next standard (and short term for the next GNAT version).

Furthermore, they are not really "built-in", so there is no need to go 
back to the standards authority if you feel they are lacking. The source 
for Ada.Numerics.Generic_Real_Array is right there, and contains no magic.

Adding "**" to suit your own purposes should be no more complicated than
the implementation of "+" in Ada.Numerics.Generic_Real_Array.

So, how does that work?

Ada.Numerics.Generic_Real_Array is found (on my machine) at
/usr/lib/gcc/i586-suse-linux/4.6/adainclude/a-ngrear.ads (specification) 
and /usr/lib/gcc/i586-suse-linux/4.6/adainclude/a-ngrear.adb (body)

The specification is
   function "+"   (Right : Real_Vector)       return Real_Vector;
...
   pragma Inline_Always ("+");

The body - you might expect a loop , performing addition on each element.
Well, not quite...

----------------------------------------------------------------------
with System.Generic_Array_Operations; use System.Generic_Array_Operations;

...
   --  Instantiating the following subprograms directly would lead to
   --  name clashes, so use a local package.

   package Instantiations is
...
      function "+" is new
        Vector_Vector_Elementwise_Operation
          (Left_Scalar   => Real'Base,
           Right_Scalar  => Real'Base,
           Result_Scalar => Real'Base,
           Left_Vector   => Real_Vector,
           Right_Vector  => Real_Vector,
           Result_Vector => Real_Vector,
           Operation     => "+");
...
   end Instantiations;

   function "+" (Left, Right : Real_Vector) return Real_Vector
      renames Instantiations."+";

----------------------------------------------------------------------

So it uses a generic function from System.Generic_Array_Operations to 
loop over the vector, applying the Operation argument to each element.

And this generic function is available for you to use in the same manner.
(There are similar ones for unary operations like sin(), and others for 
2D matrix operations, etc)

So you can create a package along the same lines, and add **, sin, etc to 
it. Once you have the package in place, adding more matrix operations of 
your own is trivial.

What does this generic function look like?
Excerpt from System.Generic_Array_Operations - on this system, at
/usr/lib/gcc/i586-suse-linux/4.6/adainclude/s-gearop.ads(specification)
--------------------------------
   generic
      type Left_Scalar is private;
      type Right_Scalar is private;
      type Result_Scalar is private;
      type Left_Vector is array (Integer range <>) of Left_Scalar;
      type Right_Vector is array (Integer range <>) of Right_Scalar;
      type Result_Vector is array (Integer range <>) of Result_Scalar;
      with function Operation
             (Left  : Left_Scalar;
              Right : Right_Scalar) return Result_Scalar;
   function Vector_Vector_Elementwise_Operation
     (Left  : Left_Vector;
      Right : Right_Vector) return Result_Vector;
--------------------------------
which should be enough information to let you use it...

If you are curious about its implementation, look at
/usr/lib/gcc/i586-suse-linux/4.6/adainclude/s-gearop.adb (body)
--------------------------------
   -----------------------------------------
   -- Vector_Vector_Elementwise_Operation --
   -----------------------------------------

   function Vector_Vector_Elementwise_Operation
     (Left  : Left_Vector;
      Right : Right_Vector) return Result_Vector
   is
      R : Result_Vector (Left'Range);

   begin
      if Left'Length /= Right'Length then
         raise Constraint_Error with
            "vectors are of different length in elementwise operation";
      end if;

      for J in R'Range loop
         R (J) := Operation (Left (J), Right (J - R'First + Right'First));
      end loop;

      return R;
   end Vector_Vector_Elementwise_Operation;
--------------------------------
which is of course the expected loop, plus a sanity check.

- Brian