<html><head><title>[ref] 15 Number Theory</title></head>
<body text="#000000" bgcolor="#ffffff">
[<a href="../index.htm">Top</a>] [<a href = "chapters.htm">Up</a>] [<a href ="CHAP014.htm">Previous</a>] [<a href ="CHAP016.htm">Next</a>] [<a href = "theindex.htm">Index</a>]
<h1>15 Number Theory</h1><p>
<P>
<H3>Sections</H3>
<oL>
<li> <A HREF="CHAP015.htm#SECT001">Prime Residues</a>
<li> <A HREF="CHAP015.htm#SECT002">Primitive Roots and Discrete Logarithms</a>
<li> <A HREF="CHAP015.htm#SECT003">Roots Modulo Integers</a>
<li> <A HREF="CHAP015.htm#SECT004">Multiplicative Arithmetic Functions</a>
<li> <A HREF="CHAP015.htm#SECT005">Continued Fractions</a>
<li> <A HREF="CHAP015.htm#SECT006">Miscellaneous</a>
</ol><p>
<p>
<a name = "I0"></a>
<font face="Gill Sans,Helvetica,Arial">GAP</font> provides a couple of elementary number theoretic functions.
Most of these deal with the group of integers coprime to <i>m</i>,
called the <strong>prime residue group</strong>.
φ(<i>m</i>) (see <a href="CHAP015.htm#SSEC001.2">Phi</a>) is the order of this group,
λ(<i>m</i>) (see <a href="CHAP015.htm#SSEC001.3">Lambda</a>) the exponent.
If and only if <i>m</i> is 2, 4, an odd prime power <i>p</i><sup><i>e</i></sup>,
or twice an odd prime power 2 <i>p</i><sup><i>e</i></sup>, this group is cyclic.
In this case the generators of the group, i.e., elements of order
φ(<i>m</i>), are called <strong>primitive roots</strong> (see <a href="CHAP015.htm#SSEC002.3">PrimitiveRootMod</a>,
<a href="CHAP015.htm#SSEC002.4">IsPrimitiveRootMod</a>).
<p>
Note that neither the arguments nor the return values of the functions
listed below are groups or group elements in the sense of <font face="Gill Sans,Helvetica,Arial">GAP</font>.
The arguments are simply integers.
<p>
<a name = ""></a>
<li><code>InfoNumtheor V</code>
<p>
<code>InfoNumtheor</code> is the info class (see <a href="CHAP007.htm#SECT004">Info Functions</a>) for the
functions in the number theory chapter.
<p>
<p>
<h2><a name="SECT001">15.1 Prime Residues</a></h2>
<p><p>
<a name = "SSEC001.1"></a>
<li><code>PrimeResidues( </code><var>m</var><code> ) F</code>
<p>
<code>PrimeResidues</code> returns the set of integers from the range <code>0..Abs(</code><var>m</var><code>)-1</code>
that are coprime to the integer <var>m</var>.
<p>
<code>Abs(</code><var>m</var><code>)</code> must be less than 2<sup>28</sup>, otherwise the set would probably be
too large anyhow.
<p>
<a name = "I1"></a>
<pre>
gap> PrimeResidues( 0 ); PrimeResidues( 1 ); PrimeResidues( 20 );
[ ]
[ 0 ]
[ 1, 3, 7, 9, 11, 13, 17, 19 ]
</pre>
<p>
<a name = "SSEC001.2"></a>
<li><code>Phi( </code><var>m</var><code> ) O</code>
<p>
<code>Phi</code> returns the number φ(<i>m</i> ) of positive integers less than the
positive integer <var>m</var> that are coprime to <var>m</var>.
<p>
Suppose that <i>m</i> = <i>p</i><sub>1</sub><sup><i>e</i><sub>1</sub></sup> <i>p</i><sub>2</sub><sup><i>e</i><sub>2</sub></sup> …<i>p</i><sub><i>k</i></sub><sup><i>e</i><sub><i>k</i></sub></sup>. Then φ(<i>m</i>)
is <i>p</i><sub>1</sub><sup><i>e</i><sub>1</sub>−1</sup> (<i>p</i><sub>1</sub>−1) <i>p</i><sub>2</sub><sup><i>e</i><sub>2</sub>−1</sup> (<i>p</i><sub>2</sub>−1) …<i>p</i><sub><i>k</i></sub><sup><i>e</i><sub><i>k</i></sub>−1</sup> (<i>p</i><sub><i>k</i></sub>−1).
<p>
<a name = "I2"></a>
<a name = "I3"></a>
<a name = "I4"></a>
<pre>
gap> Phi( 12 );
4
gap> Phi( 2^13-1 ); # this proves that 2^(13)-1 is a prime
8190
gap> Phi( 2^15-1 );
27000
</pre>
<p>
<a name = "SSEC001.3"></a>
<li><code>Lambda( </code><var>m</var><code> ) O</code>
<p>
<code>Lambda</code> returns the exponent λ(<i>m</i> ) of the group of prime
residues modulo the integer <var>m</var>.
<p>
λ(<i>m</i> ) is the smallest positive integer <i>l</i> such that for every
<i>a</i> relatively prime to <var>m</var> we have <i>a</i><sup><i>l</i></sup> ≡ 1 mod <i>m</i> .
Fermat's theorem asserts <i>a</i><sup>φ(<i>m</i> )</sup> ≡ 1 mod <i>m</i> ;
thus λ(<i>m</i>) divides φ(<i>m</i>) (see <a href="CHAP015.htm#SSEC001.2">Phi</a>).
<p>
Carmichael's theorem states that λ can be computed as follows:
λ(2) = 1, λ(4) = 2 and λ(2<sup><i>e</i></sup>) = 2<sup><i>e</i>−2</sup>
if 3 ≤ <i>e</i>,
λ(<i>p</i><sup><i>e</i></sup>) = (<i>p</i>−1) <i>p</i><sup><i>e</i>−1</sup> (i.e. φ(<i>m</i>)) if <i>p</i> is an odd prime
and
λ(<i>n</i>*<i>m</i>) = <tt>Lcm</tt>( λ(<i>n</i>), λ(<i>m</i>) ) if <i>n</i>, <i>m</i> are coprime.
<p>
Composites for which λ(<i>m</i>) divides <i>m</i> − 1 are called Carmichaels.
If 6<i>k</i>+1, 12<i>k</i>+1 and 18<i>k</i>+1 are primes their product is such a number.
There are only 1547 Carmichaels below 10<sup>10</sup> but 455052511 primes.
<p>
<a name = "I5"></a>
<a name = "I6"></a>
<a name = "I7"></a>
<pre>
gap> Lambda( 10 );
4
gap> Lambda( 30 );
4
gap> Lambda( 561 ); # 561 is the smallest Carmichael number
80
</pre>
<p>
<a name = "SSEC001.4"></a>
<li><code>GeneratorsPrimeResidues( </code><var>n</var><code> ) F</code>
<p>
Let <var>n</var> be a positive integer.
<code>GeneratorsPrimeResidues</code> returns a description of generators of the
group of prime residues modulo <var>n</var>.
The return value is a record with components
<p>
<dl compact>
<dt><code>primes</code>: <dd>
a list of the prime factors of <var>n</var>,
<p>
<dt><code>exponents</code>: <dd>
a list of the exponents of these primes in the factorization of <var>n</var>,
and
<p>
<dt><code>generators</code>: <dd>
a list describing generators of the group of prime residues;
for the prime factor 2, either a primitive root or a list of two
generators is stored,
for each other prime factor of <var>n</var>, a primitive root is stored.
</dl>
<p>
<pre>
gap> GeneratorsPrimeResidues( 1 );
rec( primes := [ ], exponents := [ ], generators := [ ] )
gap> GeneratorsPrimeResidues( 4*3 );
rec( primes := [ 2, 3 ], exponents := [ 2, 1 ], generators := [ 7, 5 ] )
gap> GeneratorsPrimeResidues( 8*9*5 );
rec( primes := [ 2, 3, 5 ], exponents := [ 3, 2, 1 ],
generators := [ [ 271, 181 ], 281, 217 ] )
</pre>
<p>
<p>
<h2><a name="SECT002">15.2 Primitive Roots and Discrete Logarithms</a></h2>
<p><p>
<a name = "SSEC002.1"></a>
<li><code>OrderMod( </code><var>n</var><code>, </code><var>m</var><code> ) F</code>
<p>
<code>OrderMod</code> returns the multiplicative order of the integer <var>n</var> modulo the
positive integer <var>m</var>.
If <var>n</var> and <var>m</var> are not coprime the order of <var>n</var> is not defined
and <code>OrderMod</code> will return 0.
<p>
If <var>n</var> and <var>m</var> are relatively prime the multiplicative order of <var>n</var>
modulo <var>m</var> is the smallest positive integer <i>i</i> such that
<i>n</i> <sup><i>i</i></sup> ≡ 1 mod <i>m</i> .
If the group of prime residues modulo <var>m</var> is cyclic then each element
of maximal order is called a primitive root modulo <var>m</var>
(see <a href="CHAP015.htm#SSEC002.4">IsPrimitiveRootMod</a>).
<p>
<code>OrderMod</code> usually spends most of its time factoring <var>m</var> and φ(<i>m</i> )
(see <a href="CHAP014.htm#SSEC003.6">FactorsInt</a>).
<p>
<a name = "I8"></a>
<pre>
gap> OrderMod( 2, 7 );
3
gap> OrderMod( 3, 7 ); # 3 is a primitive root modulo 7
6
</pre>
<p>
<a name = "SSEC002.2"></a>
<li><code>LogMod( </code><var>n</var><code>, </code><var>r</var><code>, </code><var>m</var><code> ) F</code>
<a name = "SSEC002.2"></a>
<li><code>LogModShanks( </code><var>n</var><code>, </code><var>r</var><code>, </code><var>m</var><code> ) F</code>
<p>
computes the discrete <var>r</var>-logarithm of the integer <var>n</var> modulo the integer
<var>m</var>.
It returns a number <var>l</var> such that <i>r</i> <sup><i>l</i> </sup> ≡ <i>n</i> mod <i>m</i>
if such a number exists.
Otherwise <code>fail</code> is returned.
<p>
<code>LogModShanks</code> uses the Baby Step - Giant Step Method of Shanks (see for
example section 5.4.1 in <a href="biblio.htm#Coh93"><cite>Coh93</cite></a> and in general
requires more memory than a call to <code>LogMod</code>.
<p>
<a name = "I9"></a>
<pre>
gap> l:= LogMod( 2, 5, 7 ); 5^l mod 7 = 2;
4
true
gap> LogMod( 1, 3, 3 ); LogMod( 2, 3, 3 );
0
fail
</pre>
<p>
<a name = "SSEC002.3"></a>
<li><code>PrimitiveRootMod( </code><var>m</var><code>[, </code><var>start</var><code>] ) F</code>
<p>
<code>PrimitiveRootMod</code> returns the smallest primitive root modulo the
positive integer <var>m</var> and <code>fail</code> if no such primitive root exists.
If the optional second integer argument <var>start</var> is given
<code>PrimitiveRootMod</code> returns the smallest primitive root that is strictly
larger than <var>start</var>.
<p>
<a name = "I10"></a>
<a name = "I11"></a>
<a name = "I12"></a>
<pre>
gap> PrimitiveRootMod( 409 ); # largest primitive root for a prime less than 2000
21
gap> PrimitiveRootMod( 541, 2 );
10
gap> PrimitiveRootMod( 337, 327 ); # 327 is the largest primitive root mod 337
fail
gap> PrimitiveRootMod( 30 ); # there exists no primitive root modulo 30
fail
</pre>
<p>
<a name = "SSEC002.4"></a>
<li><code>IsPrimitiveRootMod( </code><var>r</var><code>, </code><var>m</var><code> ) F</code>
<p>
<code>IsPrimitiveRootMod</code> returns <code>true</code> if the integer <var>r</var> is a primitive
root modulo the positive integer <var>m</var> and <code>false</code> otherwise. If <var>r</var> is
less than 0 or larger than <var>m</var> it is replaced by its remainder.
<p>
<a name = "I13"></a>
<pre>
gap> IsPrimitiveRootMod( 2, 541 );
true
gap> IsPrimitiveRootMod( -539, 541 ); # same computation as above;
true
gap> IsPrimitiveRootMod( 4, 541 );
false
gap> ForAny( [1..29], r -> IsPrimitiveRootMod( r, 30 ) );
false
gap> # there is no a primitive root modulo 30
</pre>
<p>
<p>
<h2><a name="SECT003">15.3 Roots Modulo Integers</a></h2>
<p><p>
<a name = "SSEC003.1"></a>
<li><code>Jacobi( </code><var>n</var><code>, </code><var>m</var><code> ) F</code>
<p>
<code>Jacobi</code> returns the value of the <strong>Jacobi symbol</strong> of the integer <var>n</var>
modulo the integer <var>m</var>.
<p>
Suppose that <i>m</i> = <i>p</i><sub>1</sub> <i>p</i><sub>2</sub> …<i>p</i><sub><i>k</i></sub> is a product of primes, not
necessarily distinct. Then for <var>n</var> coprime to <var>m</var> the Jacobi symbol is
defined by <i>J</i>(<i>n</i> /<i>m</i> ) = <i>L</i>(<i>n</i> /<i>p</i><sub>1</sub>) <i>L</i>(<i>n</i> /<i>p</i><sub>2</sub>) …<i>L</i>(<i>n</i> /<i>p</i><sub><i>k</i></sub>), where
<i>L</i>(<i>n</i> /<i>p</i>) is the Legendre symbol (see <a href="CHAP015.htm#SSEC003.2">Legendre</a>). By convention
<i>J</i>(<i>n</i> /1) = 1. If the gcd of <var>n</var> and <var>m</var> is larger than 1 we define
<i>J</i>(<i>n</i> /<i>m</i> ) = 0.
<p>
If <var>n</var> is a <strong>quadratic residue</strong> modulo <var>m</var>, i.e., if there exists an <i>r</i>
such that <i>r</i><sup>2</sup> ≡ <i>n</i> mod <i>m</i> then <i>J</i>(<i>n</i> /<i>m</i> ) = 1. However,
<i>J</i>(<i>n</i> /<i>m</i> ) = 1 implies the existence of such an <i>r</i> only if <var>m</var> is a
prime.
<p>
<code>Jacobi</code> is very efficient, even for large values of <var>n</var> and <var>m</var>, it is
about as fast as the Euclidean algorithm (see <a href="CHAP054.htm#SSEC007.1">Gcd</a>).
<p>
<a name = "I14"></a>
<a name = "I15"></a>
<pre>
gap> Jacobi( 11, 35 ); # 9^2 = 11 mod 35
1
gap> Jacobi( 6, 35 ); # it is -1, thus there is no r such that r^2 = 6 mod 35
-1
gap> Jacobi( 3, 35 ); # it is 1 even though there is no r with r^2 = 3 mod 35
1
</pre>
<p>
<a name = "SSEC003.2"></a>
<li><code>Legendre( </code><var>n</var><code>, </code><var>m</var><code> ) F</code>
<p>
<code>Legendre</code> returns the value of the <strong>Legendre symbol</strong> of the integer <var>n</var>
modulo the positive integer <var>m</var>.
<p>
The value of the Legendre symbol <i>L</i>(<i>n</i> /<i>m</i> ) is 1 if <var>n</var> is a
<strong>quadratic residue</strong> modulo <var>m</var>, i.e., if there exists an integer <i>r</i> such
that <i>r</i><sup>2</sup> ≡ <i>n</i> mod <i>m</i> and −1 otherwise.
<p>
If a root of <var>n</var> exists it can be found by <code>RootMod</code> (see <a href="CHAP015.htm#SSEC003.3">RootMod</a>).
<p>
While the value of the Legendre symbol usually is only defined for <var>m</var> a
prime, we have extended the definition to include composite moduli too.
The Jacobi symbol (see <a href="CHAP015.htm#SSEC003.1">Jacobi</a>) is another generalization of the
Legendre symbol for composite moduli that is much cheaper to compute,
because it does not need the factorization of <var>m</var> (see <a href="CHAP014.htm#SSEC003.6">FactorsInt</a>).
<p>
A description of the Jacobi symbol, the Legendre symbol, and related
topics can be found in <a href="biblio.htm#Baker84"><cite>Baker84</cite></a>.
<p>
<pre>
gap> Legendre( 5, 11 ); # 4^2 = 5 mod 11
1
gap> Legendre( 6, 11 ); # it is -1, thus there is no r such that r^2 = 6 mod 11
-1
gap> Legendre( 3, 35 ); # it is -1, thus there is no r such that r^2 = 3 mod 35
-1
</pre>
<p>
<a name = "SSEC003.3"></a>
<li><code>RootMod( </code><var>n</var><code>[, </code><var>k</var><code>], </code><var>m</var><code> ) F</code>
<p>
<code>RootMod</code> computes a <var>k</var>th root of the integer <var>n</var> modulo the positive
integer <var>m</var>, i.e., a <i>r</i> such that <i>r</i><sup><i>k</i> </sup> ≡ <i>n</i> mod <i>m</i> .
If no such root exists <code>RootMod</code> returns <code>fail</code>.
If only the arguments <var>n</var> and <var>m</var> are given,
the default value for <var>k</var> is 2.
<p>
In the current implementation <var>k</var> must be a prime.
<p>
A square root of <var>n</var> exists only if <code>Legendre(</code><var>n</var><code>,</code><var>m</var><code>) = 1</code>
(see <a href="CHAP015.htm#SSEC003.2">Legendre</a>).
If <var>m</var> has <i>r</i> different prime factors then there are 2<sup><i>r</i></sup> different
roots of <var>n</var> mod <var>m</var>. It is unspecified which one <code>RootMod</code> returns.
You can, however, use <code>RootsMod</code> (see <a href="CHAP015.htm#SSEC003.4">RootsMod</a>) to compute the full set
of roots.
<p>
<code>RootMod</code> is efficient even for large values of <var>m</var>, in fact the most
time is usually spent factoring <var>m</var> (see <a href="CHAP014.htm#SSEC003.6">FactorsInt</a>).
<p>
<a name = "I16"></a>
<pre>
gap> RootMod( 64, 1009 ); # note 'RootMod' does not return 8 in this case but -8;
1001
gap> RootMod( 64, 3, 1009 );
518
gap> RootMod( 64, 5, 1009 );
656
gap> List( RootMod( 64, 1009 ) * RootsUnityMod( 1009 ),
> x -> x mod 1009 ); # set of all square roots of 64 mod 1009
[ 1001, 8 ]
</pre>
<p>
<a name = "SSEC003.4"></a>
<li><code>RootsMod( </code><var>n</var><code>[, </code><var>k</var><code>], </code><var>m</var><code> ) F</code>
<p>
<code>RootsMod</code> computes the set of <var>k</var>th roots of the integer <var>n</var>
modulo the positive integer <var>m</var>, i.e., a <i>r</i> such that
<i>r</i><sup><i>k</i> </sup> ≡ <i>n</i> mod <i>m</i> .
If only the arguments <var>n</var> and <var>m</var> are given,
the default value for <var>k</var> is 2.
<p>
In the current implementation <var>k</var> must be a prime.
<p>
<pre>
gap> RootsMod( 1, 7*31 ); # the same as `RootsUnityMod( 7*31 )'
[ 1, 92, 125, 216 ]
gap> RootsMod( 7, 7*31 );
[ 21, 196 ]
gap> RootsMod( 5, 7*31 );
[ ]
gap> RootsMod( 1, 5, 7*31 );
[ 1, 8, 64, 78, 190 ]
</pre>
<p>
<a name = "SSEC003.5"></a>
<li><code>RootsUnityMod( [</code><var>k</var><code>, ] </code><var>m</var><code> ) F</code>
<p>
<code>RootsUnityMod</code> returns the set of <var>k</var>-th roots of unity modulo the
positive integer <var>m</var>, i.e.,
the list of all solutions <i>r</i> of <i>r</i><sup><i>k</i> </sup> ≡ <i>n</i> mod <i>m</i> .
If only the argument <var>m</var> is given, the default value for <var>k</var> is 2.
<p>
In general there are <i>k</i> <sup><i>n</i></sup> such roots if the modulus <var>m</var> has <i>n</i>
different prime factors <i>p</i> such that <i>p</i> ≡ 1 mod <i>k</i> . If <i>k</i> <sup>2</sup>
divides <var>m</var> then there are <i>k</i> <sup><i>n</i>+1</sup> such roots; and especially if
<i>k</i> = 2 and 8 divides <var>m</var> there are 2<sup><i>n</i>+2</sup> such roots.
<p>
In the current implementation <var>k</var> must be a prime.
<p>
<a name = "I17"></a>
<a name = "I18"></a>
<pre>
gap> RootsUnityMod( 7*31 ); RootsUnityMod( 3, 7*31 );
[ 1, 92, 125, 216 ]
[ 1, 25, 32, 36, 67, 149, 156, 191, 211 ]
gap> RootsUnityMod( 5, 7*31 );
[ 1, 8, 64, 78, 190 ]
gap> List( RootMod( 64, 1009 ) * RootsUnityMod( 1009 ),
> x -> x mod 1009 ); # set of all square roots of 64 mod 1009
[ 1001, 8 ]
</pre>
<p>
<p>
<h2><a name="SECT004">15.4 Multiplicative Arithmetic Functions</a></h2>
<p><p>
<a name = "SSEC004.1"></a>
<li><code>Sigma( </code><var>n</var><code> ) O</code>
<p>
<code>Sigma</code> returns the sum of the positive divisors of the nonzero integer
<var>n</var>.
<p>
<code>Sigma</code> is a multiplicative arithmetic function, i.e., if <var>n</var> and <i>m</i> are
relatively prime we have σ(<i>n</i> <i>m</i>) = σ(<i>n</i> ) σ(<i>m</i>).
<p>
Together with the formula σ(<i>p</i><sup><i>e</i></sup>) = (<i>p</i><sup><i>e</i>+1</sup>−1) / (<i>p</i>−1) this allows
us to compute σ(<i>n</i> ).
<p>
Integers <var>n</var> for which σ(<i>n</i> )=2 <i>n</i> are called perfect. Even
perfect integers are exactly of the form 2<sup><i>n</i> −1</sup>(2<sup><i>n</i> </sup>−1) where
2<sup><i>n</i> </sup>−1 is prime. Primes of the form 2<sup><i>n</i> </sup>−1 are called <strong>Mersenne
primes</strong>, the known ones are obtained for <var>n</var> = 2, 3, 5, 7, 13, 17, 19,
31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689,
9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091,
756839, and 859433. It is not known whether odd perfect integers exist,
however <a href="biblio.htm#BC89"><cite>BC89</cite></a> show that any such integer must have at least 300
decimal digits.
<p>
<code>Sigma</code> usually spends most of its time factoring <var>n</var> (see <a href="CHAP014.htm#SSEC003.6">FactorsInt</a>).
<p>
<pre>
gap> Sigma( 1 );
1
gap> Sigma( 1009 ); # 1009 is a prime
1010
gap> Sigma( 8128 ) = 2*8128; # 8128 is a perfect number
true
</pre>
<p>
<a name = "SSEC004.2"></a>
<li><code>Tau( </code><var>n</var><code> ) O</code>
<p>
<code>Tau</code> returns the number of the positive divisors of the nonzero integer
<var>n</var>.
<p>
<code>Tau</code> is a multiplicative arithmetic function, i.e., if <var>n</var> and <i>m</i> are
relative prime we have τ(<i>n</i> <i>m</i>) = τ(<i>n</i> ) τ(<i>m</i>).
Together with the formula τ(<i>p</i><sup><i>e</i></sup>) = <i>e</i>+1 this allows us to compute
τ(<i>n</i> ).
<p>
<code>Tau</code> usually spends most of its time factoring <var>n</var> (see <a href="CHAP014.htm#SSEC003.6">FactorsInt</a>).
<p>
<pre>
gap> Tau( 1 );
1
gap> Tau( 1013 ); # thus 1013 is a prime
2
gap> Tau( 8128 );
14
gap> Tau( 36 ); # the result is odd if and only if the argument is a perfect square
9
</pre>
<p>
<a name = "SSEC004.3"></a>
<li><code>MoebiusMu( </code><var>n</var><code> ) F</code>
<p>
<code>MoebiusMu</code> computes the value of Moebius inversion function for the
nonzero integer <var>n</var>.
This is 0 for integers which are not squarefree, i.e.,
which are divided by a square <i>r</i><sup>2</sup>. Otherwise it is 1 if <var>n</var> has a even
number and −1 if <var>n</var> has an odd number of prime factors.
<p>
The importance of μ stems from the so called inversion formula.
Suppose <i>f</i>(<i>n</i> ) is a multiplicative arithmetic function defined on the
positive integers and let <i>g</i>(<i>n</i> )=∑<sub><i>d</i> | <i>n</i> </sub><i>f</i>(<i>d</i>). Then
<i>f</i>(<i>n</i> )=∑<sub><i>d</i> | <i>n</i> </sub>μ(<i>d</i>) <i>g</i>(<i>n</i> /<i>d</i>). As a special case we have
φ(<i>n</i> ) = ∑<sub><i>d</i> | <i>n</i> </sub>μ(<i>d</i>) <i>n</i> /<i>d</i> since
<i>n</i> = ∑<sub><i>d</i> | <i>n</i> </sub>φ(<i>d</i>) (see <a href="CHAP015.htm#SSEC001.2">Phi</a>).
<p>
<code>MoebiusMu</code> usually spends all of its time factoring <var>n</var> (see
<a href="CHAP014.htm#SSEC003.6">FactorsInt</a>).
<p>
<pre>
gap> MoebiusMu( 60 ); MoebiusMu( 61 ); MoebiusMu( 62 );
0
-1
1
</pre>
<p>
<p>
<h2><a name="SECT005">15.5 Continued Fractions</a></h2>
<p><p>
<a name = "SSEC005.1"></a>
<li><code>ContinuedFractionExpansionOfRoot( </code><var>P</var><code>, </code><var>n</var><code> ) F</code>
<p>
The first <var>n</var> terms of the continued fraction expansion of the only
positive real root of the polynomial <var>P</var> with integer coefficients.
The leading coefficient of <var>P</var> must be positive and the value of <var>P</var> at 0
must be negative. If the degree of <var>P</var> is 2 and <var>n</var> = 0, the function
computes one period of the continued fraction expansion of the root in
question. Anything may happen if <var>P</var> has three or more positive real
roots.
<p>
<pre>
gap> x := Indeterminate(Integers);;
gap> ContinuedFractionExpansionOfRoot(x^2-7,20);
[ 2, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1 ]
gap> ContinuedFractionExpansionOfRoot(x^2-7,0);
[ 2, 1, 1, 1, 4 ]
gap> ContinuedFractionExpansionOfRoot(x^3-2,20);
[ 1, 3, 1, 5, 1, 1, 4, 1, 1, 8, 1, 14, 1, 10, 2, 1, 4, 12, 2, 3 ]
gap> ContinuedFractionExpansionOfRoot(x^5-x-1,50);
[ 1, 5, 1, 42, 1, 3, 24, 2, 2, 1, 16, 1, 11, 1, 1, 2, 31, 1, 12, 5, 1, 7, 11,
1, 4, 1, 4, 2, 2, 3, 4, 2, 1, 1, 11, 1, 41, 12, 1, 8, 1, 1, 1, 1, 1, 9, 2,
1, 5, 4 ]
</pre>
<p>
<a name = "SSEC005.2"></a>
<li><code>ContinuedFractionApproximationOfRoot( </code><var>P</var><code>, </code><var>n</var><code> ) F</code>
<p>
The <var>n</var>th continued fraction approximation of the only positive real root
of the polynomial <var>P</var> with integer coefficients. The leading coefficient
of <var>P</var> must be positive and the value of <var>P</var> at 0 must be negative.
Anything may happen if <var>P</var> has three or more positive real roots.
<p>
<pre>
gap> ContinuedFractionApproximationOfRoot(x^2-2,10);
3363/2378
gap> 3363^2-2*2378^2;
1
gap> z := ContinuedFractionApproximationOfRoot(x^5-x-1,20);
499898783527/428250732317
gap> z^5-z-1;
486192462527432755459620441970617283/
14404247382319842421697357558805709031116987826242631261357
</pre>
<p>
<p>
<h2><a name="SECT006">15.6 Miscellaneous</a></h2>
<p><p>
<a name = "SSEC006.1"></a>
<li><code>TwoSquares( </code><var>n</var><code> ) F</code>
<p>
<code>TwoSquares</code> returns a list of two integers <i>x</i> ≤ <i>y</i> such that the sum of
the squares of <i>x</i> and <i>y</i> is equal to the nonnegative integer <var>n</var>, i.e.,
<i>n</i> = <i>x</i><sup>2</sup>+<i>y</i><sup>2</sup>. If no such representation exists <code>TwoSquares</code> will return
<code>fail</code>. <code>TwoSquares</code> will return a representation for which the gcd of
<i>x</i> and <i>y</i> is as small as possible. It is not specified which
representation <code>TwoSquares</code> returns, if there is more than one.
<p>
Let <i>a</i> be the product of all maximal powers of primes of the form 4<i>k</i>+3
dividing <var>n</var>. A representation of <var>n</var> as a sum of two squares exists if
and only if <i>a</i> is a perfect square. Let <i>b</i> be the maximal power of 2
dividing <var>n</var> or its half, whichever is a perfect square. Then the minimal
possible gcd of <i>x</i> and <i>y</i> is the square root <i>c</i> of <i>a</i> <i>b</i>. The number
of different minimal representation with <i>x</i> ≤ <i>y</i> is 2<sup><i>l</i>−1</sup>, where <i>l</i>
is the number of different prime factors of the form 4<i>k</i>+1 of <var>n</var>.
<p>
The algorithm first finds a square root <i>r</i> of −1 modulo <i>n</i> / (<i>a</i> <i>b</i>),
which must exist, and applies the Euclidean algorithm to <i>r</i> and <var>n</var>. The
first residues in the sequence that are smaller than √{<i>n</i> /(<i>a</i> <i>b</i>)}
times <i>c</i> are a possible pair <i>x</i> and <i>y</i>.
<p>
Better descriptions of the algorithm and related topics can be found in
<a href="biblio.htm#Wagon90"><cite>Wagon90</cite></a> and <a href="biblio.htm#Zagier90"><cite>Zagier90</cite></a>.
<p>
<a name = "I19"></a>
<pre>
gap> TwoSquares( 5 );
[ 1, 2 ]
gap> TwoSquares( 11 ); # there is no representation
fail
gap> TwoSquares( 16 );
[ 0, 4 ]
gap> TwoSquares( 45 ); # 3 is the minimal possible gcd because 9 divides 45
[ 3, 6 ]
gap> TwoSquares( 125 ); # it is not [5,10] because their gcd is not minimal
[ 2, 11 ]
gap> TwoSquares( 13*17 ); # [10,11] would be the other possible representation
[ 5, 14 ]
gap> TwoSquares( 848654483879497562821 ); # 848654483879497562821 is prime
#I IsPrimeInt: probably prime, but not proven: 848654483879497562821
#I FactorsInt: used the following factor(s) which are probably primes:
#I 848654483879497562821
[ 6305894639, 28440994650 ]
</pre>
<p>
<p>
[<a href="../index.htm">Top</a>] [<a href = "chapters.htm">Up</a>] [<a href ="CHAP014.htm">Previous</a>] [<a href ="CHAP016.htm">Next</a>] [<a href = "theindex.htm">Index</a>]
<P>
<font face="Gill Sans,Helvetica,Arial">GAP 4 manual<br>December 2008
</font></body></html>
