diff --git a/cytools/__init__.py b/cytools/__init__.py
index d56fff2..a353ac7 100644
--- a/cytools/__init__.py
+++ b/cytools/__init__.py
@@ -19,7 +19,7 @@
 from cytools.utils import read_polytopes, fetch_polytopes
 
 # Latest version
-version = "0.3.3"
+version = "0.3.4"
 versions_with_serious_bugs = []
 
 # Check for more recent versions of CYTools
diff --git a/cytools/calabiyau.py b/cytools/calabiyau.py
index fa58456..0dcafbc 100644
--- a/cytools/calabiyau.py
+++ b/cytools/calabiyau.py
@@ -918,8 +918,7 @@ def glsm_charge_matrix(self, include_origin=True):
         """
         if self._glsm_charge_matrix is not None:
             return np.array(self._glsm_charge_matrix)[:,(0 if include_origin else 1):]
-        toric_divs = [0]+list(self.prime_toric_divisors())
-        pts = self.polytope().points_to_indices(self.polytope().points()[toric_divs])
+        pts = [0]+list(self.prime_toric_divisors())
         self._glsm_charge_matrix = self.polytope().glsm_charge_matrix(
                                             include_origin=True,
                                             points=pts)
@@ -956,8 +955,7 @@ def glsm_linear_relations(self, include_origin=True):
         """
         if self._glsm_linrels is not None:
             return np.array(self._glsm_linrels)[(0 if include_origin else 1):,(0 if include_origin else 1):]
-        toric_divs = [0]+list(self.prime_toric_divisors())
-        pts = self.polytope().points_to_indices(self.polytope().points()[toric_divs])
+        pts = [0]+list(self.prime_toric_divisors())
         self._glsm_linrels = self.polytope().glsm_linear_relations(
                                 include_origin=True,
                                 points=pts)
@@ -1004,12 +1002,11 @@ def divisor_basis(self, include_origin=True, as_matrix=False):
         ```
         """
         if self._divisor_basis is None:
+            pts = [0]+list(self.prime_toric_divisors())
             self.set_divisor_basis(self.polytope().glsm_basis(
                                    integral=True,
                                    include_origin=True,
-                                   points=self.polytope().points_to_indices(
-                                        self.polytope().points()[[0]+list(self.prime_toric_divisors())]))
-                                    )
+                                   points=pts))
         if len(self._divisor_basis.shape) == 1:
             if 0 in self._divisor_basis and not include_origin:
                 raise Exception("The basis was requested not including the "
diff --git a/cytools/polytope.py b/cytools/polytope.py
index 4324a2e..41d6be2 100644
--- a/cytools/polytope.py
+++ b/cytools/polytope.py
@@ -1517,13 +1517,14 @@ def facets(self):
         """
         return self.faces(self._dim-1)
 
-    def vertices(self):
+    def vertices(self, as_indices=False):
         """
         **Description:**
         Returns the vertices of the polytope.
 
         **Arguments:**
-        None.
+        - `as_indices` *(bool)*: Return the points as indices of the full
+          list of points of the polytope.
 
         **Returns:**
         *(numpy.ndarray)* The list of vertices of the polytope.
@@ -1542,6 +1543,8 @@ def vertices(self):
         ```
         """
         if self._vertices is not None:
+            if as_indices:
+                return self.points_to_indices(self._vertices)
             return np.array(self._vertices)
         if self._dim == 0:
             self._vertices = np.array([self._input_pts[0]])
@@ -1607,6 +1610,8 @@ def vertices(self):
                 if pt_tup not in input_pts:
                     input_pts.append(pt_tup)
             self._vertices = np.array([list(pt) for pt in input_pts if pt in tmp_vert])
+        if as_indices:
+            return self.points_to_indices(self._vertices)
         return np.array(self._vertices)
 
     def dual(self):
@@ -1787,24 +1792,26 @@ def glsm_charge_matrix(self, include_origin=True,
         # If the result is not cached we do the computation
         # We start by finding a basis of columns
         linrel = self.points()[list(pts_ind)].T
-        if int(round(np.linalg.det(np.array(fmpz_mat(linrel.tolist()).snf().tolist(), dtype=int)[:,:4]))) != 1:
-            raise Exception("The points do not generate the lattice.")
+        sublat_ind =  int(round(np.linalg.det(np.array(fmpz_mat(linrel.tolist()).snf().tolist(), dtype=int)[:,:linrel.shape[0]])))
         norms = [np.linalg.norm(p,1) for p in linrel.T]
         linrel = np.insert(linrel, 0, np.ones(linrel.shape[1], dtype=int), axis=0)
         good_exclusions = 0
         basis_exc = []
         indices = np.argsort(norms)
+        indices[:linrel.shape[0]] = np.sort(indices[:linrel.shape[0]])
         for n_try in range(4):
             if n_try == 1:
                 indices[:] = np.array(range(linrel.shape[1]))
             elif n_try > 1:
-                np.shuffle(indices[1:])
+                np.random.shuffle(indices[1:])
+                indices[:linrel.shape[0]] = np.sort(indices[:linrel.shape[0]])
             for ctr in range(np.prod(linrel.shape)+1):
                 found_good_basis=True
                 ctr += 1
                 if ctr > 0:
                     st = max([good_exclusions,1])
                     indices[st:] = np.roll(indices[st:], -1)
+                    indices[:linrel.shape[0]] = np.sort(indices[:linrel.shape[0]])
                 linrel_rand = np.array(linrel[:,indices])
                 try:
                     linrel_hnf = fmpz_mat(linrel_rand.tolist()).hnf()
@@ -1813,12 +1820,14 @@ def glsm_charge_matrix(self, include_origin=True,
                 linrel_rand = np.array(linrel_hnf.tolist(), dtype=int)
                 good_exclusions = 0
                 basis_exc = []
+                tmp_sublat_ind = 1
                 for v in linrel_rand:
                     for i,ii in enumerate(v):
                         if ii != 0:
                             if integral:
-                                if abs(ii) == 1:
-                                    v *= ii
+                                tmp_sublat_ind *= abs(ii)
+                                if sublat_ind % tmp_sublat_ind == 0:
+                                    v *= ii//abs(ii)
                                     good_exclusions += 1
                                 else:
                                     found_good_basis = False
@@ -1829,7 +1838,7 @@ def glsm_charge_matrix(self, include_origin=True,
                 if found_good_basis:
                     break
             if found_good_basis:
-                    break
+                break
         if not found_good_basis:
             raise Exception("failed")
             if ctr == np.prod(linrel.shape):
@@ -1845,16 +1854,16 @@ def glsm_charge_matrix(self, include_origin=True,
                                                points=points, integral=False)
         linrel_dict = {ii:i for i,ii in enumerate(indices)}
         linrel = np.array(linrel_rand[:,[linrel_dict[i] for i in range(linrel_rand.shape[1])]])
-        basis_ind = np.array(sorted([i for i in range(linrel.shape[1]) if linrel_dict[i] not in basis_exc]), dtype=int)
-        basis_exc = [i for i in range(linrel.shape[1]) if i not in basis_ind]
+        basis_ind = np.array([i for i in range(linrel.shape[1]) if linrel_dict[i] not in basis_exc], dtype=int)
+        basis_exc = np.array([indices[i] for i in basis_exc])
         glsm = np.zeros((linrel.shape[1]-linrel.shape[0],linrel.shape[1]), dtype=int)
-        glsm[:,basis_ind] = np.eye(len(basis_ind),dtype=int)
-        for nb in basis_exc:
+        glsm[:,basis_ind] = np.eye(len(basis_ind), dtype=int)
+        for nb in basis_exc[::-1]:
             tup = [(k,kk) for k,kk in enumerate(linrel[:,nb]) if kk]
-            if len(tup) != 1 or abs(tup[0][1]) != 1:
+            if sublat_ind % tup[-1][1] != 0:
                 raise Exception("Problem with linear relations")
-            i,ii = tup[0]
-            glsm[:,nb] = -ii*glsm.dot(linrel[i])
+            i,ii = tup[-1]
+            glsm[:,nb] = -glsm.dot(linrel[i])//ii
         # Check that everything was computed correctly
         if (np.linalg.matrix_rank(glsm[:,basis_ind]) != len(basis_ind)
                 or any(glsm.dot(linrel.T).flat)
diff --git a/cytools/utils.py b/cytools/utils.py
index f340d9c..6d58e5c 100644
--- a/cytools/utils.py
+++ b/cytools/utils.py
@@ -598,12 +598,13 @@ def set_divisor_basis(tv_or_cy, basis, include_origin=True):
         linrels_new[:,b] = linrels_tmp[:,len(nobasis):]
         self._curve_basis_mat = np.zeros(glsm_cm.shape, dtype=int)
         self._curve_basis_mat[:,b] = np.eye(len(b),dtype=int)
-        for nb in nobasis:
+        sublat_ind =  int(round(np.linalg.det(np.array(fmpz_mat(linrels.tolist()).snf().tolist(), dtype=int)[:,:linrels.shape[0]])))
+        for nb in nobasis[::-1]:
             tup = [(k,kk) for k,kk in enumerate(linrels_new[:,nb]) if kk]
-            if len(tup) != 1 or abs(tup[0][1]) != 1:
+            if sublat_ind % tup[-1][1] != 0:
                 raise Exception("Problem with linear relations")
-            i,ii = tup[0]
-            self._curve_basis_mat[:,nb] = -ii*self._curve_basis_mat.dot(linrels_new[i])
+            i,ii = tup[-1]
+            self._curve_basis_mat[:,nb] = -self._curve_basis_mat.dot(linrels_new[i])//ii
     # Else if input is a matrix
     elif len(b.shape) == 2:
         if not config._exp_features_enabled:
@@ -638,11 +639,12 @@ def set_divisor_basis(tv_or_cy, basis, include_origin=True):
                                 points=self.polytope().points_to_indices(self.triangulation().points()))
         self._divisor_basis_mat = np.array(new_b)
         nobasis = np.array([i for i in range(glsm_cm.shape[1]) if i not in standard_basis])
-        for nb in nobasis:
+        sublat_ind =  int(round(np.linalg.det(np.array(fmpz_mat(linrels.tolist()).snf().tolist(), dtype=int)[:,:linrels.shape[0]])))
+        for nb in nobasis[::-1]:
             tup = [(k,kk) for k,kk in enumerate(linrels[:,nb]) if kk]
-            if len(tup) != 1 or abs(tup[0][1]) != 1:
+            if sublat_ind % tup[-1][1] != 0:
                 raise Exception("Problem with linear relations")
-            i,ii = tup[0]
+            i,ii = tup[-1]
             for j in range(self._divisor_basis_mat.shape[0]):
                 self._divisor_basis_mat[j] -= self._divisor_basis_mat[j,nb]*linrels[i]
         # Finally, we invert the matrix and construct the dual curve basis
@@ -655,12 +657,12 @@ def set_divisor_basis(tv_or_cy, basis, include_origin=True):
             inv_mat *= -1;
         self._curve_basis_mat = np.zeros(glsm_cm.shape, dtype=int)
         self._curve_basis_mat[:,standard_basis] = np.array(inv_mat).T
-        for nb in nobasis:
+        for nb in nobasis[::-1]:
             tup = [(k,kk) for k,kk in enumerate(linrels[:,nb]) if kk]
-            if len(tup) != 1 or abs(tup[0][1]) != 1:
+            if sublat_ind % tup[-1][1] != 0:
                 raise Exception("Problem with linear relations")
-            i,ii = tup[0]
-            self._curve_basis_mat[:,nb] = -ii*self._curve_basis_mat.dot(linrels[i])
+            i,ii = tup[-1]
+            self._curve_basis_mat[:,nb] = -self._curve_basis_mat.dot(linrels[i])//ii
         self._curve_basis = np.array(self._curve_basis_mat)
     else:
         raise Exception("Input must be either a vector or a matrix.")
diff --git a/scripts/linux/cytools b/scripts/linux/cytools
index ac1cc24..c663080 100755
--- a/scripts/linux/cytools
+++ b/scripts/linux/cytools
@@ -106,7 +106,7 @@ cat << EOF
  ░░█████████     █████       █████   ░░██████ ░░██████  █████ ██████
   ░░░░░░░░░     ░░░░░       ░░░░░     ░░░░░░   ░░░░░░  ░░░░░ ░░░░░░
 
-        Developed by Liam McAllister's Group | Version 0.3.3
+        Developed by Liam McAllister's Group | Version 0.3.4
              https://cytools.liammcallistergroup.com
 
 EOF
diff --git a/scripts/macos/cytools b/scripts/macos/cytools
index 4f448b0..b6fb687 100755
--- a/scripts/macos/cytools
+++ b/scripts/macos/cytools
@@ -100,7 +100,7 @@ cat << EOF
  ░░█████████     █████       █████   ░░██████ ░░██████  █████ ██████
   ░░░░░░░░░     ░░░░░       ░░░░░     ░░░░░░   ░░░░░░  ░░░░░ ░░░░░░
 
-        Developed by Liam McAllister's Group | Version 0.3.3
+        Developed by Liam McAllister's Group | Version 0.3.4
              https://cytools.liammcallistergroup.com
 
 EOF
diff --git a/scripts/windows/launcher.ps1 b/scripts/windows/launcher.ps1
index 015b3ee..87ff5e5 100644
--- a/scripts/windows/launcher.ps1
+++ b/scripts/windows/launcher.ps1
@@ -38,7 +38,7 @@ $banner=@"
  ░░█████████     █████       █████   ░░██████ ░░██████  █████ ██████
   ░░░░░░░░░     ░░░░░       ░░░░░     ░░░░░░   ░░░░░░  ░░░░░ ░░░░░░
 
-        Developed by Liam McAllister's Group | Version 0.3.3
+        Developed by Liam McAllister's Group | Version 0.3.4
              https://cytools.liammcallistergroup.com
 
 "@
diff --git a/setup.py b/setup.py
index 2c8f6d0..1f85daa 100644
--- a/setup.py
+++ b/setup.py
@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="cytools",
-    version="0.3.3",
+    version="0.3.4",
     author="Liam McAllister Group",
     author_email="",
     description="A software package for analyzing Calabi-Yau hypersurfaces in toric varieties.",