www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.learn - Vectors and matrices in D

reply Artyom Shalkhakov <artyom.shalkhakov gmail.com> writes:
Hello!

First of all, thanks to WB for this awesome language! :)

I've got a couple of silly newbie questions.
I started implementing a little raytracer in D (raytracing is addictive, yeah
;)), and faced problems with vector/matrix math :(

struct Vec3 {
    float opIndex( int index ) {
        assert( index >= 0 && index < 3 );
        return ( &x )[index];
    }

    void opIndexAssign( float f, int index ) {
        assert( index >= 0 && index < 3 );
        ( &x )[index] = f;
    }

    Vec3 opAddAssign( inout Vec3 v ) {
        x += v.x;
        y += v.y;
        z += v.z;

        return *this;
    }
/*
    // on C++
    Vec3 &operator+=( const Vec3 &v ) {
        x += v.x;
        y += v.y;
        z += v.z;

        return *this;
    }
*/

    float x, y, z;
}

struct Mat3 {
    // I think that copying is costly here..
    Vec3 opIndex( int index ) {
        return rows[index];
    }
/*
    // this is how it looked like on C++
    Vec3 &operator[]( const int index ) {
        assert( index >= 0 && index < 3 );
        return rows[index];
    }
    const Vec3 &operator( const int index ) {
        assert( index >= 0 && index < 3 );
        return rows[index];
    }
*/

    void opIndexAssign( Vec3 v, int index ) {
        rows[index] = v;
    }

    private Vec3 rows[3];
}

unittest {
    Mat3   m;

    m[0][0] = 1.0f;
    assert( m[0][0] == 1.0f );
    // The above assertion fails all the time.
    // I think that DMD interprets m[0][0] as m.opIndex( 0 ).opIndexAssign( 0 ),
    // so it ends up operating on copy of row.
    // Any suggestions?
}

I tried other representations for vectors and matrices, but faced another
problem:
what happens when a static array is passed to function?

For example:
alias float[3]  Vec3;

float dot( Vec3 v1, Vec3 v2 ) {
    return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
}

If I get it correctly, the code above is identical to:

float dot( size_t length1, float *v1, size_t length2, float *v1 ) {
    return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
}

Lengths of two arrays will never ever change, therefore two parameters are
redundant. This is an overhead I'm trying to avoid. Can anyone tell me should I
really care?

Thanks in advance!
Mar 12 2007
next sibling parent mario pernici <mario.pernici mi.infn.it> writes:
Artyom Shalkhakov Wrote

 struct Vec3 {
     float opIndex( int index ) {
         assert( index >= 0 && index < 3 );
         return ( &x )[index];
     }
 
     void opIndexAssign( float f, int index ) {
         assert( index >= 0 && index < 3 );
         ( &x )[index] = f;
     }
 
     Vec3 opAddAssign( inout Vec3 v ) {
         x += v.x;
         y += v.y;
         z += v.z;
 
         return *this;
     }
 /*
     // on C++
     Vec3 &operator+=( const Vec3 &v ) {
         x += v.x;
         y += v.y;
         z += v.z;
 
         return *this;
     }
 */
 
     float x, y, z;
 }
 
 struct Mat3 {
     // I think that copying is costly here..
     Vec3 opIndex( int index ) {
         return rows[index];
     }
 /*
     // this is how it looked like on C++
     Vec3 &operator[]( const int index ) {
         assert( index >= 0 && index < 3 );
         return rows[index];
     }
     const Vec3 &operator( const int index ) {
         assert( index >= 0 && index < 3 );
         return rows[index];
     }
 */
 
     void opIndexAssign( Vec3 v, int index ) {
         rows[index] = v;
     }
 
     private Vec3 rows[3];
 }
 
 unittest {
     Mat3   m;
 
     m[0][0] = 1.0f;
     assert( m[0][0] == 1.0f );
     // The above assertion fails all the time.

Thsi works, though I am not sure it answers your question. struct Mat3 { float opIndex(size_t i, size_t j) { return rows[i][j]; } void opIndexAssign(float v, size_t i, size_t j ) { rows[i][j] = v; } private Vec3 rows[3]; } Mario
Mar 13 2007
prev sibling next sibling parent reply Bill Baxter <dnewsgroup billbaxter.com> writes:
Artyom Shalkhakov wrote:
 Hello!
 
 First of all, thanks to WB for this awesome language! :)
 
 I've got a couple of silly newbie questions.
 I started implementing a little raytracer in D (raytracing is addictive, yeah
;)), and faced problems with vector/matrix math :(
 
 struct Vec3 {
     float opIndex( int index ) {
         assert( index >= 0 && index < 3 );
         return ( &x )[index];
     }
 
     void opIndexAssign( float f, int index ) {
         assert( index >= 0 && index < 3 );
         ( &x )[index] = f;
     }
 
     Vec3 opAddAssign( inout Vec3 v ) {
         x += v.x;
         y += v.y;
         z += v.z;
 
         return *this;
     }
 /*
     // on C++
     Vec3 &operator+=( const Vec3 &v ) {
         x += v.x;
         y += v.y;
         z += v.z;
 
         return *this;
     }
 */
 
     float x, y, z;
 }
 
 struct Mat3 {
     // I think that copying is costly here..
     Vec3 opIndex( int index ) {
         return rows[index];
     }
 /*
     // this is how it looked like on C++
     Vec3 &operator[]( const int index ) {
         assert( index >= 0 && index < 3 );
         return rows[index];
     }
     const Vec3 &operator( const int index ) {
         assert( index >= 0 && index < 3 );
         return rows[index];
     }
 */
 
     void opIndexAssign( Vec3 v, int index ) {
         rows[index] = v;
     }
 
     private Vec3 rows[3];
 }
 
 unittest {
     Mat3   m;
 
     m[0][0] = 1.0f;
     assert( m[0][0] == 1.0f );
     // The above assertion fails all the time.
     // I think that DMD interprets m[0][0] as m.opIndex( 0 ).opIndexAssign( 0
),
     // so it ends up operating on copy of row.
     // Any suggestions?
 }
 
 I tried other representations for vectors and matrices, but faced another
problem:
 what happens when a static array is passed to function?
 
 For example:
 alias float[3]  Vec3;
 
 float dot( Vec3 v1, Vec3 v2 ) {
     return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
 }
 
 If I get it correctly, the code above is identical to:
 
 float dot( size_t length1, float *v1, size_t length2, float *v1 ) {
     return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
 }
 
 Lengths of two arrays will never ever change, therefore two parameters are
redundant. This is an overhead I'm trying to avoid. Can anyone tell me should I
really care?
 
 Thanks in advance!

Some general suggestions: * Search the NG archives for "raytracer". There was a thread about a simple raytracer implementation in C++ and D and its performance in the two. There was code for a working raytracer that might be useful for you as an example, but also in the thread numerous performance issues were discussed, like for instance the cost passing Vec3 type structs around by value. * Look at helix for one way to implement vector/matrix classes. http://www.dsource.org/projects/helix --bb
Mar 13 2007
parent KlausO <oberhofer users.sourceforge.net> writes:
Bill Baxter wrote:
 
 Some general suggestions:
 * Search the NG archives for "raytracer".  There was a thread about a 
 simple raytracer implementation in C++ and D and its performance in the 
 two.  There was code for a working raytracer that might be useful for 
 you as an example, but also in the thread numerous performance issues 
 were discussed, like for instance the cost passing Vec3 type structs 
 around by value.
 
 * Look at helix for one way to implement vector/matrix classes. 
 http://www.dsource.org/projects/helix
 
 --bb

There is also the Sol project on dsource which implements a raytracer. See http://www.dsource.org/projects/sol and http://www.mrshoe.org/projects/sol/ KlausO
Mar 13 2007
prev sibling parent reply Andreas Zwinkau <beza1e1 web.de> writes:
 I started implementing a little raytracer in D (raytracing is addictive, yeah
;)), and faced problems with vector/matrix math :(

Me too :) I wanted to learn D and to write a raytracer for some time. Currently i have some time to do so. I got a little sidetracked into D templates by the vector implementation. Every implementation (Artyom, helix, some C++ raytracers) i've seen so far uses a struct for the vectors. Why? Performance? Does Bill Baxter hint at this? http://www.digitalmars.com/pnews/read.php?server=news.digitalmars.com&group=digitalmars.D&artnum=46958 I wrote a generic vector class. Maybe someone would like to comment on that? I it's just addition and subtraction so far, but the pattern is clear. Here comes my code: import std.stdarg; import std.string; class Vector(T, int S) { T[S] content; this(...) { for (int i = 0; i < _arguments.length; i++) { content[i] = va_arg!(T)(_argptr); } } T opIndex(int i) { assert(i < content.length, std.string.format("%d < %d !!!\n", i, content.length)); return content[i]; } void opIndexAssign(T n, int i) { content[i] = n; } Vector!(T,S) opAdd( Vector!(T,S) other ) { Vector!(T,S) n = new Vector!(T,S)(); foreach(int i, T c; content) { n[i] = c + other[i]; } return n; } Vector!(T,S) opSub( Vector!(T,S) other ) { Vector!(T,S) n = new Vector!(T,S)(); foreach(int i, T c; content) { n[i] = c - other[i]; } return n; } } unittest { Vector!(int, 2) a = new Vector!(int, 2)(1,3); Vector!(int, 2) b = new Vector!(int, 2)(1,2); // check addition Vector!(int, 2) c = a + b; assert(c[0] == 2); assert(c[1] == 5); // check subtraction Vector!(int, 2) d = a - b; assert(d[0] == 0); assert(d[1] == 1); } alias Vector!(float,3) Vector3D; alias Vector!(float,2) Vector2D; unittest { Vector3D x = new Vector3D(2,3,4); Vector3D y = x + x; assert(y[2] - 8 < 0.001); }
Mar 15 2007
parent reply Bill Baxter <dnewsgroup billbaxter.com> writes:
Andreas Zwinkau wrote:
 I started implementing a little raytracer in D (raytracing is addictive, yeah
;)), and faced problems with vector/matrix math :(

Me too :) I wanted to learn D and to write a raytracer for some time. Currently i have some time to do so. I got a little sidetracked into D templates by the vector implementation. Every implementation (Artyom, helix, some C++ raytracers) i've seen so far uses a struct for the vectors. Why? Performance? Does Bill Baxter hint at this? http://www.digitalmars.com/pnews/read.php?server=news.digitalmars.com&group=digitalmars.D&artnum=46958 I wrote a generic vector class. Maybe someone would like to comment on that? I it's just addition and subtraction so far, but the pattern is clear. Here comes my code:

I think the reason people use structs for Vector types is because structs are for Plain Old Data, and Vectors are Plain Old Data. Mainly the differences that they're constructed on the stack by default which is generally what you want with something like a Vector (for 4-vectors and lower at least). Other things also make them more struct like -- you want A = B to make a copy of B, not make a second reference to the same B. You don't need to be able to subclass them. The main problem is that passing big structs by value as parameters is kind of slow. And it's something you want a lot. In C++ you'd use const& Vector for the argument types, but in D there's no equivalent. The good thing about using classes is that they're passed by reference by default. But that's about the only good thing about them in this context. Allocations on the heap can be slow. Using 'scope' can get you around that, but then you'll be using "scope Vector this", "scope Vector that" all the time and it will get tiresome. Anyway, if you're coming from Java then doing a heap allocation for every vector may not seem like a big deal. But Java's allocations of small objects are pretty heavily optimized from what I've heard. D isn't really there. Sorry if this is incoherent... I probably should have gone to be an hour ago. --bb
Mar 15 2007
next sibling parent Sean Kelly <sean f4.ca> writes:
Bill Baxter wrote:
 
 Anyway, if you're coming from Java then doing a heap allocation for 
 every vector may not seem like a big deal.  But Java's allocations of 
 small objects are pretty heavily optimized from what I've heard.  D 
 isn't really there.

I think some Java compilers will actually allocate classes on the stack if they can prove the object does not escape scope. This is similar to what we do manually in D with the 'scope' keyword. That said, there are still many instances where automated stack allocation is typically not possible in Java (returning data from a function, for example), but can still be done in D using structs. As an example of how much of an impact dynamic allocation has on application run times, I did a short performance analysis using Java a while back. I'll include the summary plus some timings below. ---------------------------------------------------------------------- From a performance perspective, the array-based and link-based queue implementations are effectively identical but for the enqueue method. A link-based enqueue operation allocates a new node object then performs some reference manipulation to append this node to an internal linked-list. By comparison, the array-based implementation will resize its internal array by doubling its length and copying existing elements into the new area if the array is full, and will then copy the new value into the array directly. Therefore, assuming all operations have an equal cost, the performance difference between these two implementations rests in how often the array-queue must grow to accommodate new elements. Assuming a worst-case scenario where elements are repeatedly inserted into the queue and are never removed, the array-based implementation must perform log(N) resize operations, and (N - 1) elements copied, assuming N represents the maximum number of elements inserted into the queue. This is a cost not incurred by the link-queue, where no copy operation is performed regardless of queue growth. Therefore, in terms of N, the array-queue will perform (N + log(N)) more operations than the link-queue in a worst-case scenario. So assuming the queue grows to 1000 elements, ten resize operations will be performed to produce a maximum array size of 1024 elements, and 1023 element copy operations will be performed, for a total of 1033 additional operations, or approximately twice as many operations as the linked-queue implementation. Therefore, it is reasonable to assume that the array-queue will perform approximately twice as slowly as the link-queue, assuming this worst-case scenario. Performance results are rarely as simple as such a code analysis suggests however. Memory allocations are typically fairly expensive operations, as are any garbage collection runs required to clean up discarded memory. Also, arrays typically have far better locality of reference than linked-lists (linked-list locality depends on allocator design and on how the container is used within the application), so fewer cache misses typically occur with array-based containers than with linked-list containers. Finally, the array-queue described above will actually resize very rarely compared to the number of insertion operations performed, giving what is sometimes referred to as "amortized constant time" performance for the enqueue operation. Given all these considerations, I think real-world use of these two queue implementations will show approximately equal performance, and the array-queue has the potential to be even faster than the link-queue in some situations, or on some systems. For three runs inserting one million elements in succession and then removing them (similar to the worst case for array-queue as described above), my average times were: array-queue: 153.3 milliseconds link-queue: 641 milliseconds
Mar 15 2007
prev sibling parent reply Jascha Wetzel <"[firstname]" mainia.de> writes:
Bill Baxter wrote:
 The main problem is that passing big structs by value as parameters is
 kind of slow.  And it's something you want a lot.  In C++ you'd use
 const& Vector for the argument types, but in D there's no equivalent.

what about inout? it does exactly that. DMD will even pass single inout arguments in a register.
Mar 15 2007
next sibling parent reply Jascha Wetzel <"[firstname]" mainia.de> writes:
ok, the const part is missing, but the performance part can be dealt
with. maybe "const inout" is a worthwhile feature for future D versions.
it should have a less contradictory name, though.

Jascha Wetzel wrote:
 Bill Baxter wrote:
 The main problem is that passing big structs by value as parameters is
 kind of slow.  And it's something you want a lot.  In C++ you'd use
 const& Vector for the argument types, but in D there's no equivalent.

what about inout? it does exactly that. DMD will even pass single inout arguments in a register.

Mar 15 2007
parent reply Frits van Bommel <fvbommel REMwOVExCAPSs.nl> writes:
Jascha Wetzel wrote:
 ok, the const part is missing, but the performance part can be dealt
 with. maybe "const inout" is a worthwhile feature for future D versions.
 it should have a less contradictory name, though.

Hey, no fair correcting yourself before I finish my post! ;) Seriously: as I mentioned in my post there was some discussion a while back about adding this sort of stuff, but so far nothing concrete has shown up.
Mar 15 2007
parent reply Jascha Wetzel <"[firstname]" mainia.de> writes:
;)

maybe we could have an additional "inref" modifier. on the other hand
"out" doesn't force write-only either, atm.
therefore i think having a way to return references/lvalues is more
important on that front.
i'd really like to do matrix[row,col] += bla;

Frits van Bommel wrote:
 Jascha Wetzel wrote:
 ok, the const part is missing, but the performance part can be dealt
 with. maybe "const inout" is a worthwhile feature for future D versions.
 it should have a less contradictory name, though.

Hey, no fair correcting yourself before I finish my post! ;) Seriously: as I mentioned in my post there was some discussion a while back about adding this sort of stuff, but so far nothing concrete has shown up.

Mar 15 2007
parent Frits van Bommel <fvbommel REMwOVExCAPSs.nl> writes:
[corrected upside-down post]

Jascha Wetzel wrote:
 Frits van Bommel wrote:
 Jascha Wetzel wrote:
 ok, the const part is missing, but the performance part can be dealt
 with. maybe "const inout" is a worthwhile feature for future D versions.
 it should have a less contradictory name, though.

Seriously: as I mentioned in my post there was some discussion a while back about adding this sort of stuff, but so far nothing concrete has shown up.

;) maybe we could have an additional "inref" modifier. on the other hand "out" doesn't force write-only either, atm.

'out' does force initialization though. Before anything can be read from or written to the parameter, it's automatically set to the initial value of its type. So you can read it, but it won't have the value it had in the context of the caller (unless, of course, that value was the initial value :) ).
 therefore i think having a way to return references/lvalues is more
 important on that front.

This was IIRC also a topic in the discussion I mentioned.
 i'd really like to do matrix[row,col] += bla;

So would I. Though I think that syntax is currently possible with proxy objects, that's not the most elegant solution.
Mar 15 2007
prev sibling parent Frits van Bommel <fvbommel REMwOVExCAPSs.nl> writes:
Jascha Wetzel wrote:
 Bill Baxter wrote:
 The main problem is that passing big structs by value as parameters is
 kind of slow.  And it's something you want a lot.  In C++ you'd use
 const& Vector for the argument types, but in D there's no equivalent.

what about inout? it does exactly that.

No, it doesn't do "exactly that". 'inout Vector' is more like 'Vector&' (i.e. a mutable reference instead of a non-mutable reference). From a performance viewpoint they are equivalent, but not from a semantics viewpoint. A while back Andrei and Walter were talking about adding, among other things, the concept of non-mutable references to the language. It's been a bit quiet on that front though.
 DMD will even pass single inout arguments in a register.

That it will, unless the function is a member function (or delegate). Though in that case the this/context pointer may technically be considered an extra parameter. IIRC GDC for amd64 will put several parameters in registers though, since the calling convention seems to be based on the C one which uses several registers before spilling over to the stack on that architecture.
Mar 15 2007