Asking SnapPy for core curves after surgery Suppose I give SnapPy a cusped hyperbolic 3-manifold (using, say, the link editor) and specify some filling. SnapPy can then provide a presentation of the fundamental group of the filled manifold. Can it tell me what the core curve of the added solid torus is, as a word in the fundamental group?
 A: I think you want to use Snappy's 'fillings_may_affect_generators=False' flag for your purposes. Although you can also extract the information directly if for some reason this becomes inefficient.
In[1]: M = Manifold('m004')

not needed but M.fundamental_group? will give all possible flags
In[2]: G = M.fundamental_group?
Docstring:
Manifold.fundamental_group(self, simplify_presentation=True, fillings_may_affect_generators=True, minimize_number_of_generators=True, try_hard_to_shorten_relators=True)

there is more but it is redacted here
In[3]: G = M.fundamental_group()

In[4]: G
Out[4]: 
Generators:
   a,b
Relators:
   aaabABBAb

In[5]: m=G.meridian(); G.meridian()
Out[5]: 'ab'

In[6]: l=G.longitude(); G.longitude()
Out[6]: 'aBAbABab'

In[7]: M.dehn_fill((5,1),0)

In[8]: H = M.fundamental_group(fillings_may_affect_generators=False)

In[9]: H
Out[9]: 
Generators:
   a,b
Relators:
   aaabABBAb
   ababababababaBAbAB

Here, the second relation is m^5*l. 
The core curve will be isotopic to any curve p*[m]+q*[l] (in boundary M) such that |5q-1p|=1. More generally for filling along r,s, we want |rq-sp|=1.
There are a number of simple python scripts to do that for example the extended gcd script taken from https://www.kkhaydarov.com/greatest-common-divisor-python/
def egcd(r, s):
   if r == 0:
     return (s, 0, 1)
   else:
     g, x, y = egcd(s % r, r)
     return (g, y - (b // a) * x, x)

Here the g is the gcd of r and s, p=y - (s // r) * x and q=x.  
To complete the example where we will fill along (5,1), (p,q)=(1,0) so m is sufficient.  
