zoukankan      html  css  js  c++  java
  • 使用maxima解决初等数论中的问题:


                                                   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:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    (%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:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    (%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 nprimes, where n is a positive integer. Programming in the Maxima language is different from programming in other languages like CC++, 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.

    1

    2

    3

    4

    5

    6

    7

    8

    9

    (%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:

    1

    2

    3

    4

    5

    6

    7

    8

    (%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.

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    (%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:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    15

    16

    17

    18

    19

    20

    21

    22

    /* 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():

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    15

    16

    (%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:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    (%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:

    1

    2

    3

    4

    5

    6

    (%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):

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    15

    16

    17

    /* 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:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    (%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:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    15

    16

    17

    18

    19

    20

    21

    22

    23

    24

    /* 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:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    (%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.

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    (%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

     
  • 相关阅读:
    Core2.0知识整理
    bootbox.js官方文档
    MongoDB分片详解
    MongoDB高可用集群+MMS集群监控搭建
    使用 AcceptTcpClientAsync 进行 异步 操作
    Windows证书操作
    SQL Server查询所有的表名、字段名、注释
    【转载】ASP.NET Core Web 支付功能接入 微信-扫码支付篇
    【转载】ASP.NET Core Web 支付功能接入 支付宝-电脑网页支付篇
    C#网络编程系列文章
  • 原文地址:https://www.cnblogs.com/enjoy233/p/2998405.html
Copyright © 2011-2022 走看看