diff --git a/quest/src/api/calculations.cpp b/quest/src/api/calculations.cpp
index 47e5d8a63..e18227cd5 100644
--- a/quest/src/api/calculations.cpp
+++ b/quest/src/api/calculations.cpp
@@ -428,6 +428,11 @@ qreal calcProbOfMultiQubitOutcome(Qureg qureg, vector<int> qubits, vector<int> o
 
 vector<qreal> calcProbsOfAllMultiQubitOutcomes(Qureg qureg, vector<int> qubits) {
 
+    // validate before exponentiating the caller-controlled vector size
+    validate_quregFields(qureg, __func__);
+    validate_targets(qureg, qubits.data(), qubits.size(), __func__);
+    validate_measurementOutcomesFitInGpuMem(qureg, qubits.size(), __func__);
+
     // allocate temp vector, and validate successful (since it's exponentially large!)
     vector<qreal> out;
     qindex numOut = powerOf2(qubits.size());
diff --git a/quest/src/api/channels.cpp b/quest/src/api/channels.cpp
index c6702438a..72aa22772 100644
--- a/quest/src/api/channels.cpp
+++ b/quest/src/api/channels.cpp
@@ -108,7 +108,7 @@ void freeAllMemoryIfAnyAllocsFailed(T& obj) {
     // determine whether any node experienced a failure
     bool anyFail = didAnyLocalAllocsFail(obj);
     if (comm_isActive())
-        anyFail = comm_isTrueOnAllNodes(anyFail);
+        anyFail = comm_isTrueOnAnyNode(anyFail);
 
     // if so, free all memory before subsequent validation
     if (anyFail)
diff --git a/quest/src/api/decoherence.cpp b/quest/src/api/decoherence.cpp
index 4e1901f25..3bd393e20 100644
--- a/quest/src/api/decoherence.cpp
+++ b/quest/src/api/decoherence.cpp
@@ -143,6 +143,7 @@ void mixQureg(Qureg outQureg, Qureg inQureg, qreal inProb) {
 
 void mixSuperOp(Qureg qureg, int* targets, int numTargets, SuperOp superop) {
     validate_quregFields(qureg, __func__);
+    validate_quregIsDensityMatrix(qureg, __func__);
     validate_targets(qureg, targets, numTargets, __func__);
     validate_superOpFields(superop, __func__);
     validate_superOpIsSynced(superop, __func__);
diff --git a/quest/src/api/environment.cpp b/quest/src/api/environment.cpp
index e6fd2903a..37c523646 100644
--- a/quest/src/api/environment.cpp
+++ b/quest/src/api/environment.cpp
@@ -375,7 +375,7 @@ void printQuregAutoDeployments(bool isDensMatr) {
     // auto-deploy will still try to deploy there to (then subsequent validation will fail)
     int maxQubits = mem_getMaxNumQuregQubitsBeforeGlobalMemSizeofOverflow(isDensMatr, global_envPtr->numNodes);
 
-    for (int numQubits=1; numQubits<maxQubits; numQubits++) {
+    for (int numQubits=1; numQubits<=maxQubits; numQubits++) {
 
         // re-choose auto deployment
         useDistrib  = modeflag::USE_AUTO;
diff --git a/quest/src/api/initialisations.cpp b/quest/src/api/initialisations.cpp
index 36f63910c..a8be7ae3e 100644
--- a/quest/src/api/initialisations.cpp
+++ b/quest/src/api/initialisations.cpp
@@ -293,7 +293,8 @@ void setDensityQuregAmps(Qureg qureg, qindex startRow, qindex startCol, vector<v
         ptrs[i] = amps[i].data();
 
     // C function performs main validation
-    setDensityQuregAmps(qureg, startRow, startCol, ptrs.data(), len, (len>0)? len : 0); // avoid seg-fault
+    size_t numCols = (len > 0)? amps[0].size() : 0; // avoid seg-fault
+    setDensityQuregAmps(qureg, startRow, startCol, ptrs.data(), len, numCols);
 }
 
 void setQuregAmps(Qureg qureg, qindex startInd, vector<qcomp> amps) {
diff --git a/quest/src/api/matrices.cpp b/quest/src/api/matrices.cpp
index 07e37025c..10a889341 100644
--- a/quest/src/api/matrices.cpp
+++ b/quest/src/api/matrices.cpp
@@ -166,7 +166,7 @@ void freeAllMemoryIfAnyAllocsFailed(T matr) {
     // ascertain whether any allocs failed on any node
     bool anyFail = didAnyLocalAllocsFail(matr);
     if (comm_isActive())
-        anyFail = comm_isTrueOnAllNodes(anyFail);
+        anyFail = comm_isTrueOnAnyNode(anyFail);
 
     // if so, free all heap fields
     if (anyFail)
@@ -310,7 +310,7 @@ FullStateDiagMatr validateAndCreateCustomFullStateDiagMatr(int numQubits, int us
     out.cpuElems = cpu_allocArray(numElemsPerNode);
     out.gpuElems = (useGpuAccel)? gpu_allocArray(numElemsPerNode) : nullptr;
 
-    validateMatrixAllocs(out, __func__);
+    validateMatrixAllocs(out, caller);
     setInitialHeapFlags(out);
     return out;
 }
@@ -743,7 +743,7 @@ extern "C" void setFullStateDiagMatrFromMultiDimLists(FullStateDiagMatr out, voi
 template<class T> 
 void validateAndPrintMatrix(T matr, const char* caller) {
     validate_matrixFields(matr, caller);
-    validate_numReportedNewlinesAboveZero(__func__); // because trailing newline mandatory
+    validate_numReportedNewlinesAboveZero(caller); // because trailing newline mandatory
 
     // syncable matrices must be synced before reporting (though only CPU elems are printed)
     if constexpr (util_isHeapMatrixType<T>())
diff --git a/quest/src/api/paulis.cpp b/quest/src/api/paulis.cpp
index d5ac4d8e7..4875bb542 100644
--- a/quest/src/api/paulis.cpp
+++ b/quest/src/api/paulis.cpp
@@ -39,7 +39,7 @@ bool didAnyAllocsFailOnAnyNode(PauliStrSum sum) {
         ! mem_isAllocated(sum.isApproxHermitian) );
     
     if (comm_isActive())
-        anyFail = comm_isTrueOnAllNodes(anyFail);
+        anyFail = comm_isTrueOnAnyNode(anyFail);
 
     return anyFail;
 }
diff --git a/quest/src/api/qureg.cpp b/quest/src/api/qureg.cpp
index 683ca183d..57de85c0b 100644
--- a/quest/src/api/qureg.cpp
+++ b/quest/src/api/qureg.cpp
@@ -117,7 +117,7 @@ bool didAnyAllocsFailOnAnyNode(Qureg qureg) {
 
     bool anyFail = didAnyLocalAllocsFail(qureg);
     if (comm_isActive())
-        anyFail = comm_isTrueOnAllNodes(anyFail);
+        anyFail = comm_isTrueOnAnyNode(anyFail);
 
     return anyFail;
 }
@@ -130,7 +130,7 @@ void freeAllMemoryIfAnyAllocsFailed(Qureg qureg) {
         return;
 
     // otherwise, free everything that was successfully allocated (freeing nullptr is legal)
-    cpu_deallocArray(qureg.cpuAmps);
+    cpu_deallocNumaArray(qureg.cpuAmps, qureg.numAmpsPerNode);
     cpu_deallocArray(qureg.cpuCommBuffer);
 
     // although we avoid calling GPU deallocation in non-GPU mode
@@ -166,7 +166,7 @@ Qureg validateAndCreateCustomQureg(int numQubits, int isDensMatr, int useDistrib
 
     // if any of the above mallocs failed, below validation will memory leak; so free first (but don't set to nullptr)
     freeAllMemoryIfAnyAllocsFailed(qureg);
-    validate_newQuregAllocs(qureg, __func__);
+    validate_newQuregAllocs(qureg, caller);
 
     // TODO / DEBUG:
     // TEMPORARILY FORBID HIP GPU ACCELERATION OF > 2^32 AMPS PER NODE,
@@ -534,6 +534,11 @@ extern "C" void _wrap_getDensityQuregAmp(qcomp* out, Qureg qureg, qindex row, qi
 
 vector<qcomp> getQuregAmps(Qureg qureg, qindex startInd, qindex numAmps) {
 
+    // validate before interpreting numAmps as an allocation size
+    validate_quregFields(qureg, __func__);
+    validate_quregIsStateVector(qureg, __func__);
+    validate_basisStateIndices(qureg, startInd, numAmps, __func__);
+
     // allocate the output vector, and validate successful
     vector<qcomp> out;
     auto callback = [&]() { validate_tempListAllocSucceeded(false, numAmps, sizeof(qcomp), __func__); };
@@ -547,9 +552,14 @@ vector<qcomp> getQuregAmps(Qureg qureg, qindex startInd, qindex numAmps) {
 
 vector<vector<qcomp>> getDensityQuregAmps(Qureg qureg, qindex startRow, qindex startCol, qindex numRows, qindex numCols) {
 
+    // validate before multiplying dimensions or interpreting them as allocation sizes
+    validate_quregFields(qureg, __func__);
+    validate_quregIsDensityMatrix(qureg, __func__);
+    validate_basisStateRowCols(qureg, startRow, startCol, numRows, numCols, __func__);
+
     // allocate the output matrix, and validate successful
     vector<vector<qcomp>> out;
-    qindex numElems = numRows * numCols; // never overflows (else Qureg alloc would fail)
+    qindex numElems = numRows * numCols;
     auto callback1 = [&]() { validate_tempListAllocSucceeded(false, numElems, sizeof(qcomp), __func__); };
     util_tryAllocMatrix(out, numRows, numCols, callback1);
 
diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp
index 6fd5781ba..56c09792b 100644
--- a/quest/src/api/trotterisation.cpp
+++ b/quest/src/api/trotterisation.cpp
@@ -36,7 +36,7 @@ void internal_applyFirstOrderTrotterRepetition(
 ) {
     // apply each sum term as a gadget, in forward or reverse order
     for (qindex i=0; i<sum.numTerms; i++) {
-        int j = reverse? sum.numTerms - i - 1 : i;
+        qindex j = reverse? sum.numTerms - i - 1 : i;
         qindex k = (sumOrdering.empty())? j : sumOrdering[j];
         qcomp coeff = sum.coeffs[k];
         PauliStr str = sum.strings[k];
diff --git a/quest/src/comm/comm_routines.cpp b/quest/src/comm/comm_routines.cpp
index cf6956454..900ece98c 100644
--- a/quest/src/comm/comm_routines.cpp
+++ b/quest/src/comm/comm_routines.cpp
@@ -168,42 +168,54 @@ int getMaxNumMessages() {
 }
 
 
-std::array<qindex,2> dividePow2PayloadIntoMessages(qindex numAmps) {
-    assert_commPayloadIsPowerOf2(numAmps);
+std::array<qindex,2> dividePow2PayloadIntoMessages(qindex messageLength) {
+    assert_commPayloadIsPowerOf2(messageLength);
 
     // use single message if possible
-    if (numAmps < MAX_MESSAGE_LENGTH)
-        return {numAmps, 1};
+    if (messageLength < MAX_MESSAGE_LENGTH) {
+        assert_commMessageParamFitsIntoMpiInt(messageLength);
+        return {messageLength, 1};
+    }
     
-    // else, payload divides evenly between max-size messages (always fits in int)
-    qindex numMessages = numAmps / MAX_MESSAGE_LENGTH;
+    // else, payload divides evenly between max-size messages
+    qindex numMessages = messageLength / MAX_MESSAGE_LENGTH;
 
     // which we must be able to uniquely tag
     if (numMessages > getMaxNumMessages())
         error_commNumMessagesExceedTagMax();
 
-    // outputs always fit in 'int' but we force them to be 'qindex' since 
-    // caller will multiply them with ints and could easily overflow
+    // ensure all outputs can safely narrow at the MPI boundary
+    // (they're kept qindex to avoid overflows in algebra)
+    assert_commMessageParamFitsIntoMpiInt(MAX_MESSAGE_LENGTH);
+    assert_commMessageParamFitsIntoMpiInt(numMessages);
+
     return {MAX_MESSAGE_LENGTH, numMessages};
 }
 
 
-std::array<qindex,3> dividePayloadIntoMessages(qindex numAmps) {
+std::array<qindex,3> dividePayloadIntoMessages(qindex messageLength) {
 
     // use single message if possible
-    if (numAmps < MAX_MESSAGE_LENGTH)
-        return {numAmps, 1, 0};
+    if (messageLength < MAX_MESSAGE_LENGTH) {
+        assert_commMessageParamFitsIntoMpiInt(messageLength);
+        return {messageLength, 1, 0};
+    }
 
     // else, use as many max-size messages as possible, and one smaller msg
-    qindex numMaxSizeMsgs   = numAmps / MAX_MESSAGE_LENGTH; // floors
-    qindex remainingMsgSize = numAmps - numMaxSizeMsgs * MAX_MESSAGE_LENGTH;
+    qindex numMaxSizeMsgs   = messageLength / MAX_MESSAGE_LENGTH; // floors
+    qindex remainingMsgSize = messageLength - numMaxSizeMsgs * MAX_MESSAGE_LENGTH;
 
     // all of which we must be able to uniquely tag
-    if (numMaxSizeMsgs + 1 > getMaxNumMessages())
+    qindex numMessages = numMaxSizeMsgs + (remainingMsgSize > 0);
+    if (numMessages > getMaxNumMessages())
         error_commNumMessagesExceedTagMax();
 
-    // outputs always fit in 'int' but we force them to be 'qindex' since 
-    // caller will multiply them with ints and could easily overflow
+    // ensure all outputs can safely narrow at the MPI boundary
+    // (they're kept qindex to avoid overflows in algebra)
+    assert_commMessageParamFitsIntoMpiInt(MAX_MESSAGE_LENGTH);
+    assert_commMessageParamFitsIntoMpiInt(numMaxSizeMsgs);
+    assert_commMessageParamFitsIntoMpiInt(remainingMsgSize);
+
     return {MAX_MESSAGE_LENGTH, numMaxSizeMsgs, remainingMsgSize};
 }
 
@@ -221,20 +233,21 @@ void exchangeArrays(qcomp* send, qcomp* recv, qindex numElems, int pairRank) {
 
     // each message is asynchronously dispatched with a final wait, as per arxiv.org/abs/2308.07402
 
-    // we will send payload in multiple asynch messages (create two requests per msg for subsequent synch)
+    // we will send payload in multiple asynch messages
     auto [messageSize, numMessages] = dividePow2PayloadIntoMessages(numElems);
-    vector<MPI_Request> requests(2*numMessages, MPI_REQUEST_NULL);
+    vector<MPI_Request> recvRequests(numMessages, MPI_REQUEST_NULL);
+    vector<MPI_Request> sendRequests(numMessages, MPI_REQUEST_NULL);
 
     // asynchronously exchange the messages (effecting MPI_Isendrecv), using unique tags 
     // so that messages are permitted to arrive out-of-order (supporting UCX adaptive-routing)
     for (qindex m=0; m<numMessages; m++) {
-        int tag = static_cast<int>(m); // gauranteed int, but m*messageSize needs qindex
-        MPI_Irecv(&recv[m*messageSize], messageSize, MPI_QCOMP, pairRank, tag, mpiComm, &requests[2*m]);
-        MPI_Isend(&send[m*messageSize], messageSize, MPI_QCOMP, pairRank, tag, mpiComm, &requests[2*m+1]);
+        MPI_Irecv(&recv[m * messageSize], messageSize, MPI_QCOMP, pairRank, /*tag=*/m, mpiComm, &recvRequests[m]);
+        MPI_Isend(&send[m * messageSize], messageSize, MPI_QCOMP, pairRank, /*tag=*/m, mpiComm, &sendRequests[m]);
     }
 
     // wait for all exchanges to complete (MPI will automatically free the request memory)
-    MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
+    MPI_Waitall(numMessages, recvRequests.data(), MPI_STATUSES_IGNORE);
+    MPI_Waitall(numMessages, sendRequests.data(), MPI_STATUSES_IGNORE);
 
 #else
     error_commButEnvNotDistributed();
@@ -253,16 +266,17 @@ void asynchSendArray(qcomp* send, qindex numElems, int pairRank) {
 
     MPI_Comm mpiComm = comm_getMpiComm();
 
-    // we will not track nor wait for the asynch send; instead, the caller will later comm_sync()
-    MPI_Request nullReq = MPI_REQUEST_NULL;
+    // we will not track nor wait for the asynch sends; the receiver completes every
+    // matching receive before the caller's later comm_sync(), so req begins NULL
+    MPI_Request detachedReq = MPI_REQUEST_NULL;
 
     // divide the data into multiple messages
     auto [messageSize, numMessages] = dividePow2PayloadIntoMessages(numElems);
 
-    // asynchronously send the uniquely-tagged messages
+    // asynchronously send the uniquely-tagged messages, detaching requests so MPI can later release
     for (qindex m=0; m<numMessages; m++) {
-        int tag = static_cast<int>(m); // gauranteed int, but m*messageSize needs qindex
-        MPI_Isend(&send[m*messageSize], messageSize, MPI_QCOMP, pairRank, tag, mpiComm, &nullReq);
+        MPI_Isend(&send[m*messageSize], messageSize, MPI_QCOMP, pairRank, /*tag=*/m, mpiComm, &detachedReq);
+        MPI_Request_free(&detachedReq);
     }
 
 #else
@@ -283,13 +297,11 @@ void receiveArray(qcomp* dest, qindex numElems, int pairRank) {
     vector<MPI_Request> requests(numMessages, MPI_REQUEST_NULL);
 
     // listen to receive each uniquely-tagged message asynchronously (as per arxiv.org/abs/2308.07402)
-    for (qindex m=0; m<numMessages; m++) {
-        int tag = static_cast<int>(m); // gauranteed int, but m*messageSize needs qindex
-        MPI_Irecv(&dest[m*messageSize], messageSize, MPI_QCOMP, pairRank, tag, mpiComm, &requests[m]);
-    }
+    for (qindex m=0; m<numMessages; m++)
+        MPI_Irecv(&dest[m*messageSize], messageSize, MPI_QCOMP, pairRank, /*tag=*/m, mpiComm, &requests[m]);
 
     // receivers wait for all messages to be received (while sender asynch proceeds)
-    MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
+    MPI_Waitall(numMessages, requests.data(), MPI_STATUSES_IGNORE);
 
 #else
     error_commButEnvNotDistributed();
@@ -332,6 +344,7 @@ void globallyCombineNonUniformSubArrays(
     // but unlikely in other routines, their payloads can differ significantly in size,
     // so we do not know the total number of requests needed in advance
     vector<MPI_Request> requests;
+    qindex numRequests = 0;
 
     // each node broadcasts their partition (in-turn, but each is asynch)...
     for (int sendRank=0; sendRank<numNodes; sendRank++) {
@@ -340,8 +353,12 @@ void globallyCombineNonUniformSubArrays(
         // (they almost definitely send only one, but we divide for defensive design)
         auto [bigMsgSize, numBigMsgs, remMsgSize] = dividePayloadIntoMessages(numSendPerRank[sendRank]);
 
+        // ensure the accumulating request count remains representable by MPI_Waitall()
+        numRequests += numBigMsgs + (remMsgSize > 0);
+        assert_commMessageParamFitsIntoMpiInt(numRequests);
+
         // these involve big asynch messages (could be 'numSend', or the max message size)
-        for (int m=0; m<numBigMsgs; m++) {
+        for (qindex m=0; m<numBigMsgs; m++) {
             qindex recvInd = globalRecvIndPerRank[sendRank] + (m * bigMsgSize);
             requests.push_back(MPI_REQUEST_NULL);
             MPI_Ibcast(&recv[recvInd], bigMsgSize, MPI_QCOMP, sendRank, mpiComm, &requests.back());
@@ -356,7 +373,7 @@ void globallyCombineNonUniformSubArrays(
     }
 
     // wait for all broadcasts to complete (MPI will automatically free the request memory)
-    MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
+    MPI_Waitall(numRequests, requests.data(), MPI_STATUSES_IGNORE);
 
 #else
     error_commButEnvNotDistributed();
@@ -659,28 +676,36 @@ void comm_sendAmpsToRoot(int sendRank, qcomp* send, qcomp* recv, qindex numAmps)
 #if QUEST_COMPILE_MPI
 
     MPI_Comm mpiComm = comm_getMpiComm();
-
-    // only the sender and root nodes need to continue
-    int recvRank = ROOT_RANK;
     int myRank = comm_getRank();
+    int recvRank = ROOT_RANK;
+
+    // only the two sender and recv nodes need to continue
     if (myRank != sendRank && myRank != recvRank)
         return;
 
+    // avoid issuing an unmatched self-send when root already owns the payload
+    if (sendRank == recvRank) {
+
+        // ensure recvRank has its 'send' overwritten with 'recv'
+        if (myRank == recvRank && send != recv)
+            cpu_copyArray(recv, send, numAmps);
+
+        return;
+    }
+
     // create an MPI_Request for every asynch MPI call
     auto [messageSize, numMessages] = dividePow2PayloadIntoMessages(numAmps);
     vector<MPI_Request> requests(numMessages, MPI_REQUEST_NULL);
 
     // asynchronously copy 'send' in sendRank over to 'recv' in recvRank, using 
     // uniquely-tagged messages such that they may arrive out-of-order, enabling AR
-    for (qindex m=0; m<numMessages; m++) {
-        int tag = static_cast<int>(m);
+    for (qindex m=0; m<numMessages; m++)
         (myRank == sendRank)?
-            MPI_Isend(&send[m*messageSize], messageSize, MPI_QCOMP, recvRank, tag, mpiComm, &requests[m]): // sender
-            MPI_Irecv(&recv[m*messageSize], messageSize, MPI_QCOMP, sendRank, tag, mpiComm, &requests[m]); // root
-    }
+            MPI_Isend(&send[m*messageSize], messageSize, MPI_QCOMP, recvRank, /*tag=*/m, mpiComm, &requests[m]): // sender
+            MPI_Irecv(&recv[m*messageSize], messageSize, MPI_QCOMP, sendRank, /*tag=*/m, mpiComm, &requests[m]); // root
 
     // wait for all exchanges to complete (MPI will automatically free the request memory)
-    MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
+    MPI_Waitall(numMessages, requests.data(), MPI_STATUSES_IGNORE);
 
 #else
     error_commButEnvNotDistributed();
@@ -691,6 +716,9 @@ void comm_sendAmpsToRoot(int sendRank, qcomp* send, qcomp* recv, qindex numAmps)
 void comm_broadcastIntsFromRoot(int* arr, qindex length) {
 #if QUEST_COMPILE_MPI
 
+    // Bcast can only handle up to int-max; we presently only ever pass length=1
+    assert_commMessageParamFitsIntoMpiInt(length);
+
     int sendRank = ROOT_RANK;
     MPI_Bcast(arr, length, MPI_INT, sendRank, comm_getMpiComm());
 
@@ -703,6 +731,9 @@ void comm_broadcastIntsFromRoot(int* arr, qindex length) {
 void comm_broadcastUnsignedsFromRoot(unsigned* arr, qindex length) {
 #if QUEST_COMPILE_MPI
 
+    // Bcast can only handle up to int-max; we presently only ever pass length=1
+    assert_commMessageParamFitsIntoMpiInt(length);
+
     int sendRank = ROOT_RANK;
     MPI_Bcast(arr, length, MPI_UNSIGNED, sendRank, comm_getMpiComm());
 
@@ -753,7 +784,16 @@ void comm_reduceReal(qreal* localReal) {
 void comm_reduceReals(qreal* localReals, qindex numLocalReals) {
 #if QUEST_COMPILE_MPI
 
-    MPI_Allreduce(MPI_IN_PLACE, localReals, numLocalReals, MPI_QREAL, MPI_SUM, comm_getMpiComm());
+    // since numLocalReals can exceed max-int, we must divide into max-length messages
+    auto [bigMsgSize, numBigMsgs, remMsgSize] = dividePayloadIntoMessages(numLocalReals);
+    
+    // send max-length messages
+    for (qindex m=0; m<numBigMsgs; m++)
+        MPI_Allreduce(MPI_IN_PLACE, &localReals[m*bigMsgSize], bigMsgSize, MPI_QREAL, MPI_SUM, comm_getMpiComm());
+
+    // send remainder message
+    if (remMsgSize > 0)
+        MPI_Allreduce(MPI_IN_PLACE, &localReals[numBigMsgs*bigMsgSize], remMsgSize, MPI_QREAL, MPI_SUM, comm_getMpiComm());
 
 #else
     error_commButEnvNotDistributed();
@@ -777,6 +817,22 @@ bool comm_isTrueOnAllNodes(bool val) {
 }
 
 
+bool comm_isTrueOnAnyNode(bool val) {
+#if QUEST_COMPILE_MPI
+
+    // perform global OR and broadcast result back to all nodes
+    int local = (int) val;
+    int global;
+    MPI_Allreduce(&local, &global, 1, MPI_INT, MPI_LOR, comm_getMpiComm());
+    return (bool) global;
+
+#else
+    error_commButEnvNotDistributed();
+    return false;
+#endif
+}
+
+
 bool comm_isTrueOnRootNode(bool val) {
     #if QUEST_COMPILE_MPI
 
@@ -800,25 +856,25 @@ bool comm_isTrueOnRootNode(bool val) {
  */
 
 
-vector<string> comm_gatherStringsToRoot(char* localChars, int maxNumLocalChars) {
+vector<string> comm_gatherStringsToRoot(const char* localChars, int numLocalChars) {
 #if QUEST_COMPILE_MPI
 
     // no need to validate array sizes and memory alloc successes;
-    // these are trivial O(#nodes)-size arrays containing <20 chars
+    // these are trivial O(#nodes)-size arrays containing few bytes
     int numNodes = comm_getNumNodes();
 
-    // root makes one big contiguous array to receive all chars contiguously
-    vector<char> allChars(maxNumLocalChars * numNodes);
+    // root makes one big contiguous array to receive all bytes contiguously
+    vector<char> allChars(numLocalChars * numNodes);
 
-    // all nodes send root all their local chars
+    // all nodes send root their fixed-width byte sequence
     int recvRank = ROOT_RANK;
-    MPI_Gather(localChars, maxNumLocalChars, MPI_CHAR, allChars.data(),
-        maxNumLocalChars, MPI_CHAR, recvRank, comm_getMpiComm());
+    MPI_Gather(localChars, numLocalChars, MPI_CHAR, allChars.data(),
+        numLocalChars, MPI_CHAR, recvRank, comm_getMpiComm());
 
-    // divide allChars into stings, delimited by each node's terminal char
+    // preserve every byte, including embedded nulls, in fixed-width strings
     vector<string> out(numNodes);
     for (int r=0; r<numNodes; r++)
-        out[r] = std::string(&allChars[r*maxNumLocalChars]);
+        out[r] = std::string(&allChars[r*numLocalChars], numLocalChars);
 
     return out;
 
diff --git a/quest/src/comm/comm_routines.hpp b/quest/src/comm/comm_routines.hpp
index e75e889f6..8e8aa45ba 100644
--- a/quest/src/comm/comm_routines.hpp
+++ b/quest/src/comm/comm_routines.hpp
@@ -69,6 +69,8 @@ void comm_reduceReals(qreal* localReals, qindex numLocalReals);
 
 bool comm_isTrueOnAllNodes(bool val);
 
+bool comm_isTrueOnAnyNode(bool val);
+
 bool comm_isTrueOnRootNode(bool val);
 
 
@@ -77,7 +79,7 @@ bool comm_isTrueOnRootNode(bool val);
  * GATHER METHODS
  */
 
-vector<std::string> comm_gatherStringsToRoot(char* localChars, int maxNumLocalChars);
+vector<std::string> comm_gatherStringsToRoot(const char* localChars, int numLocalChars);
 
 
 
diff --git a/quest/src/core/bitwise.hpp b/quest/src/core/bitwise.hpp
index 1c079bdfe..8501b54b5 100644
--- a/quest/src/core/bitwise.hpp
+++ b/quest/src/core/bitwise.hpp
@@ -200,7 +200,7 @@ INLINE qindex getValueOfBits(qindex number, const int* bitIndices, int numIndice
     qindex value = 0;
 
     for (int i=0; i<numIndices; i++)
-        value |= getBit(number, bitIndices[i]) << i;
+        value |= ((qindex) getBit(number, bitIndices[i])) << i;
 
     return value;
 }
@@ -446,7 +446,7 @@ INLINE qindex getIntegerFromBits(const int* bits, int numBits) {
     qindex value = 0;
 
     for (int i=0; i<numBits; i++)
-        value |= bits[i] << i;
+        value |= ((qindex) bits[i]) << i;
 
     return value;
 }
diff --git a/quest/src/core/errors.cpp b/quest/src/core/errors.cpp
index c78f1c727..d5a91ab29 100644
--- a/quest/src/core/errors.cpp
+++ b/quest/src/core/errors.cpp
@@ -20,6 +20,7 @@
 
 #include <iostream>
 #include <string>
+#include <limits>
 
 using std::string;
 
@@ -236,6 +237,12 @@ void assert_commPayloadIsPowerOf2(qindex numAmps) {
         raiseInternalError("A communication function was given a payload which was unexpectedly not a power of 2, as implied by preconditions.");
 }
 
+void assert_commMessageParamFitsIntoMpiInt(qindex value) {
+
+    if (value > std::numeric_limits<int>::max())
+        raiseInternalError("A function attempted to communicate with a message parameter which exceeded the maximum int accepted by MPI.");
+}
+
 void assert_commQuregIsDistributed(Qureg qureg) {
 
     if (!qureg.isDistributed)
@@ -419,7 +426,7 @@ void assert_exponentMatchesTemplateParam(qcomp exponent, bool hasPower, bool use
     assert_exponentMatchesTemplateParam(exponent, hasPower);
 
     if (!hasPower && useRealPow)
-        raiseInternalError("A CPU or GPU subroutine received an invalid combination of HasPower=false and UseRealPow=false template parameters");
+        raiseInternalError("A CPU or GPU subroutine received an invalid combination of HasPower=false and UseRealPow=true template parameters");
 }
 
 
@@ -756,6 +763,11 @@ void error_cuQuantumTempCpuAllocFailed() {
     raiseInternalError("Attempted allocation of temporary host-memory for a cuQuantum routine failed.");
 }
 
+void error_cuQuantumIncompatibleWithMaximallyDistributedDephasing() {
+
+    raiseInternalError("The current extreme distribution scheme (a single density matrix column per node) is incompatible with our optimised implementation of one-qubit dephasing using cuQuantum.");
+}
+
 
 
 /*
diff --git a/quest/src/core/errors.hpp b/quest/src/core/errors.hpp
index 7db4e69f5..753a6d35b 100644
--- a/quest/src/core/errors.hpp
+++ b/quest/src/core/errors.hpp
@@ -21,6 +21,7 @@
 #include "quest/include/qureg.h"
 #include "quest/include/matrices.h"
 
+#include <cstddef>
 #include <string>
 
 using std::string;
@@ -111,6 +112,8 @@ void assert_commBoundsAreValid(Qureg qureg, qindex sendInd, qindex recvInd, qind
 
 void assert_commPayloadIsPowerOf2(qindex numAmps);
 
+void assert_commMessageParamFitsIntoMpiInt(qindex value);
+
 void assert_commQuregIsDistributed(Qureg qureg);
 
 void assert_commFullStateDiagMatrIsDistributed(FullStateDiagMatr matr);
@@ -307,6 +310,8 @@ void error_cuQuantumInitOrFinalizedButNotCompiled();
 
 void error_cuQuantumTempCpuAllocFailed();
 
+void error_cuQuantumIncompatibleWithMaximallyDistributedDephasing();
+
 
 
 /*
diff --git a/quest/src/core/fastmath.hpp b/quest/src/core/fastmath.hpp
index 79884659f..da4ef8d57 100644
--- a/quest/src/core/fastmath.hpp
+++ b/quest/src/core/fastmath.hpp
@@ -96,7 +96,7 @@ INLINE void fast_getSubQuregValues(qindex basisStateIndex, int* numQubitsPerSubQ
     if (areSigned)
         for (int v=0; v<numSubQuregs; v++)
             if (getBit(outValues[v], numQubitsPerSubQureg[v]-1))
-                outValues[v] -= powerOf2(numQubitsPerSubQureg[v] - 1);
+                outValues[v] -= powerOf2(numQubitsPerSubQureg[v]);
 }
 
 
diff --git a/quest/src/core/localiser.cpp b/quest/src/core/localiser.cpp
index 83a23b921..8a89fea38 100644
--- a/quest/src/core/localiser.cpp
+++ b/quest/src/core/localiser.cpp
@@ -2079,7 +2079,7 @@ qcomp localiser_statevec_calcExpecPauliStrSum(Qureg qureg, PauliStrSum sum) {
     std::unordered_map<Key, vector<Term>> groups;
 
     // group sum's terms into those with identical communication patterns
-    for (int i=0; i<sum.numTerms; i++) {
+    for (qindex i=0; i<sum.numTerms; i++) {
         Term term = tuple{sum.strings[i], sum.coeffs[i]};
         Key totalKey = paulis_getKeyOfSameMixedAmpsGroup(sum.strings[i]);
         Key prefixKey = getBitsLeftOfIndex(totalKey, qureg.logNumAmpsPerNode - 1); // 0 if !qureg.isDistributed
diff --git a/quest/src/core/memory.cpp b/quest/src/core/memory.cpp
index 7f11494a1..9ff65b05a 100644
--- a/quest/src/core/memory.cpp
+++ b/quest/src/core/memory.cpp
@@ -48,17 +48,23 @@
     if (hasBuffer && numNodes > 1)
         maxLocalNumAmps /= 2; // floors
 
+    // avoid log2(0) when even one amplitude cannot fit
+    if (maxLocalNumAmps == 0)
+        return -1;
+
+    // distribution increases the global amplitude exponent before converting it to logical
+    // qubits; this order matters for density matrices, whose amplitude exponent is doubled
+    int maxGlobalAmpQubits = std::floor(std::log2(maxLocalNumAmps)) + logBase2(numNodes);
+
     // density matrices require square more memory, so halve (flooring) the number of qubits
-    int maxLocalNumQubits = std::floor(std::log2(maxLocalNumAmps));
+    int maxGlobalNumQubits = maxGlobalAmpQubits;
     if (isDensMatr)
-        maxLocalNumQubits /= 2; // floors
+        maxGlobalNumQubits /= 2; // floors
 
     // superoperators require square more memory still, so halve (flooring) the number of qubits
     if (isSuperOp)
-        maxLocalNumQubits /= 2; // floors
+        maxGlobalNumQubits /= 2; // floors
 
-    // doubling nodes permits 1 additional qubit
-    int maxGlobalNumQubits = maxLocalNumQubits + logBase2(numNodes);
     return maxGlobalNumQubits;
 }
 
@@ -120,7 +126,7 @@ int mem_getMaxNumSuperOpQubitsBeforeIndexOverflow() {
 
 qindex mem_getMaxNumKrausMapMatricesBeforeIndexOverflow(int numQubits) {
 
-    qindex numElemPerMatrix = powerOf2(2 * numQubits);
+    qindex numElemPerMatrix = powerOf2(2 * numQubits); // caller prior validates no overflow
     qindex maxNumTotalElems = std::numeric_limits<qindex>::max();
     qindex maxNumMatrices = maxNumTotalElems / numElemPerMatrix; // floors
     return maxNumMatrices;
@@ -196,9 +202,10 @@ int mem_getMaxNumSuperOpQubitsBeforeGlobalMemSizeofOverflow() {
 
 qindex mem_getMaxNumKrausMapMatricesBeforeLocalMemSizeofOverflow(int numQubits) {
 
-    qindex numMatrWithMaxTotalElems = mem_getMaxNumKrausMapMatricesBeforeIndexOverflow(numQubits);
-    qindex numMatrWithMaxTotalMem = numMatrWithMaxTotalElems / sizeof(qcomp); // floors
-    return numMatrWithMaxTotalMem;
+    qindex numElemsPerMatrix = powerOf2(2 * numQubits); // caller prior validates no overflow
+    size_t maxNumElems = std::numeric_limits<size_t>::max() / sizeof(qcomp); // floors
+    size_t maxNumMatrices = maxNumElems / static_cast<size_t>(numElemsPerMatrix); // floors
+    return static_cast<qindex>(maxNumMatrices);
 }
 
 
@@ -368,6 +375,11 @@ bool mem_canMatrixFitInMemory(int numQubits, bool isDense, int numNodes, qindex
 
     // work out the maximum "local" qubits that can fit in memory
     qindex maxLocalNumElems = memBytesPerNode / sizeof(qcomp); // floors
+
+    // max distribution (useful for tests) has 1 amp per node, = 0 local qubits!
+    if (maxLocalNumElems == 0)
+        return false;
+
     int maxLocalNumQubits  = std::floor(std::log2(maxLocalNumElems));
 
     // dense matrices (as opposed to diagonals) require square more memory
diff --git a/quest/src/core/parser.cpp b/quest/src/core/parser.cpp
index 9d9194a3f..86ca4f611 100644
--- a/quest/src/core/parser.cpp
+++ b/quest/src/core/parser.cpp
@@ -118,7 +118,7 @@ namespace regexes {
 bool isWhiteSpace(char ch) {
 
     // matches patterns:space and more (additionally; newlines)
-    return isspace(ch); 
+    return std::isspace(static_cast<unsigned char>(ch));
 }
 
 bool isNotWhiteSpace(char ch) {
diff --git a/quest/src/core/paulilogic.cpp b/quest/src/core/paulilogic.cpp
index 58ddc39f8..a56b61d50 100644
--- a/quest/src/core/paulilogic.cpp
+++ b/quest/src/core/paulilogic.cpp
@@ -244,7 +244,7 @@ PAULI_MASK_TYPE paulis_getKeyOfSameMixedAmpsGroup(PauliStr str) {
     for (int i=0; i<=maxInd; i++) {
         int pauli = paulis_getPauliAt(str, i);
         int isXY = (pauli == 1 || pauli == 2);
-        key |= (isXY << i);
+        key |= ((PAULI_MASK_TYPE) isXY) << i;
     }
 
     return key;
@@ -304,7 +304,7 @@ qindex paulis_getTargetBitMask(PauliStrSum sum) {
     qindex mask = 0;
 
     // mask has 1 where any str has a != Id
-    for (int t=0; t<sum.numTerms; t++)
+    for (qindex t=0; t<sum.numTerms; t++)
         mask |= paulis_getTargetBitMask(sum.strings[t]);
 
     return mask;
@@ -379,7 +379,7 @@ qindex paulis_getNumTermsInPauliStrSumProdOfAdjointWithSelf(PauliStrSum in) {
     // adj(in).in has fewer terms than the numTerms^2 bound, since 
     // a.a = I (causing -n and +1 below) and a.b ~ b.a (causing /2);
     // we do not however consider any cancellations of coefficients
-    int n = in.numTerms;
+    qindex n = in.numTerms;
     return 1 + (n*n - n)/2;
 }
 
diff --git a/quest/src/core/printer.cpp b/quest/src/core/printer.cpp
index 863317c8f..a7207dcb2 100644
--- a/quest/src/core/printer.cpp
+++ b/quest/src/core/printer.cpp
@@ -473,18 +473,22 @@ void print_oneFewerNewlines() {
 string printer_getMemoryWithUnitStr(size_t numBytes) {
 
     // retain "bytes" instead of "B" since more recognisable
-    vector<string> units = {"bytes", "KiB", "MiB", "GiB", "TiB", "PiB", "EiB", "ZiB", "YiB"};
+    vector<string> units = {"bytes", "KiB", "MiB", "GiB", "TiB", "PiB", "EiB"};
 
-    // unit sizes (in bytes) may overflow; but then so too will numBytes
+    // unit sizes in bytes
     vector<size_t> sizes(units.size());
     for (size_t i=0; i<units.size(); i++)
         sizes[i] = powerOf2(10 * i); // = 1024^i
 
-    // find the biggest unit for which numBytes/sizes[ind] > 1
-    int ind = 1;
-    while (numBytes > sizes[ind])
+    // find the biggest unit no larger than numBytes
+    size_t ind = 0;
+    while (ind + 1 < sizes.size() && numBytes >= sizes[ind + 1])
         ind++;
-    ind--;
+
+    // at four significant figures, values at least halfway from 1023 to 1024
+    // round to 1024 of the current unit, which is more cleanly shown as 1 of the next
+    if (ind + 1 < sizes.size() && numBytes >= sizes[ind + 1] - sizes[ind] / 2)
+        return "1 " + units[ind + 1];
 
     // express numBytes in terms of new unit, forcefully rounding to 4 sig-figs max
     qreal frac = numBytes / static_cast<qreal>(sizes[ind]);
@@ -974,7 +978,7 @@ void populateManyColumnQcompmatr(qcompmatr &matr, qindex startRow, qindex startC
     // note this achieves node consensus; every node populates their matr, which is much more
     // demanding than only the root note, however we expect this is not a big deal; live printing
     // in distributed settings is exceptionally rare and will print small, tractable subsets
-    localiser_densmatr_getAmps(ptrs.data(), qureg, startRow, startCol, matr.size(), matr[0].capacity());
+    localiser_densmatr_getAmps(ptrs.data(), qureg, startRow, startCol, matr.size(), matr[0].size());
 }
 
 
diff --git a/quest/src/core/randomiser.cpp b/quest/src/core/randomiser.cpp
index 7b35a29fc..c3c47751f 100644
--- a/quest/src/core/randomiser.cpp
+++ b/quest/src/core/randomiser.cpp
@@ -135,7 +135,7 @@ int rand_getRandomSingleQubitOutcome(qreal probOfZero) {
     qreal sample = distrib(mainGenerator);
 
     // produces 0 with probOfZero probability
-    return sample > probOfZero;
+    return sample >= probOfZero;
 }
 
 
diff --git a/quest/src/core/utilities.cpp b/quest/src/core/utilities.cpp
index 71facdf64..14c4bb1b7 100644
--- a/quest/src/core/utilities.cpp
+++ b/quest/src/core/utilities.cpp
@@ -398,6 +398,11 @@ bool util_willProdOverflow(vector<qindex> terms) {
     qindex max = std::numeric_limits<qindex>::max();
     qindex prod = 1;
 
+    // any zero factor makes the full product zero and cannot overflow
+    for (auto x : terms)
+        if (x == 0)
+            return false;
+
     for (auto x : terms) {
 
         // division floors, so prod strictly exceed
@@ -417,7 +422,7 @@ bool util_willSumOverflow(vector<qindex> terms) {
 
     for (auto x : terms) {
 
-        if (sum >= max - x)
+        if (sum > max - x)
             return true;
         
         sum += x;
@@ -635,9 +640,9 @@ bool getHermiticity(T elems, qindex dim, qreal eps) {
     /// consider multithreading or GPU-accelerating this
     /// when caller is big and e.g. has GPU memory
 
-    // check adjoint(elems) == elems
+    // check adjoint(elems) == elems, including that each diagonal is real
     for (qindex r=0; r<dim; r++)
-        for (qindex c=0; c<r; c++)
+        for (qindex c=0; c<=r; c++)
             if (!isApprox(elems[r][c], std::conj(elems[c][r]), eps))
                 return false;
 
@@ -857,6 +862,7 @@ bool util_isCPTP(KrausMap map, qreal eps) {
     /// use a multithreaded routine
 
     *(map.isApproxCPTP) = 1;
+    qreal epsSquared = eps * eps;
 
     // check whether each element satisfies Identity = sum dagger(m)*m
     for (qindex r=0; r<map.numRows; r++) {
@@ -870,7 +876,7 @@ bool util_isCPTP(KrausMap map, qreal eps) {
 
             // fail if too distant from Identity element...
             qreal distSquared = std::norm(elem - (r==c));
-            if (distSquared > eps) {
+            if (distSquared > epsSquared) {
 
                 // by recording the result and returning immediately
                 *(map.isApproxCPTP) = 0;
diff --git a/quest/src/core/validation.cpp b/quest/src/core/validation.cpp
index beeae12f1..d21988dc3 100644
--- a/quest/src/core/validation.cpp
+++ b/quest/src/core/validation.cpp
@@ -1642,7 +1642,7 @@ void validate_randomSeeds(unsigned* seeds, int numSeeds, const char* caller) {
     }
 
     assertThat(!isNull, report::RANDOM_SEEDS_PTR_IS_NULL, caller);
-    assertThat(numRootSeeds > 0, report::INVALID_NUM_RANDOM_SEEDS, {{"${NUM_SEEDS}", numSeeds}}, caller);
+    assertThat(numRootSeeds > 0, report::INVALID_NUM_RANDOM_SEEDS, {{"${NUM_SEEDS}", numRootSeeds}}, caller);
 }
 
 void validate_newEpsilonValue(qreal eps, const char* caller) {
@@ -2585,21 +2585,19 @@ void validate_matrixRowsAllSameSize(vector<vector<qcomp>> matrix, const char* ca
 
     size_t dim = matrix[0].size();
 
-    size_t row=0;
-    bool allSame = true;
-    for (row=0; row<matrix.size() && allSame; row++)
-        allSame = (matrix[row].size() == dim);
+    size_t outlierInd = 0;
+    while (outlierInd < matrix.size() && matrix[outlierInd].size() == dim)
+        outlierInd++;
 
-    // extremely lazily avoiding seg-fault from premature matrix[row==end])
-    if (allSame)
+    if (outlierInd == matrix.size())
         return;
 
     tokenSubs vars = {
-        {"{FIRST_LEN}",    dim},
-        {"${OUTLIER_IND}", row},
-        {"${OUTLIER_LEN}", matrix[row].size()}};
+        {"${FIRST_LEN}",    dim},
+        {"${OUTLIER_IND}", outlierInd},
+        {"${OUTLIER_LEN}", matrix[outlierInd].size()}};
 
-    assertThat(allSame, report::NESTED_VECTOR_MATRIX_HAS_INCONSISTENT_NUM_COLUMNS, vars, caller);
+    assertThat(false, report::NESTED_VECTOR_MATRIX_HAS_INCONSISTENT_NUM_COLUMNS, vars, caller);
 }
 
 
@@ -3406,7 +3404,7 @@ void validate_superOpDimMatchesTargs(SuperOp op, int numTargets, const char* cal
     if (!global_isValidationEnabled)
         return;
 
-    tokenSubs vars = {{"${OP_QUBITS}", op.numQubits}, {"{NUM_TARGS}", numTargets}};
+    tokenSubs vars = {{"${OP_QUBITS}", op.numQubits}, {"${NUM_TARGS}", numTargets}};
     assertThat(op.numQubits == numTargets, report::SUPER_OP_SIZE_MISMATCHES_NUM_TARGETS, vars, caller);
 }
 
@@ -3815,7 +3813,7 @@ void validate_newPauliStrSumParams(qindex numTerms, const char* caller) {
     size_t memPerTerm = sizeof(qcomp) + sizeof(PauliStr);
     size_t maxNumTerms = std::numeric_limits<size_t>::max() / memPerTerm;
     tokenSubs vars = {{"${NUM_TERMS}", numTerms}, {"${NUM_BYTES_PER_TERM}", memPerTerm}, {"${MAX_NUM_TERMS}", maxNumTerms}};
-    assertThat(numTerms < (qindex) maxNumTerms, report::NEW_PAULI_STR_SUM_MEM_WOULD_OVERFLOW, vars, caller);
+    assertThat(numTerms <= (qindex) maxNumTerms, report::NEW_PAULI_STR_SUM_MEM_WOULD_OVERFLOW, vars, caller);
 
     // attempt to fetch RAM, and simply return if we fail; if we unknowingly
     // didn't have enough RAM, then alloc validation will trigger later
@@ -4407,7 +4405,7 @@ void validate_measurementOutcomesFitInGpuMem(Qureg qureg, int numQubits, const c
     size_t memNeeded = sizeof(qreal) * numOutcomes;
 
     tokenSubs vars = {
-        {"{NUM_QUBITS}",    numQubits}, 
+        {"${NUM_QUBITS}",   numQubits},
         {"${NUM_OUTCOMES}", numOutcomes}, 
         {"${MEM_NEEDED}",   memNeeded}, 
         {"${MEM_AVAIL}",    memAvail}
diff --git a/quest/src/cpu/cpu_config.cpp b/quest/src/cpu/cpu_config.cpp
index bd51236bb..34b234b89 100644
--- a/quest/src/cpu/cpu_config.cpp
+++ b/quest/src/cpu/cpu_config.cpp
@@ -136,7 +136,7 @@ int cpu_getCurrentNumThreads() {
 qindex getNumPagesToContainArray(long pageLen, qindex arrLen) {
 
     // round up to the nearest page
-    return static_cast<qindex>(std::ceil(arrLen / (qreal) pageLen));
+    return arrLen / pageLen + (arrLen % pageLen != 0);
 }
 
 
@@ -223,6 +223,11 @@ qcomp* cpu_allocNumaArray(qindex length) {
     // use integer type for safe address arithmetic below
     uintptr_t offsetAddr = reinterpret_cast<uintptr_t>(rawAddr);
 
+    // mbind() accepts a bit array, which must span every configured NUMA node
+    size_t bitsPerMaskWord = 8 * sizeof(unsigned long);
+    size_t numMaskWords = (numNodes + bitsPerMaskWord - 1) / bitsPerMaskWord;
+    vector<unsigned long> nodeMask(numMaskWords, 0);
+
     for (int node=0, shift=numNodes; node < numNodes; ++node) {
 
         // decide number of pages to bind to NUMA node
@@ -230,16 +235,24 @@ qcomp* cpu_allocNumaArray(qindex length) {
         qindex numPagesInNode = baseNumPagesPerNode + (shift <= 0);
         qindex numBytesInNode = numPagesInNode * pageSize; // validation prevents overflow
 
-        // bind those pages from the offset address to the node (identified by mask)
-        unsigned long nodeMask = 1UL << node;
-        unsigned long numBitsInMask = 8 * sizeof(nodeMask);
-        void* nodeAddr = reinterpret_cast<void*>(offsetAddr);
-        long success = mbind(nodeAddr, numBytesInNode, MPOL_BIND, &nodeMask, numBitsInMask, 0);
+        // skip nodes assigned no pages (which occurs when there are fewer pages than nodes)
+        if (numPagesInNode > 0) {
+            
+            // bind those pages from the offset address to the node (identified by mask)
+            size_t maskWord = node / bitsPerMaskWord;
+            unsigned long maskBit = 1UL << (node % bitsPerMaskWord);
+            nodeMask[maskWord] = maskBit;
+
+            void* nodeAddr = reinterpret_cast<void*>(offsetAddr);
+            long success = mbind(nodeAddr, numBytesInNode, MPOL_BIND, nodeMask.data(), numMaskWords * bitsPerMaskWord, 0);
+
+            nodeMask[maskWord] = 0;
 
-        // treat bind failure as internal error (even though it can result from insufficient kernel mem),
-        // rather than permitting silent fallback to non-NUMA awareness which might be astonishingly slow
-        if (success == -1)
-            error_numaBindingFailed();
+            // treat bind failure as internal error (even though it can result from insufficient kernel mem),
+            // rather than permitting silent fallback to non-NUMA awareness which might be astonishingly slow
+            if (success == -1)
+                error_numaBindingFailed();
+        }
 
         // prepare next node's address
         offsetAddr += numPagesInNode * pageSize;
@@ -451,7 +464,7 @@ void cpu_copyMatrix(qcomp** dest, vector<vector<qcomp>> src, qindex dim) {
 void cpu_copyPauliStrSum(PauliStrSum out, PauliStr* strings, qcomp* coeffs) {
 
     // serially copy data over to new heap memory
-    for (int i=0; i<out.numTerms; i++) {
+    for (qindex i=0; i<out.numTerms; i++) {
         out.strings[i] = strings[i];
         out.coeffs[i] = coeffs[i];
     }
diff --git a/quest/src/cpu/cpu_subroutines.cpp b/quest/src/cpu/cpu_subroutines.cpp
index c1ba9b3aa..6968d403e 100644
--- a/quest/src/cpu/cpu_subroutines.cpp
+++ b/quest/src/cpu/cpu_subroutines.cpp
@@ -829,7 +829,7 @@ void cpu_statevec_allTargDiagMatr_sub(Qureg qureg, FullStateDiagMatr matr, qcomp
     // every iteration modifies one amp, using one element
     qindex numIts = qureg.numAmpsPerNode;
 
-    #pragma omp parallel for if(qureg.isMultithreaded||qureg.isMultithreaded)
+    #pragma omp parallel for if(qureg.isMultithreaded||matr.isMultithreaded)
     for (qindex n=0; n<numIts; n++) {
 
         cpu_qcomp elem = elems[n];
@@ -1169,6 +1169,9 @@ void cpu_statevec_setQuregToWeightedSum_sub(Qureg outQureg, vector<qcomp> coeffs
     // use template param to compile-time unroll inner loop below
     int numQuregs = accel_tryUseCompileTimeValue<NumQuregs>(inQuregs.size());
 
+    // below, we choose to multithread only if outQureg is multithreaded,
+    // regardless of whether any of the inQuregs are multithreaded
+
     #pragma omp parallel for if(outQureg.isMultithreaded)
     for (qindex n=0; n<numIts; n++) {
 
@@ -1192,7 +1195,7 @@ void cpu_densmatr_mixQureg_subA(qreal outProb, Qureg outQureg, qreal inProb, Qur
     
     qindex numIts = outQureg.numAmpsPerNode;
 
-    #pragma omp parallel for if(outQureg.isMultithreaded)
+    #pragma omp parallel for if(outQureg.isMultithreaded||inDensMatr.isMultithreaded)
     for (qindex n=0; n<numIts; n++)
         outAmps[n] = (outProb * outAmps[n]) + (inProb * inAmps[n]);
 }
@@ -1207,6 +1210,9 @@ void cpu_densmatr_mixQureg_subB(qreal outProb, Qureg outQureg, qreal inProb, Qur
     qindex numIts = outQureg.numAmpsPerNode;
     qindex dim = inStateVec.numAmps;
 
+    // only multithread if the larger density matrix outQureg is multithreaded;
+    // disregard whether inStateVec is multithreaded 
+
     #pragma omp parallel for if(outQureg.isMultithreaded)
     for (qindex n=0; n<numIts; n++) {
 
@@ -2102,7 +2108,7 @@ void cpu_statevec_calcProbsOfAllMultiQubitOutcomes_sub(qreal* outProbs, Qureg qu
 
     // clear amps (may be compile-time unrolled, or parallelised)
     #pragma omp parallel for if(parallelise)
-    for (int i=0; i<numOutcomes; i++)
+    for (qindex i=0; i<numOutcomes; i++)
         outProbs[i] = 0;
 
     #pragma omp parallel for if(qureg.isMultithreaded)
@@ -2149,7 +2155,7 @@ void cpu_densmatr_calcProbsOfAllMultiQubitOutcomes_sub(qreal* outProbs, Qureg qu
     
     // clear amps; be compile-time unrolled, and/or parallelised (independent of qureg)
     #pragma omp parallel for if(parallelise)
-    for (int i=0; i<numOutcomes; i++)
+    for (qindex i=0; i<numOutcomes; i++)
         outProbs[i] = 0;
 
     #pragma omp parallel for if(qureg.isMultithreaded)
@@ -2219,7 +2225,7 @@ qreal cpu_densmatr_calcHilbertSchmidtDistance_sub(Qureg quregA, Qureg quregB) {
     // every local amp contributes to the reduction
     qindex numIts = quregA.numAmpsPerNode;
 
-    #pragma omp parallel for reduction(+:dist) if(quregA.isMultithreaded||quregA.isMultithreaded)
+    #pragma omp parallel for reduction(+:dist) if(quregA.isMultithreaded||quregB.isMultithreaded)
     for (qindex n=0; n<numIts; n++)
         dist += norm(ampsA[n] - ampsB[n]); // |A-B|^2
 
@@ -2241,6 +2247,9 @@ qcomp cpu_densmatr_calcFidelityWithPureState_sub(Qureg rho, Qureg psi) {
     // every local density matrix amp contributes to the reduction
     qindex numIts = rho.numAmpsPerNode;
 
+    // below, we parallelise only if rho is multithreaded; we disregard
+    // whether the smaller statevector psi is multithreaded
+
     #pragma omp parallel for reduction(+:fidRe,fidIm) if(rho.isMultithreaded)
     for (qindex n=0; n<numIts; n++) {
 
diff --git a/quest/src/gpu/gpu_config.cpp b/quest/src/gpu/gpu_config.cpp
index 001cc62c0..5feff91a5 100644
--- a/quest/src/gpu/gpu_config.cpp
+++ b/quest/src/gpu/gpu_config.cpp
@@ -370,14 +370,11 @@ int gpu_getMaxNumThreadsPerBlock() {
 }
 
 
-std::array<char,17> getBoundGpuUuid() {
+std::array<char,16> getBoundGpuUuid() {
 #if QUEST_COMPILE_CUDA
     assert_gpuHasBeenBound(hasGpuBeenBound);
 
-    constexpr int numUuidChars = 16;
-    constexpr int numOutChars = numUuidChars + 1; // terminal char
-
-    std::array<char,numOutChars> out;
+    constexpr int numUuidBytes = 16;
     cudaUUID_t uuid;
 
     // ROCm v5's cudaDeviceProp doesn't have a uuid field
@@ -391,12 +388,8 @@ std::array<char,17> getBoundGpuUuid() {
         uuid = prop.uuid;
     #endif
 
-    // copy char[16] to out[17] (not human readable)
-    std::copy_n(uuid.bytes, numUuidChars, out.begin());
-
-    // include terminal char so that subsequent to-string
-    // operations will succeed without knowing string length
-    out[numOutChars-1] = '\0';
+    std::array<char,numUuidBytes> out;
+    std::copy_n(uuid.bytes, numUuidBytes, out.begin());
     return out;
 
 #else
@@ -437,11 +430,11 @@ bool gpu_areAnyNodesBoundToSameGpu() {
     if (!comm_isActive())
         return false;
 
-    // obtain bound GPU's UUID; a unique identifier 16-char identifier
-    auto uuidStr = getBoundGpuUuid();
+    // obtain the bound GPU's fixed-width binary UUID
+    auto uuid = getBoundGpuUuid();
 
-    // we can repurpose string-to-root sending to collect all uuids
-    auto allUuids = comm_gatherStringsToRoot(uuidStr.data(), uuidStr.size());
+    // gather every byte, including embedded nulls, and compare the identifiers
+    auto allUuids = comm_gatherStringsToRoot(uuid.data(), static_cast<int>(uuid.size()));
     auto uniqueUuids = std::set<std::string>(allUuids.begin(), allUuids.end());
 
     // and assess whether they're all unique (non-root's bools are overwritten)
diff --git a/quest/src/gpu/gpu_cuquantum.cuh b/quest/src/gpu/gpu_cuquantum.cuh
index 6323f549f..44e43a97f 100644
--- a/quest/src/gpu/gpu_cuquantum.cuh
+++ b/quest/src/gpu/gpu_cuquantum.cuh
@@ -284,6 +284,11 @@ void cuquantum_densmatr_oneQubitDephasing_subB(Qureg qureg, int ketQubit, qreal
     // note it must still be a suffix qubit since cuQuantum does not know qureg is distributed
     int targ = qureg.logNumAmpsPerNode - 1; // leftmost suffix bra qubit
 
+    // only an extreme distribution (one density-matrix column per node) makes
+    // the chosen target equal the control, which cuStateVec forbids
+    if (targ == ketQubit)
+        error_cuQuantumIncompatibleWithMaximallyDistributedDephasing();
+
     bool conj = false;
     auto ctrls  = lists_getList64({ketQubit});
     auto states = lists_getList64({!braBit});
diff --git a/quest/src/gpu/gpu_kernels.cuh b/quest/src/gpu/gpu_kernels.cuh
index 42617c4a6..c75ae2801 100644
--- a/quest/src/gpu/gpu_kernels.cuh
+++ b/quest/src/gpu/gpu_kernels.cuh
@@ -61,7 +61,7 @@
 
 
 __forceinline__ __device__ qindex getThreadInd() {
-    return blockIdx.x*blockDim.x + threadIdx.x;
+    return static_cast<qindex>(blockIdx.x) * blockDim.x + threadIdx.x;
 }
 
 
@@ -71,8 +71,8 @@ __host__ qindex getNumBlocks(qindex numThreads, int numThreadsPerBlock) {
     /// improve this with cudaOccupancyMaxPotentialBlockSize(),
     /// making it function specific
 
-    // CUDA ceil
-    return ceil(numThreads / static_cast<qreal>(numThreadsPerBlock));
+    // exact integer ceiling, without losing precision for large launches
+    return numThreads / numThreadsPerBlock + (numThreads % numThreadsPerBlock != 0);
 }
 
 
diff --git a/quest/src/gpu/gpu_subroutines.cpp b/quest/src/gpu/gpu_subroutines.cpp
index e9936c9c3..0cb41c38d 100644
--- a/quest/src/gpu/gpu_subroutines.cpp
+++ b/quest/src/gpu/gpu_subroutines.cpp
@@ -498,9 +498,10 @@ void gpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co
         qindex numThreads = gpu_getMaxNumConcurrentThreads();
         int numThreadsPerBlock = gpu_getNumThreadsPerBlock();
         
-        // use strictly 2^# threads to maintain precondition of all kernels
+        // use at most a power-of-two number of concurrent threads to maintain
+        // the kernels' divisibility precondition without exceeding the hardware bound
         if (!isPowerOf2(numThreads))
-            numThreads = util_getNextPowerOf2(numThreads);
+            numThreads = util_getNextPowerOf2(numThreads) / 2;
 
         // no point in dispatching more threads than batches
         if (numThreads > numBatches)
@@ -1506,7 +1507,11 @@ INSTANTIATE_FUNC_OPTIMISED_FOR_NUM_TARGS( void, gpu_densmatr_partialTrace_sub, (
 qreal gpu_statevec_calcTotalProb_sub(Qureg qureg) {
 
 #if QUEST_COMPILE_CUQUANTUM
-    return cuquantum_statevec_calcTotalProb_sub(qureg);
+    // cuStateVec's Z-basis routine requires a target qubit, which does not exist
+    // when a maximally distributed statevector stores one amplitude per node
+    return (qureg.logNumAmpsPerNode > 0)?
+        cuquantum_statevec_calcTotalProb_sub(qureg):
+        thrust_statevec_calcTotalProb_sub(qureg);
 
 #elif QUEST_COMPILE_CUDA
     return thrust_statevec_calcTotalProb_sub(qureg);