intersectionArrayToSpectrum:=function(a,F)
 d:=#a div 2;
 b0:=a[1];
 bi:=func<i|a[i+1]>;
 ci:=func<i|a[d+i]>;
 ai:=func<i|b0-bi(i)-ci(i)>;
 R<x>:=PolynomialRing(F);
 v:=[R|1,x];
 vi:=func<i|v[i+1]>;
 n:=1;
 K:=[1,b0];
 for i in [1..d-1] do
  Append(~v,1/ci(i+1)*(x*v[i+1]-ai(i)*v[i+1]-bi(i-1)*v[i]));
  Append(~K,K[i+1]*bi(i)/ci(i+1));
 end for;
 n:=&+K;
 th:=[b0] cat [r[1]:r in Roots(&+[v[i+1]:i in [0..d]])];
 mtp:=[n/(&+[Evaluate(v[i+1],t)^2/K[i+1]:i in [0..d]]):t in th];
 return [[th[i],mtp[i]]:i in [1..d+1]];
end function;

F:=Rationals();
a:=[3,2,1,1,2,3];
intersectionArrayToSpectrum(a,F) eq [[3,1],[-3,1],[-1,3],[1,3]];
F<s>:=QuadraticField(13);
a:=[3,2,1,1,1,3];
intersectionArrayToSpectrum(a,F) eq [[3,1],[-1,3],
 [1/2*(-s-1),1/13*(-4*s+52)],[1/2*(s-1),1/13*(4*s+52)]];