[Issue 5880] New: std.numerics.dotProduct for fixed-size arrays
d-bugmail at puremagic.com
d-bugmail at puremagic.com
Sun Apr 24 11:25:16 PDT 2011
http://d.puremagic.com/issues/show_bug.cgi?id=5880
Summary: std.numerics.dotProduct for fixed-size arrays
Product: D
Version: D2
Platform: All
OS/Version: All
Status: NEW
Severity: enhancement
Priority: P2
Component: Phobos
AssignedTo: nobody at puremagic.com
ReportedBy: bearophile_hugs at eml.cc
--- Comment #0 from bearophile_hugs at eml.cc 2011-04-24 11:21:34 PDT ---
A third overload for fixed-sized arrays offers:
- compile-time errors for the length mismatch instead of run-time ones (not in
release build);
- allows compilers to optimize the code better because the lengths are known at
compile-time;
- the fixed-size argument arrays are by reference, avoiding the copy in this
case too.
/*pure*/ CommonType!(ElementType!(Range1), ElementType!(Range2))
dotProduct(Range1, Range2)(Range1 a, Range2 b)
if (isInputRange!(Range1) && isInputRange!(Range2) &&
!(isArray!(Range1) && isArray!(Range2)))
{
// can't be pure yet because of length property and enforce
enum bool haveLen = hasLength!(Range1) && hasLength!(Range2);
static if (haveLen) enforce(a.length == b.length);
typeof(return) result = 0;
for (; !a.empty; a.popFront, b.popFront)
{
result += a.front * b.front;
}
static if (!haveLen) enforce(b.empty);
return result;
}
/// Ditto
pure Unqual!(CommonType!(F1, F2))
dotProduct(F1, F2)(in F1[] avector, in F2[] bvector)
if (!isStaticArray!F1 || !isStaticArray!F2)
{
immutable n = avector.length;
assert(n == bvector.length);
auto avec = avector.ptr, bvec = bvector.ptr;
typeof(return) sum0 = 0, sum1 = 0;
const all_endp = avec + n;
const smallblock_endp = avec + (n & ~3);
const bigblock_endp = avec + (n & ~15);
for (; avec != bigblock_endp; avec += 16, bvec += 16)
{
sum0 += avec[0] * bvec[0];
sum1 += avec[1] * bvec[1];
sum0 += avec[2] * bvec[2];
sum1 += avec[3] * bvec[3];
sum0 += avec[4] * bvec[4];
sum1 += avec[5] * bvec[5];
sum0 += avec[6] * bvec[6];
sum1 += avec[7] * bvec[7];
sum0 += avec[8] * bvec[8];
sum1 += avec[9] * bvec[9];
sum0 += avec[10] * bvec[10];
sum1 += avec[11] * bvec[11];
sum0 += avec[12] * bvec[12];
sum1 += avec[13] * bvec[13];
sum0 += avec[14] * bvec[14];
sum1 += avec[15] * bvec[15];
}
for (; avec != smallblock_endp; avec += 4, bvec += 4) {
sum0 += avec[0] * bvec[0];
sum1 += avec[1] * bvec[1];
sum0 += avec[2] * bvec[2];
sum1 += avec[3] * bvec[3];
}
sum0 += sum1;
/* Do trailing portion in naive loop. */
while (avec != all_endp)
sum0 += (*avec++) * (*bvec++);
return sum0;
}
/// Ditto
pure Unqual!(CommonType!(F1, F2))
dotProduct(F1, F2, size_t n, size_t n2)(ref const F1[n] avector, ref const
F2[n2] bvector)
if (isStaticArray!(typeof(avector)) && isStaticArray!(typeof(bvector)))
{
static assert(n == n2); // do not move this to the template constraints
auto avec = avector.ptr, bvec = bvector.ptr;
typeof(return) sum0 = 0, sum1 = 0;
const all_endp = avec + n;
const smallblock_endp = avec + (n & ~3);
const bigblock_endp = avec + (n & ~15);
for (; avec != bigblock_endp; avec += 16, bvec += 16)
{
sum0 += avec[0] * bvec[0];
sum1 += avec[1] * bvec[1];
sum0 += avec[2] * bvec[2];
sum1 += avec[3] * bvec[3];
sum0 += avec[4] * bvec[4];
sum1 += avec[5] * bvec[5];
sum0 += avec[6] * bvec[6];
sum1 += avec[7] * bvec[7];
sum0 += avec[8] * bvec[8];
sum1 += avec[9] * bvec[9];
sum0 += avec[10] * bvec[10];
sum1 += avec[11] * bvec[11];
sum0 += avec[12] * bvec[12];
sum1 += avec[13] * bvec[13];
sum0 += avec[14] * bvec[14];
sum1 += avec[15] * bvec[15];
}
for (; avec != smallblock_endp; avec += 4, bvec += 4) {
sum0 += avec[0] * bvec[0];
sum1 += avec[1] * bvec[1];
sum0 += avec[2] * bvec[2];
sum1 += avec[3] * bvec[3];
}
sum0 += sum1;
/* Do trailing portion in naive loop. */
while (avec != all_endp)
sum0 += (*avec++) * (*bvec++);
return sum0;
}
unittest
{
minidot()
double[] a0, b0;
assert(dotProduct(a0, b0) == 0);
assert(dotProduct([1.0, 2.0], [4.0, 6.0]) == 16.0);
assert(dotProduct([1, 3, -5], [4, -2, -1]) == 3);
assert(dotProduct(iota(1, 5), iota(10, 41, 10)) == 300);
int[4] a1 = [1, 2, 3, 4];
int[4] b1 = [10, 20, 30, 40];
assert(dotProduct(a1, b1) == 300);
int[] c1 = [10, 20, 30, 40];
assert(dotProduct(a1, c1) == 300);
int[5] c2 = [10, 20, 30, 40, 0];
assert(!__traits(compiles, { dotProduct(a1, c2); } )); // can't compile
// more unittests needed
}
--
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
More information about the Digitalmars-d-bugs
mailing list