A sequence of much faster methods have however been developed over the past few decades, one simple example that works for most n being the so-called rho method of John Pollard (compare the quadratic residue sequences discussed below): Module[{f = Mod[# 2 + 1, n] &, a = 2, b = 5, c}, While[(c = GCD[n, a - b])  1, {a, b} = {f[a], f[f[b]]}]; c] Most existing methods depend on facts in number theory that are fairly easy to state, though implementing them for maximum efficiency tends to lead to complex programs.

The patterns are exactly repetitive only when Tan[ θ ]  u/v , where u and v are elements of a primitive Pythagorean triple (so that u , v and Sqrt[u 2 + v 2 ] are all integers, and GCD[u, v]  1 ).

[Intractability in] systems of limited size In the system x  Mod[x + m, n] from page 255 the repetition period n/GCD[m, n] can be computed using Euclid's algorithm in at most about Log[GoldenRatio, n] steps.

The most common are ones based on repetition or iteration, classic examples being Euclid's algorithm for GCD (page 915 ), Newton's method for FindRoot and the Gaussian elimination method for LinearSolve .

One based on gradual extension of work by Richard Stoneham from 1971 is that numbers of the form Sum[1/(p n b p n ), {n, ∞ }] for prime p > 2 are normal in base b (for GCD[b, p]  1 ), and are transcendental.

What is behind this is that each of the numbers in the basic sequence here must be a so-called quadratic residue of the form Mod[v 2 , m] , and given any such quadratic residue x the expression GCD[x + Mod[x 2 , m], m] turns out always to be a factor of m —and at least sometimes a non-trivial one.

Linear Diophantine equations such as a x  b y + c yield simple repetitive results, as in the pictures below, and can be handled essentially just by knowing ExtendedGCD[a, b] .

↔ a  Quotient[c b , Binomial[c, b]] a  GCD[b, c] ↔ (b c > 0 ∧ a d  b ∧ a e  c ∧ a + c f  b g) a  Floor[b/c] ↔ (a c + d  b ∧ d < c) PrimeQ[a] ↔ (GCD[(a - 1)!

As discovered by Srinivasa Ramanujan in 1918 its fluctuations (see below) can be obtained from the formula 1/6 π 2 n Sum[Apply[Plus, Cos[2 π n Select[ Range[s], GCD[s, #]  1 &]/s]]/s 2 , {s, ∞ }] (c) Squares are taken to be of positive or negative integers, or zero.

This same basic scheme can be used with any associative function h — Max , GCD , And , Dot , Join or whatever—so long as suitable forms for r and s are used.