Application of Maxima:

Elementary number theory using Maxima

Elementary number theory using Maxima

Prime numbers

You might remember that for any integer  greater than 1,  is a  if its factors are 1 and itself. The integers 2, 3, 5, and 7 are primes, but 9 is not prime because . The command primep() is useful for testing whether or not an integer is prime:

(%i1) primep(2);
(%o1)                                true
(%i2) primep(3);
(%o2)                                true
(%i3) primep(5);
(%o3)                                true
(%i4) primep(7);
(%o4)                                true
(%i5) primep(9);
(%o5)                                false

And the command next_prime(n) returns the next prime number greater than or equal to :

(%i6) next_prime(9);
(%o6)                                 11
(%i7) next_prime(11);
(%o7)                                 13
(%i8) next_prime(13);
(%o8)                                 17
(%i9) next_prime(17);
(%o9)                                 19
(%i10) next_prime(19);
(%o10)                                23

Let’s now define a function called primes_first_n() in  to return a list of the first  primes, where  is a positive integer. Programming in the Maxima language is different from programming in other languages like ,, and . For example, if your variable is be assigned a number, you don’t need to define whether your variable is of type int, long, double, or bool. All you have to do is use the colon operator “:” to assign some value to a variable, like in the following example.

(%i50) num : 13$
(%i51) num;
(%o51)                                13
(%i52) str : "my string"$
(%i53) str;
(%o53)                             my string
(%i54) L : ["a", "list", "of", "strings"]$
(%i55) L;
(%o55)                      [a, list, of, strings]

Before defining the function primes_first_n(), there are two useful built-in Maxima functions that you should know about. These are append() and last(). The function append() can be used to append an element to a list, whereaslast() can be used to return the last element of a list:

(%i56) L : ["a", "list"];
(%o56)                             [a, list]
(%i57) L : append(L, ["of", "strings"]);
(%o57)                      [a, list, of, strings]
(%i58) L;
(%o58)                      [a, list, of, strings]
(%i59) last(L);
(%o59)                              strings

Below is the function primes_first_n() which pulls together the features of next_prime(n), append(), and last(). Notice that it is defined at the Maxima command line interface.

(%i60) primes_first_n(n) := (
           if n < 1 then
               []
           else (
               L : [2],
               for i:2 thru n do (
                   L : append(L, [next_prime(last(L))])
               ),
               L
           )
       )$

You can also put the above function inside a text file called, say, /home/mvngu/primes.mac with the following content:

/* Return the first n prime numbers.
 *
 * INPUT:
 *
 * - n -- a positive integer greater than 1.
 *
 * OUTPUT:
 *
 * - A list of the first n prime numbers. If n is less than 1, then return
 *   an empty list.
 */
primes_first_n(n) := (
    if n < 1 then
        []
    else (
        L : [2],
        for i:2 thru n do (
            L : append(L, [next_prime(last(L))])
        ),
        L
    )
)$

Like C++ and Java, Maxima also uses “/*” to denote the beginning of a comment block and “*/” to denote the end of a comment block. To load the content of the file /home/mvngu/primes.mac into Maxima, you use the command load(). Let’s load the above file and experiment with the custom-defined function primes_first_n():

(%i64) load("/home/mvngu/primes.mac");
(%o64)                      /home/mvngu/primes.mac
(%i65) primes_first_n(0);
(%o65)                                []
(%i66) primes_first_n(-1);
(%o66)                                []
(%i67) primes_first_n(-2);
(%o67)                                []
(%i68) primes_first_n(1);
(%o68)                                [2]
(%i69) primes_first_n(2);
(%o69)                              [2, 3]
(%i70) primes_first_n(3);
(%o70)                             [2, 3, 5]
(%i71) primes_first_n(10);
(%o71)               [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

Factorizing integers

Integer  is about breaking up an integer into smaller components. In , these smaller components are usually the prime factors of the integer. Use the command ifactors() to compute the prime factors of positive integers:

(%i76) ifactors(10);
(%o76)                         [[2, 1], [5, 1]]
(%i77) (2^1) * (5^1);
(%o77)                                10
(%i78) ifactors(25);
(%o78)                             [[5, 2]]
(%i79) 5^2;
(%o79)                                25
(%i80) ifactors(62);
(%o80)                         [[2, 1], [31, 1]]
(%i81) ifactors(72);
(%o81)                         [[2, 3], [3, 2]]
(%i82) (2^3) * (3^2);
(%o82)                                72

The prime factors of 10 are  and . When you multiply these two prime factors together, you end up with . The expression  is called the prime factorization of 10. Similarly, the expression  is the prime factorization of 25, and  is the prime factorization of 72.

Greatest common divisors

Closely related to integer factorization is the concept of  (GCD). The Maxima commandgcd() is able to compute the GCD of two expressions  and  where this makes sense. These two expressions may be integers, polynomials, or some objects for which it makes sense to compute their GCD. For the moment, let’s just work with integers:

(%i1) gcd(9, 12);
(%o1)                                  3
(%i2) gcd(21, 49);
(%o2)                                  7
(%i3) gcd(22, 11);
(%o3)                                 11

The GCD of two integers  and  can be recursively defined as follows:

where  is the remainder when  is divided by . The above recursive definition can be easily translated to a Maxima function for integer GCD as follows (credit goes to  for the Maxima code):

/* Return the greatest common divisor (GCD) of two integers a and b.
 *
 * INPUT:
 *
 * - a -- an integer
 * - b -- an integer
 *
 * OUTPUT:
 *
 * - The greatest common divisor of a and b.
 */
igcd(a, b) := block(
    if b = 0 then
        return(a)
    else
        return( igcd(b, mod(a,b)) )
);

Save the above code to a text file and load it first before using the function. Or you can define the function from the Maxima command line interface and proceed to use it:

(%i5) igcd(a, b) := block(
          if b = 0 then
              return(a)
          else
              return( igcd(b, mod(a,b)) )
      )$
(%i6) igcd(9, 12);
(%o6)                                  3
(%i7) igcd(21, 49);
(%o7)                                  7
(%i8) igcd(22, 11);
(%o8)                                 11

The  provides an interesting relationship between , and the pair  and . Here is a Maxima function definition courtesy of amca:

/* Apply the extended Euclidean algorithm to compute integers s and t such
 * that gcd(a,b) = as + bt.
 *
 * INPUT:
 *
 * - a -- an integer
 * - b -- an integer
 *
 * OUTPUT:
 *
 * - A triple of integers (s, t, d) satisfying the relationship 
 *   d = gcd(a,b) = as + bt. This algorithm does not guarantee that s and t
 *   are unique. There may be other pairs of s and t that satisfy the requirement.
 */
igcdex(a,b) := block(
    [d, x, y, d1, x1, y1],
    if b = 0 then
        return([1, 0, a])
    else (
        [x1, y1, d1] : igcdex(b, mod(a,b)),
        [x, y, d] : [y1, x1 - quotient(a,b)*y1, d1],
        return([x, y, d])
    )
);

Or you can define it from the Maxima command line:

(%i9) igcdex(a,b) := block(
          [d, x, y, d1, x1, y1],
          if b = 0 then
              return([1, 0, a])
          else (
              [x1, y1, d1] : igcdex(b, mod(a,b)),
              [x, y, d] : [y1, x1 - quotient(a,b)*y1, d1],
              return([x, y, d])
          )
      )$

Let’s use the function igcdex() for various pairs of integers and verify the result.

(%i15) igcdex(120, 23);
(%o15)                           [- 9, 47, 1]
(%i16) 120*(-9) + 23*47;
(%o16)                                 1
(%i17) igcdex(2000, 2009);
(%o17)                          [- 893, 889, 1]
(%i18) 2000*(-893) + 2009*889;
(%o18)                                 1
(%i19) igcdex(24, 56);
(%o19)                            [- 2, 1, 8]
(%i20) 24*(-2) + 56*1;
(%o20)                                 8

			
			
			

				Sage 3.3 released | mvngu
			

			

				
				

					Technical preview of David Harvey’s zn_poly library exposed to the Sage library (Martin Albrecht).
def f(p,n):
    P = PolynomialRing(GF(next_prime(p)),'x')
    f = P.random_element(n)
    g = P.random_element(n)

    t0 = cputime()
    r0 = f*g
    t0 = cputime(t0)

    t1 = cputime()
    r1 = f._mul_zn_poly(g)
    t1 = cputime(t1)

    assert(r0 == r1)

    return p,n,t0,t1

for i in range(21): 
   f(2**47,2**i)

					

						returns on sage.math
					
# (140737488355328, 1, 0.0, 0.0)
# (140737488355328, 2, 0.0, 0.0)
# (140737488355328, 4, 0.00099999999999766942, 0.0)
# (140737488355328, 8, 0.0, 0.0)
# (140737488355328, 16, 0.0, 0.0)
# (140737488355328, 32, 0.0059990000000027521, 0.0)
# (140737488355328, 64, 0.0, 0.0)
# (140737488355328, 128, 0.0, 0.0)
# (140737488355328, 256, 0.0, 0.0)
# (140737488355328, 512, 0.0, 0.00099999999999766942)
# (140737488355328, 1024, 0.00099999999999766942, 0.0)
# (140737488355328, 2048, 0.0020000000000024443, 0.0019989999999978636)
# (140737488355328, 4096, 0.0049989999999979773, 0.005000000000002558)
# (140737488355328, 8192, 0.010998000000000729, 0.011997999999998399)
# (140737488355328, 16384, 0.023995999999996798, 0.023997000000001378)
# (140737488355328, 32768, 0.050992000000000814, 0.052991999999996153)
# (140737488355328, 65536, 0.1149820000000048, 0.10598499999999689)
# (140737488355328, 131072, 0.29195599999999189, 0.21996599999999944)
# (140737488355328, 262144, 0.6119070000000022, 0.45393199999999467)
# (140737488355328, 524288, 1.5217689999999919, 1.0278430000000043)
# (140737488355328, 1048576, 3.1365230000000111, 2.0966819999999871)

				
				

					Use Pohlig-Hellman for generic discrete logarithm (Yann Laigle-Chapuy) — This results in significant improvement in performance and less memory foot print. Here’s an example with a smooth order:
sage: factor(5^15-1)
2^2 * 11 * 31 * 71 * 181 * 1741

# BEFORE
sage: F.=GF(5^15)
sage: g=F.gen()
sage: u=g^123456789
sage: time log(u,g)
CPU times: user 271.39 s, sys: 4.72 s, total: 276.11 s
Wall time: 276.96 s
123456789
sage: get_memory_usage()
378.21875

# AFTER
sage: F.=GF(5^15)
sage: g=F.gen()
sage: u=g^123456789
sage: time log(u,g)
CPU times: user 0.14 s, sys: 0.00 s, total: 0.14 s
Wall time: 0.16 s
123456789
sage: get_memory_usage()
115.8984375

					

						And here’s another example with a not-so-smooth order:
					
sage:factor(3^13-1)
2 * 797161

# BEFORE
sage: F.=GF(3**13)
sage: g=F.gen()
sage: u=g^1234567
sage: timeit('log(u,g)')
5 loops, best of 3: 1.54 s per loop
sage: get_memory_usage()
155.11328125

# AFTER
sage: F.=GF(3**13)
sage: g=F.gen()
sage: u=g^1234567
sage: timeit('log(u,g)')
5 loops, best of 3: 931 ms per loop
sage: get_memory_usage()
139.4296875

				

			


			

				 
			
			

			
			

			
			

				

					AfterMath issue 5 is out
				
				

					

						

							It’s now the end of Semester 2 of 2008, all examinations are over and done with, and issue 5 of  is out for your reading pleasure. You can grab the , hot off the digital press.
						
						

							This issue contains articles covering topics such as:
						
						

							

								the Hairy Ball Theorem
							
							

								a differential equation variant of Gronwall’s inequality
							
							

								a survey on Maxima
							
							

								using PARI/GP to study number theory and cryptography
							
						
						

							Mark Ioppolo’s article, “The Hairy Ball Theorem,” contains a proof of a version of Henri Poincare’s conjecture popularly known as the Hairy Ball Theorem. Poincare was able to provide a proof [5] of the version he conjectured, and Brouwer [1] provides a proof of a generalization of Poincare’s version. The proofs in [1, 5] rely upon rather advanced techniques, whereas Milnor’s proof [4] uses a clever argument that relies on less sophisticated techniques. Ioppolo’s article aims to recount Milnor’s clever arguments.
						
						

							The article “Differential Inequalities” by Wilson Ong considers Gronwall’s inequality, a result originally published in [3]. Gronwall’s inequality is usually stated as follows [2]: Let ,  and  be real continuous functions defined in , with  for . Suppose that on  we have the inequality
						
						

							 
						
						

							Then
						
						

							 
						
						

							This inequality is clearly an inequality involving integrals. Ong’s exposition covers a variant of this, where the Gronwall type inequalities involve first and second order differential equations.
						
						

							In  of AfterMath, there’s an article on using the computer algebra system (CAS)  to study basic linear algebra. In this issue (issue 5), Alasdair McAndrew’s article “A Morsel of Maxima” presents an overview of this venerable CAS. McAndrew first provides a sketch of Maxima’s history, followed by an outline of Maxima’s state as of March 2007. The article then considers various interfaces to Maxima, including command line based and graphical user interfaces. The remainder of the article surveys features of Maxima that are useful in studying undergraduate mathematics: arithmetic, polynomial and rational functions, calculus, linear algebra, basic Maxima programming, solving equations, differential equations, plotting, and help and documentation.
						
						

							And there’s the article “Number Theory and Cryptography using PARI/GP” written by yours truly. It introduces various  commands that are useful for studying elementary number theory and undergraduate cryptography, in particular the RSA public key cryptosystem. But be warned that the exposition on RSA and cryptography is for educational purposes only, and readers who require more detailed discussions should consult specialized texts on the subject.
						
						

							Issue 5 also contains a section on mathematics jokes, as well as a section containing some basic mathematical problems.
						
						

							References 
						
						

							

								L. E. J. Brouwer. Uber Abbildung von Mannigfaltigkeiten. Math. Ann., pages 97-115, 1912.
							
							

								S. S. Dragomir. Some Gronwall Type Inequalities and Applications. Nova Science Publishers, Hauppauge, New York, USA, 2003.
							
							

								T. H. Gronwall. Note on the derivatives with respect to a parameter of the solutions of a system of differential equations. Ann. Math., 20(2):293-296, 1919.
							
							

								J. Milnor. Analytic proofs of the “Hairy Ball Theorem” and the Brouwer Fixed Point Theorem. The American Mathematical Monthly, 85(7):521-524, 1978.
							
							

								H. Poincare. Sur les courbes definies par les equations differentielles. J. Math. Pures. Appl., 4(1):167-244, 1885.