Array operations -- not actually so difficult

[Sean Kelly]

Bill Baxter wrote:
> Don Clugston wrote:
>> Array operations are one of the very few 2.0 features that haven't 
>> made it into 1.0. My experiments into expression templates suggest 
>> that it might not be very complicated after all.
>>
>> Consider something like a dot product,
>>
>> real a[], b[], c[];
>> real d = dot(a+b, c);
>>
>> For performance, it would be better if the compiler actually *didn't* 
>> evaluate a+b. Instead, that should compile to:
>>
>> real sum = 0;
>> for (int i=0; i<c.length; ++i) {
>>      sum+=(a[i]+b[i])*c[i];
>> }
>> d=sum;
>>
>> In fact, we could create a fantastic library implementation of array 
>> operations with a tiny bit of syntactic sugar:
>> instead of generating an error message for array operations, look for 
>> an opXXX function, and treat it the same way the implicit array 
>> properties work.
> 
> How does that result in the above code for dot product?
> Looks to me like if all you have is calling opXXX functions the above 
> dot would still become two loops with a tmp for the a+b.

Not necessarily.  Consider:

     struct AddExp(T, U)
     {
         T lhs;
         U rhs;

         static AddExp!(T,U) opCall( T lhs, U rhs )
         in
         {
             assert( lhs.length == rhs.length );
             assert( is( typeof(T[0]) == typeof(U[0]) ) );
         }
         body
         {
             AddExp!(T,U) exp;
             exp.lhs = lhs;
             exp.rhs = rhs;
             return exp;
         }

         typeof(T[0]) opIndex( size_t i )
         {
             return lhs[i] + rhs[i];
         }

         typeof(T[0])[] opSlice( size_t b, size_t e )
         {
             typeof(T[0])[] sum = new typeof(T[0])[lhs.length];

             for( size_t i = 0; i < lhs.length; ++i )
             {
                 sum[i] = lhs[i] + rhs[i];
             }
             return sum;
         }

         size_t length()
         {
             return lhs.length;
         }
     }


     AddExp!(T,U) opAdd(T, U)( T lhs, U rhs )
     {
         return AddExp!(T,U)( lhs, rhs );
     }


     T[] opAssign(T, U)( inout T[] lhs, U rhs )
     in
     {
         assert( is( typeof(T[0]) == typeof(U[0]) ) );
     }
     body
     {
        lhs = rhs[0 .. $];
        return lhs;
     }


     real dotp(T,U)( T lhs, U rhs )
     {
         real sum = 0.0;
         assert( lhs.length == rhs.length );
         for( size_t i = 0; i < lhs.length; ++i )
         {
             sum += lhs[i] * rhs[i];
         }
         return sum;
     }


     void main()
     {
         real[] a, b, c;
         a ~= 1.0; a ~= 2.0;
         b ~= 3.0; b ~= 4.0;
         c ~= 5.0; c ~= 6.0;

         // test opSlice
         auto   s = opAdd( a, b );
         real[] r = s[0 .. s.length];

         for( size_t i = 0; i < r.length; ++i )
         {
             printf( "%Lf = %Lf\n", r[i], s[i] );
         }

         // compound expression
         real res = dotp( opAdd( a, b ), c );

         printf( "\n%Lf\n", res );
     }

prints:

     4.000000 = 4.000000
     6.000000 = 6.000000

     56.000000

Thanks to the wonder of templates, a properly designed expression struct 
could be evaluated as-is rather than needing to be converted to an array 
before being passed to dotp().  If we were able to define free-floating 
operators as per Don's suggestion, the above would be darn slick :-)

By the way, I originally tried slicing using the standard syntax "s[0 .. 
$]" and got this error:

     test.d(84): Error: undefined identifier __dollar

I imagine this is a known issue?


Sean