You can use the  Student:-LinearAlgebra:-RotationMatrix  command to parameterize a circle in 3D. In the example below  M  is rotation matrix around the axis defined by the vector  <2,3,-4>  by the angle  t ,  A  is the center of the circle,  B  is the point on this circle: 

restart;
M:=Student:-LinearAlgebra:-RotationMatrix(t, <2,3,-4>);
A:=<1,2,3>: 
B:=<2,4,5>:
Eq:=A + M.(B-A); # The parametric equations of the circle
plots:-spacecurve(Eq, t=0..2*Pi, color=red, axes=normal, scaling=constrained, view=[-2..5.5,-2..5.5,-2..5.5]);

 
Note that the vectors  <2, 3, -4>  and  B - A  are perpendicular.

restart;
y := t -> 0.7*sin(300*t);
filenm := "https://www.gw-openscience.org/GW150914data/P150914/fig2-unfiltered-template-reconstruction-H.txt";
K := ImportMatrix(filenm, source = Matlab, skiplines = 1);
R:=[seq([K[i,1],abs(K[i,3]-y(K[i,1]))], i=1..3441)]:
plot([K[() .. (), [1, 3]], y(t)], t = 0.35 .. 0.4, color=[red,blue]);
plot(R, color=green, view=[0.35..0.4,0..1]);

I understood the plot of residuals as the plot of the modulus of the difference of two plots (blue and red).

Edit. If you just take the differences (no modulus), you get a graph like Carl's.

R:=[seq([K[i,1],y(K[i,1])-K[i,3]], i=1..3441)]:
plot(R, color=green, view=[0.35..0.4,-1.5..1.5]);

When you use  break, you only get out of one innermost loop (the rest continue to work).
To solve the problem, just replace two lines of code 
print(dedges);
break;
with one line
return dedges;

Edit. For any number of pairs of edges, you can always limit yourself to one cycle. Carl has already answered your second question. Below is another method that does not use the Iterator package. This method is suitable also for older versions of Maple:

restart;
with(GraphTheory):
with(SpecialGraphs):
g:=CompleteGraph(6);
g1:=DeleteEdge(g, {3,2}, inplace = false);

Conremovedges:= proc(G::Graph,Paircrossedges)
local n, P, p, dedges;
uses combinat;
n:=nops(Paircrossedges);
P:=permute([1$n,2$n],n);
for p in P do
dedges:={seq(Paircrossedges[i][p[i]], i=1..n)}; 
if EdgeConnectivity(DeleteEdge(G, dedges, inplace = false))=3 then
return dedges;
end if;
end do:
"not found";
end proc:

# Example of use:
Paircrossedges:=[[{1,6},{3,5}],[{2,5},{1,4}],[{2,6},{3,4}]];
Conremovedges(g1,Paircrossedges);

If you want the square root not to be grayed out when using the inert form, but to look as usual, then use the  InertForm  package with the option  inert=false :

InertForm:-Display(1/%sqrt(2), inert=false);

If you just need a procedure that returns the distance from a specific vertex  v  to a specific edge  e  in a graph  G , then here it is:

restart;
Dist:=proc(v,e,G)
uses GraphTheory;
min(Distance(G,v,e[1]), Distance(G,v,e[2]));
end proc:

Examples of use:

G := GraphTheory:-SpecialGraphs:-PetersenGraph():
E := convert(GraphTheory:-Edges(G),list);
Dist(1,{8,9},G); # Distance from the vertix 1 to the edge {8,9}
Dist~(1,E,G); # List of distances from vertex 1 to all edges from the list E 

In fact, the  plot is shown incorrectly  and  pp(2020) is calculated incorrectly . The coefficients of the polynomial differ too much from each other. Simply increasing  Digits  will not help here. It is unclear to what extent  Digits  should be increased. Conversion of all constants into exact fractions radically solves the problem.. Then all calculations are done exactly:

restart;
a := 0.340699639252428*10^13;
b := -0.118229737138742*10^11;
c := 0.175816922773262*10^8;
d := -14523.7138112711;
e := 7.19782673654200;
f := -0.00214008825847444;
g := 0.353463615941623*10^(-6);
h := -0.250170509163729*10^(-10);
a,b,c,d,e,f,g,h:=convert~([a,b,c,d,e,f,g,h],fraction)[];

pp := x -> a + b*x + c*x^2 + d*x^3 + e*x^4 + f*x^5 + g*x^6 + h*x^7;

plot(pp(x), x = 1960 .. 2100, gridlines = true, size = [800, 400]);
pp(2020);
evalf(%);

Now all the contradictions between the values of the function  pp  on the plot and the values of the function  pp  calculated at separate points  ( like pp(2020)) disappear.

Example:

restart;
myTable:=table(["a"=1,"b"=-1,"c"=1]);
select(i->myTable[i[]]=1,{indices(myTable)});

                                               {["a"], ["c"]}

Edit.

If we define  f  as a procedure, then we get a more convenient and compact syntax for specific calculations: 

restart;
with(orthopoly):
f:=unapply(convert(diff( P(n, cos(theta)), theta, theta),diff), n):
int( f(2), theta=0..Pi/4);
int( f(3), theta=0..Pi/4);

As an alternative you can use  diff  instead of  D :

restart;
g:=x->x^(1/3);
f:=x->convert(g(x),surd);
diff(f(x),x);

It's very simple:

restart;
A:=[1,2,6,7,9,10,15,17]:
A[3..5];

                                             [6, 7, 9]

restart;
f:=(a*x1+b*x2)/x3+c*x1^2/(e*x2+d);
normal(f);
numer(%);

Or

restart;
(a*x1+b*x2)/x3+c*x1^2/(e*x2+d)=f;
normal(%);
%*denom(lhs(%));

Code in Maple 2018.2 (I have no newer versions):

restart;
B:=sqrt( (-4*u^(1/3)+1)*u^(4/3));
A:=1/(-12*u+3*u^(2/3)-3*B);
res:=Int(A,u);
res1:=IntegrationTools:-Change(res,t=u^(1/3));
value(res1);
simplify(%) assuming t>0;
eval(%, t=u^(1/3)); 

You can get several more solutions from one solution if you use transformations that transform a cubic integer lattice into itself. They will be, for example, reflections about the coordinate planes, rotations around the coordinate axes at angles that are multiples of 90 degrees, and some others. But in terms of geometry, you end up with the same pair of vectors. New solutions can be easily obtained using the brute force method. First we find the list L and then from the list  L, the list  L1  is obtained, where none of the coordinates is equal to 0:

restart;
k:=0:
for a from 0 to 9 do
for b from a to 9 do
for c from b to 9 do
for x from -9 to 9 do
for y from -9 to 9 do
for z from -9 to 9 do
if 4*(a*x+b*y+c*z)^2=3*(a^2+b^2+c^2)*(x^2+y^2+z^2) and a^2+b^2+c^2<>0 and x^2+y^2+z^2<>0 and a*x+b*y+c*z>0 then k:=k+1;
L[k]:=[<a,b,c>,<x,y,z>] fi;
od: od: od: od: od: od:
L:=convert(L,list):  
nops(L);
L1:=select(t->`and`(seq(t[1][i]<>0,i=1..3),seq(t[2][i]<>0,i=1..3)),L):
nops(L1);

The very first pair of vectors in the lists  L  and  L1  is your original example, in which the first and third coordinates are reversed (reflection of space relative to the plane  x=z ). To reduce the search, I chose the vector  <a,b,c>  with non-decreasing coordinates.

This can be done in many ways. Here's one:

A := [1,0,3,4,5,5]: B := [0,1,2,6,3,6]:
nops(select(`>`, A - B, 0));

                                           3

The simple recursive procedure does the work:

restart;
a:=proc(n)
local m;
option remember;
if n=1 then return 1 elif n=2 then return 0 else
for m from 1 to n-2 do
if a(n-1)=a(n-1-m) then return m fi; od;
min(seq(abs(a(n-1)-a(k)),k=1..n-2)) fi;
end proc:

Examples of use:

seq(a(n), n=1..100);
plots:-listplot([%], style=point, symbol=solidcircle, color=red);

Edit. The procedure above works about 2 times faster than tomleslie's method.