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;

The sphere can be tangent to (some of) the edges in the exterior of the tetrahedron. In this case, Kitonum's conditions must be adjusted.

with(plots):with(plottools):
display(tetrahedron([[0.,0.,0.],[3.,0.,0.],[-4.,3.,0.],[-2.502782406,1.,4.628518959]],transparency=0.2),
        sphere([7.807886553,23.42365966,8.196378356],24.81629404,transparency=0.8,color=yellow),axes=none,
        orientation=[60,70,150]);

not(symbol) in this context is the type of an expression obtained by applying the function not to a symbol.
Examples of this type are not(x), not(Pi).

The type  "not a symbol"  would be Not(symbol).
type("abc", Not(symbol)); # true
type("abc", not(symbol)); # false
 

restart;
H:= k -> < cos(2*k/n*Pi), sin(2*k/n*Pi) >:
n:=100; dt:=0.005; t:=0;
P:=proc(k) global t; Threads:-Sleep(t);t:=t+dt; 
  plots:-display(plot([cos,sin,0..2*Pi]), plots:-arrow(<0,0>,H(k))) end:
Explore(P(k), k=1..n, animate, numframes=n);

Q:=proc(k) 
  plots:-display(plot([cos,sin,0..2*Pi]), plots:-arrow(<0,0>,H(k))) end:
Explore(Q(sqrt(k)), k=1..n^2, animate, numframes=5*n);

This is straightforward (unlike your previous question).

asympt(int((1-t^2)^n, t=0..1),n);

    sqrt(Pi)*sqrt(1/n)/2 + ...