Archive

Archive for the ‘Maxima’ Category

Elementary number theory using Maxima

Prime numbers

You might remember that for any integer $n$ greater than 1, $n$ is a prime number if its factors are 1 and itself. The integers 2, 3, 5, and 7 are primes, but 9 is not prime because $9 = 3^2 = 3 \times 3$. 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 $n$:

(%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 Maxima to return a list of the first $n$ primes, where $n$ is a positive integer. Programming in the Maxima language is different from programming in other languages like C, C++, and Java. 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, whereas last() 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 &lt; 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 &lt; 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 factorization is about breaking up an integer into smaller components. In number theory, 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 $2^1 = 2$ and $5^1 = 5$. When you multiply these two prime factors together, you end up with $10 = 2^1 \times 5^1$. The expression $2^1 \times 5^1$ is called the prime factorization of 10. Similarly, the expression $5^2$ is the prime factorization of 25, and $2^3 \times 3^2$ is the prime factorization of 72. Greatest common divisors Closely related to integer factorization is the concept of greatest common divisor (GCD). The Maxima command gcd() is able to compute the GCD of two expressions $e_1$ and $e_2$ 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 $a$ and $b$ can be recursively defined as follows: $\gcd(a,b) = \begin{cases} a & \text{if } b = 0 \\ \gcd(b, a \bmod b) & \text{otherwise} \end{cases}$ where $a \bmod b$ is the remainder when $a$ is divided by $b$. The above recursive definition can be easily translated to a Maxima function for integer GCD as follows (credit goes to amca 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 extended Euclidean algorithm provides an interesting relationship between $\gcd(a,b)$, and the pair $a$ and $b$. 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  String processing with Maxima 19 August 2009 Leave a comment The Maxima module stringproc.lisp provides many functions for working with and processing strings. You load the string processing module using the command load(stringproc). A string can be thought of as a sequence of characters delimited by double quotation marks. Something like “abc” is a valid string, but ‘abc’ is not a valid string because the sequence of characters abc is delimited by single quotation marks: (%i1) "abc"; (%o1) abc (%i2) 'abc'; Incorrect syntax: ' is not an infix operator 'abc'; ^  When studying cryptography, it is useful to know about the ASCII integer corresponding to a character. The ASCII integer corresponding to the capital letter “A” is 65, that of its lower case equivalent is 97. Here is a list of the letters of the English alphabet and their corresponding ASCII integers: A 65 a 97 B 66 b 98 C 67 c 99 D 68 d 100 E 69 e 101 F 70 f 102 G 71 g 103 H 72 h 104 I 73 i 105 J 74 j 106 K 75 k 107 L 76 l 108 M 77 m 109 N 78 n 110 O 79 o 111 P 80 p 112 Q 81 q 113 R 82 r 114 S 83 s 115 T 84 t 116 U 85 u 117 V 86 v 118 W 87 w 119 X 88 x 120 Y 89 y 121 Z 90 z 122  You can work with these integers, performing various arithmetic operations on them. However, it is easier to work with the upper case letters and reduce their corresponding ASCII integers to values between 0 and 25, inclusive. Let’s load the module stringproc.lisp and use the function charat() to get characters of a string at specific indices: (%i12) load(stringproc)$
(%i13) strUpper : "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
(%o13)                    ABCDEFGHIJKLMNOPQRSTUVWXYZ
(%i14) charat(strUpper, 1);
(%o14)                                 A
(%i15) charat(strUpper, 2);
(%o15)                                 B
(%i16) charat(strUpper, 3);
(%o16)                                 C


You should take note that Maxima starts counting from 1, not from 0.

The function cint() returns the ASCII integer corresponding to a character. To obtain the ASCII integers corresponding to the upper case English letters, you can use the function charlist() to convert a string to a list of characters, then apply cint() to each of those characters. Again, the function maplist() is very handy in this case. Recall that maplist() takes two arguments: the first argument is a function, in this case it’s the function cint(); and the second argument is a list of values that you want to apply that function to. Here is how you can obtain the ASCII integers corresponding to all the capital letters of the English alphabet:

(%i23) cint("A");
(%o23)                                65
(%i24) cint("B");
(%o24)                                66
(%i25) cint("C");
(%o25)                                67
(%i26) L : charlist(strUpper);
(%o26) [A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, T, U, V, W,
X, Y, Z]
(%i27) maplist(cint, L);
(%o27) [65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81,
82, 83, 84, 85, 86, 87, 88, 89, 90]


This can be done in one line:

(%i28) maplist(cint, charlist(strUpper));
(%o28) [65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81,
82, 83, 84, 85, 86, 87, 88, 89, 90]


Now to convert each of the above ASCII integers to values between 0 and 25, inclusive, you need to subtract 65 from each of those values. One way to do so is to define a function, call it to26(), that converts a character to its ASCII integer and subtract 65 from that integer. The result from to26() is an integer between 0 and 25, inclusive. Then you use maplist() to apply that function to a list of characters:

(%i32) to26(x) := cint(x) - 65$(%i33) to26("A"); (%o33) 0 (%i34) to26("B"); (%o34) 1 (%i35) to26("C"); (%o35) 2 (%i36) maplist(to26, charlist(strUpper)); (%o36) [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]  Another way to do this is to use an un-named function, which can be achieved using the command lambda(): (%i37) maplist(lambda([x], cint(x) - 65), charlist(strUpper)); (%o37) [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]  What about going the other way round? Say you have an integer between 0 and 25, inclusive. You want to convert that integer to its corresponding upper case letter. To do so, you reverse the above process as follows. If $n$ is an integer between 0 and 25, add 65 to $n$ and use the function ascii() to convert the result to an ASCII character: (%i38) n : 0; (%o38) 0 (%i39) ascii(n + 65); (%o39) A (%i40) n : 1; (%o40) 1 (%i41) ascii(n + 65); (%o41) B (%i42) n : 2; (%o42) 2 (%i43) ascii(n + 65); (%o43) C (%i44) L : maplist(lambda([x], cint(x) - 65), charlist(strUpper)); (%o44) [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] (%i45) maplist(lambda([x], ascii(x + 65)), L); (%o45) [A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, T, U, V, W, X, Y, Z]  Notice that the result above is a list of characters. To convert that into a string, use the function simplode(): (%i46) simplode( maplist(lambda([x], ascii(x + 65)), L) ); (%o46) ABCDEFGHIJKLMNOPQRSTUVWXYZ  These are just some of the commands in the string processing module stringproc.lisp of Maxima. For more string processing functions, have a look at chapter 76 of the Maxima manual. Introducing Maxima 17 August 2009 3 comments Maxima is a computer algebra system. In other words, it is a mathematics software that is capable of symbolic manipulation. It is available for Windows as a command line interface or as a graphical interface via wxMaxima: You can also use Maxima under Linux: For a short introduction to Maxima, see this 10-minute tutorial. The Maxima website contains many types of documentation for learning Maxima, ranging from short tutorials to the Maxima manual. Numbers and Arithmetics Think of Maxima as a powerful calculator. To add two numbers, you use the addition operator “+”. You subtract two numbers using the subtraction operator “-”, multiply two numbers using “*”, and divide two numbers using “/”. Here are some examples of these in action: (%i1) 3 + 5; (%o1) 8 (%i2) 4 * 6; (%o2) 24 (%i3) 9 - 7; (%o3) 2 (%i4) 12 / 6; (%o4) 2  For now, don’t worry too much about things like “(%i1)” and “(%o1)”. But you should know that “(%i1)” means input number 1, while “(%o1)” means output number 1. Exponentiation can be done using “^”. The order in which exponentiation is performed is important. You can control the order of operation using parentheses: (%i8) 2^3; (%o8) 8 (%i9) (2^3)^2; (%o9) 64 (%i10) 2^(3^2); (%o10) 512 (%i11) 2^3^2; (%o11) 512  Use the exclamation mark “!” to do factorials: (%i12) 2!; (%o12) 2 (%i13) 3!; (%o13) 6 (%i14) 4!; (%o14) 24  Lists An easy way to create lists is via the opening square bracket “[" and the closing square bracket "]“. Here is a list of numbers: (%i16) L : [2, 4, 6, 8]; (%o16) [2, 4, 6, 8]  The operator “:” is used for assigning something to a variable. In the above, we assign the list [1, 2, 3, 4] to the variable L. A list can have a mixture of numbers and strings: (%i17) L2 : [1, 2, 3, "a", "b", "c"]; (%o17) [1, 2, 3, a, b, c]  Unlike some programming languages (e.g. C, C++, Java, Python), lists defined in Maxima have their indices starting from 1. For the list L above, the first element 2 has index 1, the second element 4 has index 2, the third element 6 has index 3, and so on. The elements of a list are accessed using square brackets: (%i22) L[1]; (%o22) 2 (%i23) L[2]; (%o23) 4 (%i24) L[3]; (%o24) 6 (%i25) L[4]; (%o25) 8  Use the command length() to find out the number of elements in a list. This command can also be used to loop through all elements of a list: (%i29) length(L); (%o29) 4 (%i30) for i:1 thru length(L) do print(L[i])$
2
4
6
8


Notice that the dollar sign “$” was used to suppress any output other than those of the command print(). Sometimes the dollar sign is used to suppress output because you might not be interested in viewing the output or the output is too long. Try these: (%i31) 3^(4^5); (%o31) 37339184874102004353295975418486658822540977678373400775063693172207904\ 061726525122999368893880397722046876506543147515810872705459216085858135133698\ 280918731419174859426258093880701995195640428557181804104668128879740292551766\ 801234061729839657473161915238672304623512593489605859058828465479354050593620\ 237654780744273058214452705898875625145281779341335214192074462302751872918543\ 286237573706398548531947641692626381997288700690701389925652429719852769874927\ 4196276811060702333710356481 (%i32) 3^(4^5)$


Lists can also be created using the command makelist(). A general syntax of this command is makelist(expr, i, i_0, i_n) where expr is an expression involving the variable i, with i starting from i_0 and ends at i_n. For example, here is how you can create a list of integers from 0 to 10:

(%i34) makelist(i, i, 0, 10);
(%o34)                [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]


In this case our expression is just i, with the variable i going from 0 to 10. Here is a list of positive integer squares less than 30:

(%i36) makelist(i^2, i, 1, 5);
(%o36)                         [1, 4, 9, 16, 25]


and a list of reciprocals:

(%i37) makelist(1 / (i+2), i, 1, 10);
1  1  1  1  1  1  1  1   1   1
(%o37)                 [-, -, -, -, -, -, -, --, --, --]
3  4  5  6  7  8  9  10  11  12


In the last example, the expression is 1 / (i+2), the variable is i, which goes from 1 to 10.

The command sum() can add up all numbers in a range of values. You can sum all numbers within a range as follows:

(%i41) sum(i, i, 0, 10);
(%o41)                                55
(%i42) sum(i^3, i, 1, 5);
(%o42)                                225
(%i43) sum(1 / (i+2), i, 1, 10);
44441
(%o43)                               -----
27720


The command sum(i, i, 0, 10) implements the expression

$\displaystyle{\sum_{i=0}^{10} i}$

The command sum(i^3, i, 1, 5) implements the expression

$\displaystyle{\sum_{i=1}^{5} i^3}$

And the command sum(1 / (i+2), i, 1, 10) implements the expression

$\displaystyle{\sum_{i=1}^{10} \frac{1}{i+2}}$

The command sum() can also work with lists. The last three summations can also be done as follows:

(%i44) L : makelist(i, i, 0, 10);
(%o44)                [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
(%i45) L2 : makelist(i^3, i, 1, 5);
(%o45)                        [1, 8, 27, 64, 125]
(%i46) L3 : makelist(1 / (i+2), i, 1, 10);
1  1  1  1  1  1  1  1   1   1
(%o46)                 [-, -, -, -, -, -, -, --, --, --]
3  4  5  6  7  8  9  10  11  12
(%i47) sum(L[i], i, 1, length(L));
(%o47)                                55
(%i48) sum(L2[i], i, 1, length(L2));
(%o48)                                225
(%i49) sum(L3[i], i, 1, length(L3));
44441
(%o49)                               -----
27720


The command prod() can be used to multiply a range of numbers together. For instance, the expression

$\displaystyle{\prod_{i=1}^{10} i}$

can be computed in Maxima as

(%i58) prod(i, i, 1, 10);
(%o58)                              3628800


The expression

$\displaystyle{\prod_{i=1}^{5} i^3}$

can be evaluated as

(%i59) prod(i^3, i, 1, 5);
(%o59)                              1728000


And the expression

$\displaystyle{\prod_{i=1}^{10} \frac{1}{i+2}}$

is computed as

(%i60) prod(1 / (i+2), i, 1, 10);
1
(%o60)                             ---------
239500800


All of these three expressions can also be computed as follows:

(%i61) L : makelist(i, i, 1, 10);
(%o61)                  [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
(%i62) L2 : makelist(i^3, i, 1, 5);
(%o62)                        [1, 8, 27, 64, 125]
(%i63) L3 : makelist(1 / (i+2), i, 1, 10);
1  1  1  1  1  1  1  1   1   1
(%o63)                 [-, -, -, -, -, -, -, --, --, --]
3  4  5  6  7  8  9  10  11  12
(%i64) prod(L[i], i, 1, length(L));
(%o64)                              3628800
(%i65) prod(L2[i], i, 1, length(L2));
(%o65)                              1728000
(%i66) prod(L3[i], i, 1, length(L3));
1
(%o66)                             ---------
239500800


Functions and Maps

Maxima uses the “colon-equals” operator “:=” in order to create functions. Here is a function that returns 1 if the input is a prime number, and returns 0 otherwise. We then evaluate the function for some integer values.

(%i67) func(x) := if primep(x) then 1 else 0$(%i70) func(1); (%o70) 0 (%i71) func(2); (%o71) 1 (%i72) func(3); (%o72) 1 (%i73) func(4); (%o73) 0 (%i74) func(5); (%o74) 1 (%i75) func(6); (%o75) 0 (%i76) func(7); (%o76) 1 (%i77) func(8); (%o77) 0 (%i78) func(9); (%o78) 0  You can achieve the same effect by using the command maplist(): (%i80) maplist(func, makelist(i, i, 1, 9)); (%o80) [0, 1, 1, 0, 1, 0, 1, 0, 0]  As you can see, maplist() applies a function to every element of a list and then returns the results as another list. But you don’t have to create the function func() first. You can also use the command lambda() to achieve the same effect: (%i81) maplist(lambda([x], if primep(x) then 1 else 0), makelist(i, i, 1, 9)); (%o81) [0, 1, 1, 0, 1, 0, 1, 0, 0]  As a final example, suppose you have a list of numbers from 1 to 100. If any number in the list is prime, you change that number to 1. If the number is not prime, you change it to 0. At the end of this computation, you would end up with a list of the same length as the original list, but the elements of the new list are 1′s and 0′s. To find the number of primes between 1 and 100, you sum the new list. Here is how to do so: (%i85) L : makelist(i, i, 1, 100); (%o85) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100] (%i86) for i:1 thru length(L) do ( if primep(L[i]) then L[i] : 1 else L[i] : 0 )$
(%i87) L;
(%o87) [0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1,
0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0,
0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0]
(%i88) sum(L[i], i, 1, length(L));
(%o88)                                25


This means that there are 25 prime numbers between 1 and 100.