Geometric mean for total work calc

src/bignum.h

@@ -136,7 +136,7 @@ public:
 
         if (sn < (int64_t)0)
         {
-            // Since the minimum signed integer cannot be represented as positive so long as its type is signed, 
+            // Since the minimum signed integer cannot be represented as positive so long as its type is signed,
             // and it's not well-defined what happens if you make it unsigned before negating it,
             // we instead increment the negative integer by 1, convert it, then increment the (now positive) unsigned integer by 1 to compensate
             n = -(sn + 1);
@@ -516,6 +516,77 @@ public:
     friend inline const CBigNum operator-(const CBigNum& a, const CBigNum& b);
     friend inline const CBigNum operator/(const CBigNum& a, const CBigNum& b);
     friend inline const CBigNum operator%(const CBigNum& a, const CBigNum& b);
+
+    CBigNum nthRoot(int n) const
+    {
+        assert(n > 1);
+        if (BN_is_zero(this))
+            return 0;
+        assert(!BN_is_negative(this));
+
+        // starting approximation
+        int nRootBits = (BN_num_bits(this) + n - 1) / n;
+        int nStartingBits = std::min(8, nRootBits);
+        CBigNum bnUpper = *this;
+        bnUpper >>= (nRootBits - nStartingBits)*n;
+        CBigNum bnCur = 0;
+        for (int i = nStartingBits - 1; i >= 0; i--)
+        {
+            CBigNum bnNext = bnCur;
+            bnNext += 1 << i;
+            CBigNum bnPower(1);
+            for (int j = 0; j < n; j++)
+                bnPower *= bnNext;
+            if (BN_cmp(&bnPower, &bnUpper) <= 0)
+                bnCur = bnNext;
+        }
+        if (nRootBits == nStartingBits)
+            return bnCur;
+        bnCur <<= nRootBits - nStartingBits;
+
+        // iterate: cur = cur + (*this / cur^^(n-1) - cur)/n
+        CBigNum bnDelta;
+        const CBigNum bnRoot(n);
+        int nTerminate = 0;
+        // this should always converge in fewer steps, but limit just in case
+        for (int it = 0; it < 20; it++)
+        {
+            CBigNum bnDenominator = 1;
+            for (int i = 0; i < n - 1; i++)
+                bnDenominator *= bnCur;
+            bnDelta = *this / bnDenominator - bnCur;
+            if (BN_is_zero(&bnDelta))
+                return bnCur;
+            if (BN_is_negative(&bnDelta))
+            {
+                if (nTerminate == 1)
+                    return bnCur - 1;
+                BN_set_negative(&bnDelta, 0);
+                if (BN_cmp(&bnDelta, &bnRoot) <= 0)
+                {
+                    bnCur -= 1;
+                    nTerminate = -1;
+                    continue;
+                }
+                BN_set_negative(&bnDelta, 1);
+            }
+            else
+            {
+                if (nTerminate == -1)
+                    return bnCur;
+                if (BN_cmp(&bnDelta, &bnRoot) <= 0)
+                {
+                    bnCur += 1;
+                    nTerminate = 1;
+                    continue;
+                }
+            }
+            bnCur += bnDelta / n;
+            nTerminate = 0;
+        }
+        return bnCur;
+    }
+
 };
 
 

src/core.h

@@ -47,6 +47,7 @@ static const int AUXPOW_START_MAINNET = 9999999; //TODO change me
 static const int AUXPOW_START_TESTNET = 1;
 
 static const int Phase2Reward_Start = 9999999;
+static const int GeoAvgWork_Start = 9999999;
 
 inline int GetAlgo(int nVersion)
 {

src/main.h

@@ -965,6 +965,26 @@ public:
     {
         CBigNum bnRes;
         if ((TestNet() && (nHeight > 500)) ||
+            (nHeight >= GeoAvgWork_Start))
+        {
+            CBigNum nBlockWork = GetBlockWork();
+            int nAlgo = GetAlgo();
+            for (int algo = 0; algo < NUM_ALGOS; algo++)
+            {
+                if (algo != nAlgo)
+                {
+                    CBigNum nBlockWorkAlt = GetPrevWorkForAlgoWithDecay2(algo);
+                    if (nBlockWorkAlt != 0)
+                        nBlockWork *= nBlockWorkAlt;
+                }
+            }
+            // Compute the geometric mean
+            bnRes = bnRes.nthRoot(NUM_ALGOS);
+            // Scale to roughly match the old work calculation
+            bnRes <<= 8;
+        }
+        else
+        if (TestNet() ||
             (nHeight >= nBlockAlgoNormalisedWorkStart))
         {
             // Adjusted block work is the sum of work of this block and the
@@ -975,7 +995,7 @@ public:
             {
                 if (algo != nAlgo)
                 {
-                    if (nHeight >= nBlockAlgoNormalisedWorkDecayStart2)
+                    if (TestNet() || (nHeight >= nBlockAlgoNormalisedWorkDecayStart2))
                         nBlockWork += GetPrevWorkForAlgoWithDecay2(algo);
                     else
                         if (nHeight >= nBlockAlgoNormalisedWorkDecayStart)
@@ -1000,7 +1020,7 @@ public:
 	/*        int algo = GetAlgo();
         if (algo == ALGO_SHA256D)
             return CheckProofOfWork(GetBlockHash(), nBits, algo);
-        else*/ 
+        else*/
             return true;
     }