The authors initiate a study of discrete analogues of certain non-translation invariant singular integral operators such as singular Radon transforms whose continuous versions include operators of the form \[ R(f)(x)=\text{ p.v.} \int_{\mathbb{R}^k} f(\gamma_t(x))K(t) dt (x\in \mathbb{R}^n) \] where \(\gamma_t\) is a family of diffeomorphisms satisfying some non-flatness condition. A discrete version that reflects this non-degeneracy is \[ R_d(f)(n,n')=\sum_{m\in\mathbb{Z}^k, m\neq 0} f(n-m,n'-Q(n,m))K(m) ((n,n')\in\mathbb{Z}^k\times \mathbb{Z}^{\ell}) \] where \(Q:\mathbb{Z}^k\times\mathbb{Z}^k\to\mathbb{Z}^{\ell}\) is the restriction of a polynomial and the kernel \(K\) has a Fourier transform that agrees with a bounded function away from the origin. By Plancherel’s theorem, \(\ell^2(\mathbb{Z}^N)\)-boundedness of \(R_d\) for \(N=k+\ell\) follows from \(\ell^2(\mathbb{Z}^k)\)-boundedness of the operator \[ T(f)(n)=\sum_{m\in\mathbb{Z}^k, m\neq n} e^{2\pi i P(n,m)}K(n-m)f(m). \] To deduce the boundedness of \(R_d\) it is crucial that the bound for \(T\) does not depend on the coefficients of \(P\), since in the case of \(R_d\) \(P\) takes the form \(P(n,m)=\sum_{j=1}^l \xi_j Q_j(n,n-m)\). The main result of the paper is thus the \(\ell^2(\mathbb{Z}^k)\)-boundedness of the operator \(f\mapsto T(f)\), with a bound depending on \(k\) and the degree of \(P\), but not otherwise on its coefficients. In fact, this independence of the operator norm on the coefficients is the main difficulty. The detailed analysis requires first breaking down the kernel \(K\) by a partition of unity over dyadic annuli. A focus on the case \(k= 1\) and \(P(n,m)= \xi Q(n,n-m)\) reduces the problem to considering the possible ways in which \(\xi \in \mathbb{T}\) can be well-approximated by a rational number whose denominator is, when reduced to lowest terms, of size proportional to a power of the radius of a given annulus. The analysis then reduces to estimates of classical Gauss and Weyl sums, except that for the general problem the sums are over multidimensional objects must be regarded as being operator-valued, and must apply to appropriate shifts of \(P\). These complications aside, in the univariate case one must estimate \(S= \sum_{1\leq n\leq r} e^{2\pi i P(n)}\) where \(P(x)= \sum_{k=1}^{\text{ deg}P}\theta_k x^k\). Assuming that for some \(k\in\{1,\dots,\text{deg}P\}\) one has a rational \(p/q\) in lowest terms such that \(|\theta_k-p/q|<1/q^{2}\) and for some \(\epsilon>0\), \(r^\epsilon\leq q\leq r^{k-\epsilon}\), then there is a \(\delta>0\) such that \(|S|\leq O(r^{1-\delta})\) in which the bound depends only on the degree of \(P\). Appropriate assembly of estimates of this type ultimately yields the boundedness of \(T\).