C:=(m::nonnegint,n::integer) -> 
  coeff(coeff(product((1-q^i)*(1-q^i/z)*(1-q^(i-1)*z), i=1..m+1),q,m),z,n):

C(6,2),  C(6,4); 

      0,  1

(k is an arbitrary integer)

restart;
f:=(x,c,d)->(x-c)*(x+c)*(x-d)*(x+d):
g:=(x,c,d)->(x-c)*(x+c-1)*(x-d)*(x+d-1):
N:=10:
for c to N do  for d from c+1 to N do
  for x1 to c-1 do
    for F in [f,g] do
      if nops(factor(F(x,c,d)-F(x1,c,d)))=4 then lprint(F(x+k,c,d),m=F(x1,c,d)) fi;
    od
  od
od od:

(x+k-4)*(x+k+3)*(x+k-5)*(x+k+4), m = 180
(x+k-4)*(x+k+3)*(x+k-6)*(x+k+5), m = 360
(x+k-4)*(x+k+4)*(x+k-7)*(x+k+7), m = 720
(x+k-5)*(x+k+4)*(x+k-7)*(x+k+6), m = 504
(x+k-5)*(x+k+4)*(x+k-9)*(x+k+8), m = 1260
(x+k-5)*(x+k+4)*(x+k-10)*(x+k+9), m = 1800
(x+k-5)*(x+k+5)*(x+k-10)*(x+k+10), m = 2016
(x+k-6)*(x+k+5)*(x+k-7)*(x+k+6), m = 1260
(x+k-6)*(x+k+6)*(x+k-7)*(x+k+7), m = 1440
(x+k-6)*(x+k+5)*(x+k-9)*(x+k+8), m = 1080
(x+k-7)*(x+k+7)*(x+k-9)*(x+k+9), m = 2880
(x+k-7)*(x+k+6)*(x+k-10)*(x+k+9), m = 3780
(x+k-8)*(x+k+8)*(x+k-9)*(x+k+9), m = 5040
(x+k-9)*(x+k+8)*(x+k-10)*(x+k+9), m = 5544
(x+k-9)*(x+k+8)*(x+k-10)*(x+k+9), m = 2520
 

1. In Maple 2020 it can be done because a (statement)  is an expression

restart;
x:=10: str:="A":
str:=cat(str,  `if`(x=10,  [(x:=11)," it was 10"][2], [(x:=8)," it was not 10"][2]));  x;

2. `if`  is alias for ifelse  and a link appears in the help page of if. Or, use ?ifelse 

p:=randpoly([x,y,z]);
`%+`(sort([op(p)], key=abs@coeffs)[]);

Use value(%)  to  go back.

This isthe well known  Pell's equation.

S:=isolve(61*x^2+1=y^2):
for _Z1 from 0 to 6 do
lprint(eval([x,y],rationalize(S[4]))[])
od;

0, 1
226153980, 1766319049
798920165762330040, 6239765965720528801
2822295814832482312327709940, 22042834973108102061352541449
9970149719303180503641083029374964080, 77869358613928486808166555366140995201
35220930741174421456911021812718768924061809900, 275084262906388245923976756042747916825335226249
124422801783292138491822391332416163557158135530198606120, 971773147303355325052564141449134520779147876502526039201
 

It is easy to obtain such equalities. Here is a simple example. It can be extended to a generic one.
Why do you consider them so important?

g:=proc(z)
   local m:=subs(sqrt(3)=0,z),
         n:=z-m;
   m^2+n^2-1+ 2*m*n-sqrt(3)
end: 
F:=c->sqrt(1+sqrt(3)+c):
G:= n -> (F@@n)(g@@n)(1+sqrt(3))=1+sqrt(3):
G(3);G(4);G(5);

You do not have a value (true or false) for print_table
Note also that a better  header for your proc is:

proc(f::procedure, a::realcons, b::realcons, N::posint, print_table::truefalse)
 

You must use parallel substitutions in ex:
ex := simplify(subs([ X = cos(theta)*X-sin(theta)*Y, Y = sin(theta)*X+cos(theta)*Y ], eq)); 

(You have defined the procedure f. It is then easier to use f(X-9/7, Y+8/7)  and similar for ex, instead of subs) .

It is interesting that for diff, the first argument can be almost anything, including a mathematical nonsense. For sets, lists, rtables it acts as expected (elementwise), but even strings are accepted. 
Note that series may produce errors.

f:=sin(x) + x^3*"a";
g:=sin(x) + [cos(x), x^3];
h:=sin(x) + [cos(x)*"a"] + {x^2*"b"^3 + "c"};

You can reduce the time to 0 because if a,b,... are >0 (as you have tried) then g must be a square.

The problem is generated by simplify.

(1 + (203808*exp(-(342569*t)/506))/131537)^(131537/203808);

simplify(%, sqrt);  # or simplify(%)

For example, with a smaller exponent:

simplify((1 + (203808*exp(-(342569*t)/506))/131537)^(3/7), sqrt); 

You should try to understand and explain what happens in the next lines:

restart;
P:=proc(a,b,c)
local S:=LinearAlgebra:-RandomMatrix(3, generator=rand(-10..10),shape=skewsymmetric),
      V0:=LinearAlgebra:-RandomVector(3, generator=rand(-10..10))^+,
      A:=(S-1)^(-1) . (S+1) . <a,0,0;0,b,0;0,0,c>;
A:=ilcm(entries(denom~(A),nolist))*A^+;
<V0, V0+A[1],V0+A[2],V0+A[3]>
end:
seq(P(3,5,14), 1..10);

Only a few miliseconds are needed for the 128 solutions.

restart;
ti:=time[real]():
XY:=seq(seq(`if`(irem(x*y,x+y)=0,[x,y],NULL), y=2..x-1),x=2..30):
N:=0:
for xy in XY do
  x,y:=op(xy): z:=x*y/(x + y):
  for a from 2 to 10 do
    for b from 2 to a-1 do
      r:=a/x  + b/y; m:=numer(r); n:=denom(r);
      if n<=10 and m>=2 and m<=n and igcd(a,b,m)=1 
        then N:=N+1; sol[N]:=[x, y, z, a, b, m, n] fi
od od od: time=time[real]()-ti;
'N'=N,entries(sol);

Edit: corrected a<=30 to a<=10. Now 64 solutions in 2 ms.

The radius of convergence can be computed this way:

f:=1/(1-x^2):
x0:=1:
S:=rhs~(`union`(singular(f, x))):
r:=min(seq((abs(x-x0), x in S minus {x0}))):
f=series(f, x=x0); abs(x-x0)<r;