+/* ----------------------------------------------------------------------
+ * Small number of 64-bit integer arithmetic operations, to prevent
+ * integer overflow at the very core of cross().
+ */
+
+typedef struct {
+ long hi;
+ unsigned long lo;
+} int64;
+
+#define greater64(i,j) ( (i).hi>(j).hi || ((i).hi==(j).hi && (i).lo>(j).lo))
+#define sign64(i) ((i).hi < 0 ? -1 : (i).hi==0 && (i).lo==0 ? 0 : +1)
+
+int64 mulu32to64(unsigned long x, unsigned long y)
+{
+ unsigned long a, b, c, d, t;
+ int64 ret;
+
+ a = (x & 0xFFFF) * (y & 0xFFFF);
+ b = (x & 0xFFFF) * (y >> 16);
+ c = (x >> 16) * (y & 0xFFFF);
+ d = (x >> 16) * (y >> 16);
+
+ ret.lo = a;
+ ret.hi = d + (b >> 16) + (c >> 16);
+ t = (b & 0xFFFF) << 16;
+ ret.lo += t;
+ if (ret.lo < t)
+ ret.hi++;
+ t = (c & 0xFFFF) << 16;
+ ret.lo += t;
+ if (ret.lo < t)
+ ret.hi++;
+
+#ifdef DIAGNOSTIC_VIA_LONGLONG
+ assert(((unsigned long long)ret.hi << 32) + ret.lo ==
+ (unsigned long long)x * y);
+#endif
+
+ return ret;
+}
+
+int64 mul32to64(long x, long y)
+{
+ int sign = +1;
+ int64 ret;
+#ifdef DIAGNOSTIC_VIA_LONGLONG
+ long long realret = (long long)x * y;
+#endif
+
+ if (x < 0)
+ x = -x, sign = -sign;
+ if (y < 0)
+ y = -y, sign = -sign;
+
+ ret = mulu32to64(x, y);
+
+ if (sign < 0) {
+ ret.hi = -ret.hi;
+ ret.lo = -ret.lo;
+ if (ret.lo)
+ ret.hi--;
+ }
+
+#ifdef DIAGNOSTIC_VIA_LONGLONG
+ assert(((unsigned long long)ret.hi << 32) + ret.lo == realret);
+#endif
+
+ return ret;
+}
+
+int64 dotprod64(long a, long b, long p, long q)
+{
+ int64 ab, pq;
+
+ ab = mul32to64(a, b);
+ pq = mul32to64(p, q);
+ ab.hi += pq.hi;
+ ab.lo += pq.lo;
+ if (ab.lo < pq.lo)
+ ab.hi++;
+ return ab;
+}
+