Integer relations-I. Brun’s algorithm

Given a set of real numbers $\mathbf{x}^T=\{x_1,x_2, \cdots, x_n\}$ an integer relation $\{a_1, a_2, \cdots, a_n\}$ is a set of integers $a_i\in \mathbb{Z}$ such that

$a_1 x_1+a_2 x_2+\cdots+a_n x_n=0$

The PSLQ algorithm constructs such a set of numbers by creating a sequence of integer-valued matrices $B_j$ for which

$\mathbf{y}^T=\mathbf{x}^T\cdot B_j$

for which at least on component of $\mathbf{y}\in \mathbb{R}^n$ becomes zero to within some acceptable tolerance. If $y_k\approx 0$ for some $j$ then the $k^{th}$ column of $B_j$ is our integer vector $\mathbf{a}^T$ . The collapse of elements of $\mathbf{y}$ is rather abrupt, so one can be confident that it signals an actual relation.

Where do integer relations occur?
$\bullet$ Trig identities. Let $\theta$ be arbitrary, and construct

$\mathbf{x}^T=(\cos\theta, \cos(3\theta), \cos(5\theta), \cos^5\theta)$

These obey an integer relation with $\mathbf{a}=(10,5,1,-16)$ .
$\bullet$ Algebraic numbers. The real number $z=3.650281539872884745=\sqrt{7+\sqrt{40}}$ has a minimal polynomial. Let

$\mathbf{x}^T=(1,z,z^2,\cdots, z^n)$

then $0=9-14 z^2+z^4$ .\\
$\bullet$ Establishing formulas. Construct

$S_j=\sum_{k=0}^\infty {1\over (16)^k}{1\over 8k+j}, \quad k=1,2,\cdots,7$

and let

$\mathbf{x}^T=(\pi, S_1,S_2,S_3,S_4,S_5,S_6,S_7)$

In 1996 Simon Plouffe used the PSFQ algorithm to discover

$\pi=\sum_{k=0}^\infty {1\over (16)^k}\Big({4\over 8k+1}-{2\over 8k+4}-{1\over 8k+5}-{1\over 8k+6}\Big)$

which lets you compute any individual hexadecimal digit of $\pi$ without computing previous digits.
$\bullet$ Statistical mechanics (arXiv:cond-mat/9811264v3) . Nearest neighbor correlation functions of the hard-squares model exhibit integer relations

$\rho_2=\langle s_i s_j\rangle, \quad \rho'_2=\langle s_i s_k\rangle, \quad \rho_3=\langle s_i s_j\ s_k\rangle, \rho_4=\langle s_i s_j\ s_k s_m\rangle$

and $\rho=\langle s_p\rangle$ where $s_p$ is a lattice spin, $s_j$ is the spin adjacent to it on the right, $s_m$ is adjacent to $s_p$ on the left, $s_i$ is adjacent below, $s_k$ is adjacent above.

$6\rho-4\rho_2-2\rho'_2+3\rho_3-\rho_4=1$

There are simple algorithms such as Brun’s for finding integer relations, but the PSLQ algorithm is superior in every way, since it also gives you a limit on the size of a relation if one exists at all.
Brun’s algorithm is basically the Euclidean; let $\mathbf{x}^T=\{x_1, x_2, \cdots, x_n\}$ , and let $\mathbf{b}^T=\{\mathbf{b}_1,\mathbf{b}_2, \cdots, \mathbf{b}_n\}$ where $\mathbf{b}_i$ is a unit vector with all zero components except for the $i^{th}$ . Let $P$ be the permutation that sorts $\mathbf{x}$ into descending order

$\mathbf{x}'=P\cdot \mathbf{x}, \qquad x'_1\ge x'_2\ge \cdots \ge x'_n$

Apply this to $\mathbf{b}$ ;

$\mathbf{b}'=P\cdot \mathbf{b}$

Compute $r=[x'_1/x'_2]$ where $[ \, ]$ denotes rounding down, and redefine

$x'_1\rightarrow x'_1-r \, x'_2, \qquad \mathbf{b}'_1\rightarrow \mathbf{b}'_1-r \, \mathbf{b}'_2$

and repeat the process with a new sorting permutation. What this does is applies integer relation after integer relation (keeping track of them) until $x_n\approx 0$ , then the final $\mathbf{b}_n$ is the desired final integer relation.

; common lisp
; Brun's method for integer relations
; Make sure you use double-floats!
(defun Brun (lst iters)
(setf foo (copy-list lst))
(setf n (1- (length foo)))
(setf b (loop for i from 0 to n collect (loop for j from 0 to n collect 0)))
(loop for i from 0 to n do (setf (nth i (nth i b)) 1))

(loop for l from 1 to iters do
;; sort can destroy your list
(setf foo2 (sort (copy-list foo) #'>))
(setf y (loop for x in foo2 collect (position x foo :test #'equal)))
(setf bn (loop for x in y collect (nth x b)))
(rotatef foo foo2)
(rotatef b bn)
(setf q (truncate (/ (nth 0 foo) (nth 1 foo))))
(rplaca foo (- (nth 0 foo) (* q (nth 1 foo))))
(rplaca b (mapcar #'- (nth 0 b) (mapcar #'(lambda(x) (* x q)) (nth 1 b))))
) ;; report final list and integer relation
(list foo b))

; Try it... but use double-floats
(setf xx (list (cos 0.3d0) (cos 0.9d0) (cos 1.5d0) (expt (cos 0.3d0) 5)))
(Brun xx 9)
;; Reports ((2.1221132479043447d-7 3.9446742616799924d-7 6.493240857496918d-8
;;  5.412337245047638d-16)
;; ((-114 -12361 1324 9675) (23 2459 -264 -1925) (-73 -2833 311 2273)
;;  (10 5 1 -16)))
;;  Note 10 cos(x)+5 cos(3x)+cos(5x)-16 cos^5(x)=0

Brun’s algorithm is very sensitive to precision limitations, it may be easy to code (write it in your favorite language but use mpfr or Rmpfr) but the resulting program is not much use if the input precision is not at least long double.

You don’t like lisp? Here’s the code in julia. Note that 9 iterations gives the correct integer relation in lisp, but we need 30 in julia for the same input vector, because of lisp’s superior floating point handling.

brun=function(x,iters)
X=copy(x)
n=size(x)[1]
b =zeros(Float64,n,n)
for k in 1:n
b[k,k]=1.0
end

for j in 1:iters
p=sortperm(X,rev=true)
P=zeros(Float64,n,n)
for j in 1:n
P[j,p[j]]=1.0
end

X=P*X
b=P*b
r=fld(X[1],X[2])
X[1]=X[1] .- r*X[2]
b[1,:]=b[1,:] .- r*b[2,:]
end
b[n,:], sum(b[n,:] .*x)
end

v=[cos(3//10),cos(9//10),cos(15//10),cos(3//10)^5]
# needs 30 iterations
brun(v,30)
([10.0, 5.0, 1.0, -16.0], 1.7763568394002505e-15)

Home 2.0

Integer relations-I. Brun’s algorithm

Related Posts

Integer relations-I. Brun’s algorithm

Published by Jeff on August 2, 2021August 2, 2021

Related Posts

Linear equation solvers by LU decomposition

Algorithms: FFT-III. Julia FFTW and FITS file libraries

Algorithms: Tridiagonal Hessenberg form using Krylov subspaces