patched shifted phase functions

TysonRayJones · TysonRayJones · commit 34fb5fc5cd7d · 2021-09-02T07:32:42.000+10:00
which slightly differed to documentation.

- Previously, SCALED_INVERSE_SHIFTED_DISTANCE computed coeff/sqrt( (x2-x1-dx)^2 + ... ), but now computes coeff/sqrt( (x1-x2-dx)^2 + ... (x1 and x2 have swapped)

- Previously, SCALED_INVERSE_SHIFTED_NORM and SCALED_INVERSE_SHIFTED_DISTANCE used their divergence parameter when their denominators were precisely zero. Now, the divergence parameter is used whenever the denominator is within REAL_EPS to zero.
diff --git a/QuEST/include/QuEST.h b/QuEST/include/QuEST.h
@@ -6176,6 +6176,11 @@ void applyNamedPhaseFuncOverrides(Qureg qureg, int* qubits, int* numQubitsPerReg
  *   \f[
  *      f(\vec{r}, \theta)|_{\theta=0.5} \; = \; \begin{cases} \pi & \;\;\; \vec{r}=\vec{0} \\ \displaystyle 0.5 \left[ \sum_j^{\text{numRegs}} {r_j}^2 \right]^{-1/2} & \;\;\;\text{otherwise} \end{cases}.
  *   \f] 
+ *   Notice the order of the parameters matches the order of the words in the \p phaseFunc.
+ *   > Functions \p SCALED_INVERSE_SHIFTED_NORM and \p SCALED_INVERSE_SHIFTED_DISTANCE,
+ *   > which can have denominators arbitrarily close to zero, will invoke the 
+ *   > divergence parameter whenever the denominator is smaller than (or equal to)
+ *   > machine precision `REAL_EPS`.
  *
  * - Functions allowing the shifting of sub-register values, which are \p SCALED_INVERSE_SHIFTED_NORM
  *   and \p SCALED_INVERSE_SHIFTED_DISTANCE, need these shift values to be passed in the \p params
diff --git a/QuEST/src/CPU/QuEST_cpu.c b/QuEST/src/CPU/QuEST_cpu.c
@@ -4529,15 +4529,15 @@ void statevec_applyParamNamedPhaseFuncOverrides(
                         for (r=0; r<numRegs; r++)
                             norm += phaseInds[r]*phaseInds[r];
                     norm = sqrt(norm);
-
+ 
                     if (phaseFuncName == NORM)
                         phase = norm;
                     else if (phaseFuncName == INVERSE_NORM)
-                        phase = (norm == 0.)? params[0] : 1/norm;
+                        phase = (norm == 0.)? params[0] : 1/norm;  // smallest non-zero norm is 1
                     else if (phaseFuncName == SCALED_NORM)
                         phase = params[0] * norm;
                     else if (phaseFuncName == SCALED_INVERSE_NORM || phaseFuncName == SCALED_INVERSE_SHIFTED_NORM)
-                        phase = (norm == 0.)? params[1] : params[0] / norm;
+                        phase = (norm <= REAL_EPS)? params[1] : params[0] / norm; // unless shifted closer to zero
                 }
                 // compute product related phases
                 else if (phaseFuncName == PRODUCT || phaseFuncName == INVERSE_PRODUCT ||
@@ -4550,7 +4550,7 @@ void statevec_applyParamNamedPhaseFuncOverrides(
                     if (phaseFuncName == PRODUCT)
                         phase = prod;
                     else if (phaseFuncName == INVERSE_PRODUCT)
-                        phase = (prod == 0.)? params[0] : 1/prod;
+                        phase = (prod == 0.)? params[0] : 1/prod;  // smallest non-zero product norm is +- 1
                     else if (phaseFuncName == SCALED_PRODUCT)
                         phase = params[0] * prod;
                     else if (phaseFuncName == SCALED_INVERSE_PRODUCT)
@@ -4564,7 +4564,7 @@ void statevec_applyParamNamedPhaseFuncOverrides(
                     dist = 0;
                     if (phaseFuncName == SCALED_INVERSE_SHIFTED_DISTANCE) {
                         for (r=0; r<numRegs; r+=2)
-                            dist += (phaseInds[r+1] - phaseInds[r] - params[2+r/2])*(phaseInds[r+1] - phaseInds[r] - params[2+r/2]);
+                            dist += (phaseInds[r] - phaseInds[r+1] - params[2+r/2])*(phaseInds[r] - phaseInds[r+1] - params[2+r/2]);
                     }
                     else
                         for (r=0; r<numRegs; r+=2)
@@ -4574,11 +4574,11 @@ void statevec_applyParamNamedPhaseFuncOverrides(
                     if (phaseFuncName == DISTANCE)
                         phase = dist;
                     else if (phaseFuncName == INVERSE_DISTANCE)
-                        phase = (dist == 0.)? params[0] : 1/dist;
+                        phase = (dist == 0.)? params[0] : 1/dist; // smallest non-zero dist is 1
                     else if (phaseFuncName == SCALED_DISTANCE)
                         phase = params[0] * dist;
                     else if (phaseFuncName == SCALED_INVERSE_DISTANCE || phaseFuncName == SCALED_INVERSE_SHIFTED_DISTANCE)
-                        phase = (dist == 0.)? params[1] : params[0] / dist;
+                        phase = (dist <= REAL_EPS)? params[1] : params[0] / dist; // unless shifted closer to 0
                 }
             }
             
diff --git a/QuEST/src/GPU/QuEST_gpu.cu b/QuEST/src/GPU/QuEST_gpu.cu
@@ -3835,11 +3835,11 @@ __global__ void statevec_applyParamNamedPhaseFuncOverridesKernel(
             if (phaseFuncName == NORM)
                 phase = norm;
             else if (phaseFuncName == INVERSE_NORM)
-                phase = (norm == 0.)? params[0] : 1/norm;
+                phase = (norm == 0.)? params[0] : 1/norm; // smallest non-zero norm is 1
             else if (phaseFuncName == SCALED_NORM)
                 phase = params[0] * norm;
             else if (phaseFuncName == SCALED_INVERSE_NORM || phaseFuncName == SCALED_INVERSE_SHIFTED_NORM)
-                phase = (norm == 0.)? params[1] : params[0] / norm;
+                phase = (norm <= REAL_EPS)? params[1] : params[0] / norm; // unless shifted closer to zero
         }
         // compute product related phases
         else if (phaseFuncName == PRODUCT || phaseFuncName == INVERSE_PRODUCT ||
@@ -3852,7 +3852,7 @@ __global__ void statevec_applyParamNamedPhaseFuncOverridesKernel(
             if (phaseFuncName == PRODUCT)
                 phase = prod;
             else if (phaseFuncName == INVERSE_PRODUCT)
-                phase = (prod == 0.)? params[0] : 1/prod;
+                phase = (prod == 0.)? params[0] : 1/prod; // smallest non-zero prod is +- 1
             else if (phaseFuncName == SCALED_PRODUCT)
                 phase = params[0] * prod;
             else if (phaseFuncName == SCALED_INVERSE_PRODUCT)
@@ -3866,7 +3866,7 @@ __global__ void statevec_applyParamNamedPhaseFuncOverridesKernel(
             qreal dist = 0;
             if (phaseFuncName == SCALED_INVERSE_SHIFTED_DISTANCE) {
                 for (int r=0; r<numRegs; r+=2) {
-                    qreal dif = (phaseInds[(r+1)*stride+offset] - phaseInds[r*stride+offset] - params[2+r/2]);
+                    qreal dif = (phaseInds[r*stride+offset] - phaseInds[(r+1)*stride+offset] - params[2+r/2]);
                     dist += dif*dif;
                 }
             }
@@ -3880,11 +3880,11 @@ __global__ void statevec_applyParamNamedPhaseFuncOverridesKernel(
             if (phaseFuncName == DISTANCE)
                 phase = dist;
             else if (phaseFuncName == INVERSE_DISTANCE)
-                phase = (dist == 0.)? params[0] : 1/dist;
+                phase = (dist == 0.)? params[0] : 1/dist; // smallest non-zero dist is 1
             else if (phaseFuncName == SCALED_DISTANCE)
                 phase = params[0] * dist;
             else if (phaseFuncName == SCALED_INVERSE_DISTANCE || phaseFuncName == SCALED_INVERSE_SHIFTED_DISTANCE)
-                phase = (dist == 0.)? params[1] : params[0] / dist;
+                phase = (dist <= REAL_EPS)? params[1] : params[0] / dist; // unless shifted closer
         }
     }
     
diff --git a/tests/test_operators.cpp b/tests/test_operators.cpp
@@ -1723,7 +1723,8 @@ TEST_CASE( "applyParamNamedPhaseFunc", "[operators]" ) {
                 qreal phase = 0;
                 for (int r=0; r<numRegs; r++)
                     phase += pow(regVals[i][r] - params[2+r], 2);
-                phase = (phase == 0.)? params[1] : params[0]/sqrt(phase);
+                phase = sqrt(phase);
+                phase = (phase <= REAL_EPS)? params[1] : params[0]/phase;
                 diagMatr[i][i] = expI(phase);
             }
             
@@ -1952,8 +1953,9 @@ TEST_CASE( "applyParamNamedPhaseFunc", "[operators]" ) {
                 for (size_t i=0; i<diagMatr.size(); i++) {
                     qreal phase = 0;
                     for (int r=0; r<numRegs; r+=2)
-                        phase += pow(regVals[i][r+1]-regVals[i][r]-params[2+r/2], 2);
-                    phase = (phase == 0.)? params[1] : params[0]/sqrt(phase);
+                        phase += pow(regVals[i][r]-regVals[i][r+1]-params[2+r/2], 2);
+                    phase = sqrt(phase);
+                    phase = (phase <= REAL_EPS)? params[1] : params[0]/phase;
                     diagMatr[i][i] = expI(phase);
                 }
             }
@@ -2271,7 +2273,8 @@ TEST_CASE( "applyParamNamedPhaseFuncOverrides", "[operators]" ) {
                 qreal phase = 0;
                 for (int r=0; r<numRegs; r++)
                     phase += pow(regVals[i][r] - params[2+r], 2);
-                phase = (phase == 0.)? params[1] : params[0]/sqrt(phase);
+                phase = sqrt(phase);
+                phase = (phase <= REAL_EPS)? params[1] : params[0]/phase;
                 diagMatr[i][i] = expI(phase);
             }
             setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
@@ -2505,8 +2508,9 @@ TEST_CASE( "applyParamNamedPhaseFuncOverrides", "[operators]" ) {
                 for (size_t i=0; i<diagMatr.size(); i++) {
                     qreal phase = 0;
                     for (int r=0; r<numRegs; r+=2)
-                        phase += pow(regVals[i][r+1]-regVals[i][r]-params[2+r/2], 2);
-                    phase = (phase == 0.)? params[1] : params[0]/sqrt(phase);
+                        phase += pow(regVals[i][r]-regVals[i][r+1]-params[2+r/2], 2);
+                    phase = sqrt(phase);
+                    phase = (phase <= REAL_EPS)? params[1] : params[0]/phase;
                     diagMatr[i][i] = expI(phase);
                 }
                 
