## digitalmars.D.learn - Vectors and matrices in D

• Artyom Shalkhakov (73/73) Mar 12 2007 Hello!
• mario pernici (12/74) Mar 13 2007 Thsi works, though I am not sure it answers your question.
• Bill Baxter (11/104) Mar 13 2007 Some general suggestions:
• KlausO (7/20) Mar 13 2007 There is also the Sol project on dsource which implements a raytracer.
• Andreas Zwinkau (56/57) Mar 15 2007 Me too :)
• Bill Baxter (24/34) Mar 15 2007 I think the reason people use structs for Vector types is because
• Sean Kelly (56/61) Mar 15 2007 I think some Java compilers will actually allocate classes on the stack
• Jascha Wetzel (3/6) Mar 15 2007 what about inout? it does exactly that.
• Jascha Wetzel (4/11) Mar 15 2007 ok, the const part is missing, but the performance part can be dealt
• Frits van Bommel (5/8) Mar 15 2007 Hey, no fair correcting yourself before I finish my post! ;)
• Frits van Bommel (14/21) Mar 15 2007 No, it doesn't do "exactly that". 'inout Vector' is more like 'Vector&'
```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;
}

unittest {
Mat3   m;

m = 1.0f;
assert( m == 1.0f );
// The above assertion fails all the time.
// I think that DMD interprets m 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  Vec3;

float dot( Vec3 v1, Vec3 v2 ) {
return v1 * v2 + v1 * v2 + v1 * v2;
}

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 * v2 + v1 * v2 + v1 * v2;
}

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?

```
Mar 12 2007
```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;
}

unittest {
Mat3   m;

m = 1.0f;
assert( m == 1.0f );
// The above assertion fails all the time.

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;
}

Mario
```
Mar 13 2007
```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;
}

unittest {
Mat3   m;

m = 1.0f;
assert( m == 1.0f );
// The above assertion fails all the time.
// I think that DMD interprets m 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  Vec3;

float dot( Vec3 v1, Vec3 v2 ) {
return v1 * v2 + v1 * v2 + v1 * v2;
}

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 * v2 + v1 * v2 + v1 * v2;
}

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?

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
```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
``` 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?

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);

Vector!(int, 2) c = a + b;
assert(c == 2);
assert(c == 5);

// check subtraction
Vector!(int, 2) d = a - b;
assert(d == 0);
assert(d == 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 - 8 < 0.001);
}
```
Mar 15 2007
```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?

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
```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

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
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
```
Mar 15 2007
```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
```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
```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
```;)

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
```[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.

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.

;)

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    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