@Magma Please put the original pseudocode back in the Question! Even though it had some small errors, it was far easier to understand than the C++ code. It's far better to note the mistakes in information rather than removing it entirely. I always tell my math students, "Cross out; don't erase!"

@Carl Love By just updating the columns of H that will be changed due to the changes in R, as I mentioned earlier, I've got the time down to 3 seconds for a 64x64 matrix.

Paar:= proc(M::Matrix)
option `Author: Carl Love <carl.j.love@gmail.com> 25-Oct-2019`;
local 
     R:= rtable(M, datatype= integer[8]), n:= upperbound(M)[2], xors:= add(R)-n, 
     H:= rtable(R^%T.R, datatype= integer[8]), 
     hmax, maxij, mi, mj, lastcol, Update, H_as
;
     for lastcol from 1+n do
          H_as:= rtable(antisymmetric, H);
          hmax:= max(H_as);
          if hmax <= 1 then return (R, xors) fi;
          member(hmax, H_as, 'maxij'); 
          (mi, mj):= maxij;
          R(.., lastcol):= R[.., mi] *~ R[.., mj]; #Add new column to R.
          #Update old columns of R:
          R[.., [maxij]]:= R[.., [maxij]] *~  (1 -~ R[.., [-1$2]]); 
          xors:= xors - (hmax-1);
          H(.., lastcol):= <seq(0, 2..lastcol)>; #Add column to H.
          H(lastcol, ..):= `<|>`(seq(0, 1..lastcol)); #Add row to H.
          Update:= R^%T.R[..,[maxij, -1]]; #3 new columns/rows for H.
          H[.., [maxij, -1]]:= Update; #Update H's columns. 
          H[[maxij,-1], ..]:= Update^%T #Update H's rows.
     od
end proc
:

Please verify my results very carefully for some small matrices. As I said, I'm not familiar with this material. Nor have I ever performed such a delicate "matrix surgery": appending a row and a column, inplace, on a square matrix. This surgery is the last 5 lines of the loop, and I hope that that number of lines can be reduced.

Since the number of rows/columns that need to be updated is constant, 3, regardless of the input size, this represents a reduction of the asymptotic time complexity of the original published Paar's algorithm by a whole power. Thus, this might be worth publishing.

@Magma The Maple command RowDimension is poorly named. It is the number of rows, not the dimension of the rows, which would be the number of columns. In your code, you set n to the number of rows. I believe that that's incorrect, and it should be the number of columns. Please confirm.

@mwahab The error is a result of an obscure Maple version difference that I forgot about (so obscure that I've never seen it documented). Change the table line to

Cfs:= table([terms]=~ [Cfs]):

It turns out that doing the entire matrix product H:= R^%T.R is much faster than doing just the dot products needed to get the above-the-diagonal entries. So, this version does the same 64x64 matrix in 1 minute.

Paar:= proc(M::Matrix)
local 
     R:= rtable(M, datatype= integer[8]), n:= upperbound(M)[2], xors:= add(R)-n, 
     H, hmax, maxij, lastcol
;
     for lastcol from 1+n do
          H:= rtable(antisymmetric, R^%T.R, datatype= integer[8]);
          hmax:= max(H);
          if hmax <= 1 then return R, xors fi;
          member(hmax, H, 'maxij');
          R(.., lastcol):= R[.., maxij[1]] *~ R[.., maxij[2]];
          R[.., [maxij]]:= R[.., [maxij]] *~  (1 -~ R[.., [-1$2]]);
          xors:= xors - (hmax-1)
     od
end proc
:

Further optimization is likely possible. H does not need to be entirely recomputed each iteration because only 2 columns/rows have changed and 1 needs to be added. This type of inplace updating of a growing matrix is tricky, and it'll take me a few days to get to it.

@Magma I modified the code so that it returns the updated matrix followed by the xor count. For my first example with n=64, I started with a random binary 64x64 matrix. The result was 324 columns, the xor count was 784, and it took 16 minutes to compute. If I declare the matrices R and H with datatype= integer[8], I can do it in slightly less than half that time. What times were you seeing?

Paar:= proc(M::Matrix)
local R:= copy(M), H, hmax, maxij, lastcol, n:= upperbound(M)[2], xors:= add(M)-n;
     for lastcol from n do
          H:= rtable(antisymmetric, (1..lastcol)$2, (i,j)-> R[..,i]^%T.R[..,j]);
          hmax:= max(H);
          if hmax <= 1 then return R, xors fi;
          member(hmax, H, 'maxij');
          R(.., lastcol+1):= R[.., maxij[1]] *~ R[.., maxij[2]];
          R[.., [maxij]]:= R[.., [maxij]] *~  (1 -~ R[.., [-1$2]]);
          xors:= xors - (hmax-1)
     od
end proc
:
M:= LinearAlgebra:-RandomMatrix(2^6 $ 2, generator= rand(0..1)):
CodeTools:-Usage(Paar(M));

Note that the dot product of binary vectors equals the Hamming weight of their conjunction. I've been guessing that due to processor optimizations, dot products over the integers are faster than counting the Hamming weights the other way.

@mwahab My code needs some small modifications to handle product (or "multiple") terms. And here it is:

Jxy:= [indets(eqn, specindex(op~(0,Jets)))[]]; #just an example

Eqn:= collect((lhs-rhs)(eqn), Jxy, 'distributed'):
Cfs:= frontend(coeffs, [Eqn], [{`+`, `*`}, {}], Jxy, 'terms'):
Cfs:= table([terms=~ Cfs]):

The first line will vary depending on the type of terms that you're looking for. I recommend that you use this new version regardless of whether there will be product terms. The frontend command is needed to prevent coeffs from looking inside the F, P, etc. terms.

If you can't bracket the N as requested above, could you simply provide the smallest N that you can find for which the results are undesirable and the largest N that you can find for which the results are good?

@mwahab What you want just involves slight modifications of the indets command. I think that your question #1 means that you want the 4th derivatives of u and v. Then

Jets:= {u,v}(x,y,z,t);
Jxy:= indets(eqn, myindex~(op~(0, Jets), identical(op~(Jets)[])$4));

The only difference is $4 instead of $2. If you want all derivatives of u and v, then

Jxy:= indets(eqn, specindex(op~(0,Jets)));

@mwahab What do you mean by the phrase "and their multiples", which you've used several times? Please show a specific example where a derivative has "multiples".

@Ramanujan Could you do some manual binary search experimentation to find the largest N for which you get good results, or conversely, the smallest N for which you get bad results?

Your series command is meaningless unless eta is expressed in terms of q. In order to help you, we need to see that expression.

@Magma Change the + in the exponent to %T

@Carl Love Here's a better way to automatically generate the numeric evaluation procedure for a first-order recurrence. This doesn't need a remember table because it always makes only one reference to the previous member of the sequence.

Makeproc:= proc(Reqn::`=`, IC::function(integer)=complexcons, Depvar::function(name))
local 
     A:= op(0, Depvar), n:= op(Depvar), An,
     P:= subs(
          _R= unapply(subs(A(n-1)= An, solve(Reqn, A(n))), An),
          (n::integer)-> _R(thisproc(n-1))
     )
; 
     assign(subs(A= P, IC));
     eval(P)
end proc
: 
SP:= Makeproc(S(n)-S(n-1) = 2*S(n)^2/(2*S(n)-1), S(1)=1, S(n)):
SP(100) - SP(99); 
SP(10000) - SP(9999);

@Magma Here is code retrofitted to not use max[index], so I believe it'll work in Maple 15:

Paar:= proc(M::Matrix)
local R:= copy(M), H, hmax, maxij, lastcol;
     for lastcol from 1+upperbound(M)[2] do
          H:= rtable(antisymmetric, R^+.R);
          hmax:= max(H);
          if hmax <= 1 then return R fi;
          member(hmax, H, 'maxij');
          R(.., lastcol):= R[.., maxij[1]] *~ R[.., maxij[2]];
          R[.., [maxij]]:= R[.., [maxij]] *~  (1 -~ R[.., [-1$2]])
     od
end proc: