**Number Theory Packages**

Functions relevant to number theory are well represented in Mathematica with examples such as PrimePi, EulerPhi, MoebiusMu, and DivisorSigma. The NumberTheory packages broaden this set of functions. There are packages for proving primality, exploring the elliptic curve method for integer factorization, and finding primitive elements of multiple algebraic extensions of rationals. There are functions for approximating real numbers by rationals and approximating polynomials with real roots by polynomials with integer coefficients. You can find the continued fraction expansion of a real number or the arbitrary base expansion of a rational in terms of preperiodic and periodic parts. Number theoretic functions such as Ramujan and Siegel

are also supported.

This causes the NumberTheory packages to load as they are needed.
In[1]:= **<<NumberTheory`**

Some of the plots in this section require the Graphics packages.
In[2]:= **<<Graphics`**

This is the continued fraction representation for to order

.
In[3]:= **c = ContinuedFraction[1/Pi, 7]**

Out[3]=

You can use the representation to find a gear ratio close to

such that there are fewer than 110 teeth on the smaller gear.
In[4]:= **(While[Numerator[Normal[c]] > 110,**

c = ContinuedFractionForm[Drop[c[[1]], -1]] ];

Normal[c])

Out[4]=

The ratios of the integers given by ProjectiveRationalize approximate the ratios of the corresponding real numbers , , , and within an error of

.
In[5]:= **(reals = N[{1, 1/EulerGamma, E, Pi}];**

prec = 5;

p = ProjectiveRationalize[reals, prec])

Out[5]=

Several gear ratios can be approximated simultaneously, under the constraint that the smallest gear have fewer than 110 teeth. Here the approximation to is

.
In[6]:= **(While[p[[1]] > 110,**

prec--;

p = ProjectiveRationalize[reals, prec]];

p)

Out[6]=

Square-free integers have no factors that are perfect squares. SquareFreeQ determines whether an integer is square-free.

Two integers are coprime if their greatest common divisor (GCD) is

. This shows that there is a positive correlation between the square-free property and coprimality.
In[7]:= **(sample = Table[Random[Integer, 1000000], {100}];**

xor = Map[(squarefree = SquareFreeQ[#];

coprime = (GCD[#, Random[Integer, 1000000]] == 1);

If[Xor[squarefree, coprime], -1, 1])&, sample];

N[Apply[Plus, xor]/Length[xor]])

Out[7]=

Every positive integer can be represented as the sum of four squares, but there are relatively few perfect squares. Each of the remaining integers can be minimally represented as a sum of either two or three squares. The number of representations of an integer as a sum of squares is given by SumOfSquaresR[

d,n].

This shows the percentage of numbers from to representable as a sum of a minimum of , , , or

squares.
In[8]:= **(squares = Table[**

Which[

SumOfSquaresR[1, n] != 0, 1,

SumOfSquaresR[2, n] != 0, 2,

SumOfSquaresR[3, n] != 0, 3,

True, 4], {n, 1, 100}];

PieChart[Map[Count[squares, #]&, Range[4]]])

PrimitiveRoot[p] gives a cyclic generator of the multiplicative group mod , where is the prime

. Primitive root arrays are used for wave scattering over a range of frequencies.

This gives the primitive root array for

.
In[9]:= **(p = 19;**

r = PrimitiveRoot[p];

array = Table[PowerMod[r, n, p], {n, 0, p-2}])

Out[9]=

Here is a depiction of a reflection phase-grating based on the array. The surface is pitted with a sequence of "wells" (shown in white) having depths proportional to the elements of the array.
In[10]:= **(array5 = Flatten[Table[array, {5}]];**

StackedBarChart[(p-array5), array5,

BarLabels -> None, BarSpacing -> 0,

BarStyle -> {GrayLevel[0], GrayLevel[1]},

AspectRatio ->.1, Ticks -> None])

Using a simple model of the complex amplitude of the reflected wave, it can be shown that the spectrum is flat except at zero. This implies that the wave is reflected in all directions except the specular direction.

Equal intensity is scattered into all diffraction orders except the zero order.
In[11]:= **(reflected = Map[Exp[I 2 Pi #/p]&, array];**

powerspectrum = Abs[Fourier[reflected]]^2;

ListPlot[powerspectrum, AxesOrigin -> {0, 0},

PlotJoined -> True, PlotRange -> All])