[mathicgb] 73/393: I made the hash-table verification of a hit process two exponents at a time by comparing them as an int64 instead of two int32's. Turns out the strict aliasing rule makes this very tricky to get right, and, after a lot of effort, I finally found out that the only correct solution that also achieves good performance uses... memcpy. Yes, on individual int64 values. Turns out at least MSVC and GCC understand memcpy and won't actually call a memcpy function for something like that - it is just a way to get them to understand what's going on without the undefined-behavior-inducing cast. Strange but true.

Doug Torrance dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:34 UTC 2015


This is an automated email from the git hooks/post-receive script.

dtorrance-guest pushed a commit to branch upstream
in repository mathicgb.

commit 282e2f87246f48ba48cb47be62bc3bd67e436854
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date:   Tue Oct 23 21:43:53 2012 +0200

    I made the hash-table verification of a hit process two exponents at a time by comparing them as an int64 instead of two int32's. Turns out the strict aliasing rule makes this very tricky to get right, and, after a lot of effort, I finally found out that the only correct solution that also achieves good performance uses... memcpy. Yes, on individual int64 values. Turns out at least MSVC and GCC understand memcpy and won't actually call a memcpy function for something like that - it is [...]
---
 src/mathicgb/MonomialMap.hpp | 18 +++++++++-------
 src/mathicgb/PolyRing.cpp    |  1 +
 src/mathicgb/PolyRing.hpp    | 51 +++++++++++++++++++++++++++++++++-----------
 3 files changed, 49 insertions(+), 21 deletions(-)

diff --git a/src/mathicgb/MonomialMap.hpp b/src/mathicgb/MonomialMap.hpp
index 6b22b2f..f073e0f 100755
--- a/src/mathicgb/MonomialMap.hpp
+++ b/src/mathicgb/MonomialMap.hpp
@@ -70,11 +70,12 @@ namespace MonomialMapInternal {
       const HashValue monoHash = mRing.monomialHashValue(mono);
       Node* node = entry(hashToIndex(monoHash));
       for (; node != 0; node = node->next) {
-        if (monoHash == mRing.monomialHashValue(node->mono) &&
-          mRing.monomialEqualHintTrue(mono, node->mono)
-        ) {
+        // To my surprise, it seems to be faster to comment out this branch.
+        // I guess the hash table has too few collisions to make it worth it.
+        //if (monoHash != mRing.monomialHashValue(node->mono))
+        //  continue;
+        if (mRing.monomialEqualHintTrue(mono, node->mono))
           return &node->value;
-        }
       }
       return 0;
     }
@@ -83,11 +84,12 @@ namespace MonomialMapInternal {
       const HashValue abHash = mRing.monomialHashOfProduct(a, b);
       Node* node = entry(hashToIndex(abHash));
       for (; node != 0; node = node->next) {
-        if (abHash == mRing.monomialHashValue(node->mono) &&
-          mRing.monomialIsProductOfHintTrue(a, b, node->mono)
-        ) {
+        // To my surprise, it seems to be faster to comment out this branch.
+        // I guess the hash table has too few collisions to make it worth it.
+        //if (abHash != mRing.monomialHashValue(node->mono))
+        //  continue;
+        if (mRing.monomialIsProductOfHintTrue(a, b, node->mono))
           return &node->value;
-        }
       }
       return 0;
     }
diff --git a/src/mathicgb/PolyRing.cpp b/src/mathicgb/PolyRing.cpp
index b6eecda..2d89dcd 100755
--- a/src/mathicgb/PolyRing.cpp
+++ b/src/mathicgb/PolyRing.cpp
@@ -29,6 +29,7 @@ PolyRing::PolyRing(coefficient p0,
     mTotalDegreeGradedOnly(false)
 {
   resetCoefficientStats();
+  srand(0);
   for (size_t i=0; i<mNumVars; i++)
     mHashVals.push_back(static_cast<HashValue>(rand()));
 }
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index d683988..be6c4fa 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -10,6 +10,7 @@
 #include <vector>
 #include <memtailor.h>
 #include <cstdio>
+#include <cstring>
 
 #define LT (-1)
 #define EQ 0
@@ -117,7 +118,7 @@ class ConstMonomial
 public:
   //* Wrap a raw pointer to create a monomial
   // Assumptions:
-  //  1. This is done in the presence of an PolyRing
+  //  1. This is done in the presence of a PolyRing
   //  2. Space for the monomial has been created
   ConstMonomial() : mValue(0) {}
   ConstMonomial(const exponent *val) : mValue(val) {}
@@ -132,9 +133,9 @@ public:
   exponent component() const { return *mValue; }
 
 private:
-  exponent operator[](size_t index) const { return mValue[index]; }
+  const exponent& operator[](size_t index) const { return mValue[index]; }
+  const exponent& operator*() const { return *mValue; }
 
-  exponent operator*() const { return *mValue; }
 protected:
   exponent const* mValue;
 };
@@ -150,7 +151,7 @@ class Monomial : public ConstMonomial
 public:
   //* Wrap a raw pointer to create a monomial
   // Assumptions:
-  //  1. This is done in the presence of an PolyRing
+  //  1. This is done in the presence of a PolyRing
   //  2. Space for the monomial has been created
   Monomial() : ConstMonomial() {}
   Monomial(exponent *val) : ConstMonomial(val) {}
@@ -163,10 +164,10 @@ public:
   exponent const * unsafeGetRepresentation() const { return mValue; }
 
 private:
-  exponent operator[](size_t index) const { return mValue[index]; }
+  const exponent& operator[](size_t index) const { return mValue[index]; }
   exponent& operator[](size_t index) { return unsafeGetRepresentation()[index]; }
 
-  exponent operator*() const { return *mValue; }
+  const exponent& operator*() const { return *mValue; }
   exponent& operator*() { return * const_cast<exponent *>(mValue); }
 };
 
@@ -638,13 +639,37 @@ inline bool PolyRing::monomialIsProductOfHintTrue(
   const ConstMonomial b, 
   const ConstMonomial ab
 ) const {
-  // same idea as monomialEqualHintTrue, just applied to the sum of a and b.
-  exponent orOfXor = 0;
-  for (size_t i = mNumVars; i != 0; --i)
-    orOfXor |= ab[i] ^ (a[i] + b[i]);
-  const bool isProduct = (orOfXor == 0);
-  MATHICGB_ASSERT(isProduct == monomialIsProductOf(a, b, ab));
-  return isProduct;
+  // We compare more than one exponent at a time using 64 bit integers. This 
+  // might go one 32 bit value at the end too far, but since that space is
+  // either a degree or a hash value that is fine --- those values will also
+  // match if the monomials are equal. This does not work for negative
+  // exponents since the overflowing bit will go into the next word.
+  // It is OK that the degree field can be negative (a field we might go
+  // into without caring about it because it shares a 64 bit field with
+  // the last exponent), because it is at the end so the overflowing
+  // bit will not interfere.
+
+  // todo: ensure 8 byte alignment. Though there seem to be no ill effects
+  // for unaligned access. Performance seems to be no worse than for using
+  // 32 bit integers directly.
+
+  uint64 orOfXor = 0;
+  for (size_t i = mNumVars / 2; i != static_cast<size_t>(-1); --i) {
+    uint64 A, B, AB;
+    // We have to use std::memcpy here because just casting to a int64 breaks
+    // the strict aliasing rule which implies undefined behavior. Both MSVC and
+    // gcc don't actually call memcpy here. MSVC is a tiny big slower for this
+    // code than for casting while GCC seems to be exactly the same speed.
+    std::memcpy(&A, &a[i], 8);
+    std::memcpy(&B, &b[i], 8);
+    std::memcpy(&AB, &ab[i], 8);
+    orOfXor |= AB ^ (A + B);
+  }
+  MATHICGB_ASSERT((orOfXor == 0) == monomialIsProductOf(a, b, ab))
+
+
+  return orOfXor == 0;
+  
 }
 
 inline bool PolyRing::monomialIsProductOf(

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/mathicgb.git



More information about the debian-science-commits mailing list