:tocdepth: 2

.. _frobenius_number:

The Frobenius number – Starting to Program 
============================================================================

.. contents:: :local:

.. highlight:: python





Exploration Question: 
----------------------

::

   Suppose you have coins of value a cents and b cents.  What amounts of money can you get by combining these?


   Q: Can you think of how to formulate our question mathematically?

   A: Given positive integers a and b, for which positive integers n does there exist non-negative integers x and  y such that n=ax+by?

We will solve this problem by first trying some easy cases such as
:math:`a=2`, :math:`b=3`, or :math:`a=3`, :math:`b=4`, or :math:`a=3`,
:math:`b=5` to begin. We will explore the first two cases :math:`a=2`,
:math:`b=3`, and :math:`a=3, b=4`, together so everyone get an idea of
what we are doing. We will then group ourselves in small groups to
explore more cases.

As we shall discover in our exploration, when there is such an
:math:`n`, there is also an :math:`N` such that for all :math:`n`
greater than or equal to :math:`N` will also work, but :math:`N-1` does
not. We will call :math:`N` the **conductor** of the set
:math:`\{a, b\}` and :math:`N-1` as the **Frobenius number**. (The
problem is known as the Frobenius problem or the coin problem as well.)

The **conductor** of two positive integers :math:`a, b` is the smallest
integer :math:`N` such that all :math:`n` greater than or equal to
:math:`N` can be written as :math:`n = ax + by`.

The **Frobenius** number of two positive integers :math:`a, b` is the
largest integer :math:`C` such that there do not exist positive integers
:math:`x` and :math:`y` with :math:`C = ax + by`.

Your goal is to experiment with this. Experiment! For a given pair of
positive integers a and b can you think of some questions we can ask?
For example we can ask things like the following:

-  Is there a conductor?
-  What is the conductor?
-  How many numbers n can be written ax+by?

   ::

      Can you think of more questions to ask? We will group ourselves into groups and try to explore our questions on simple cases as we did before. 

Now, soon we’ll be using the computer to help explore more complex
examples like :math:`a=19, b=17`. However, experimental mathematics
needs only curiosity and persistence. And you don’t need a computer to
do that.

Finally we will go to the computer laboratory and introduce *SageMath* -
you can use the Sage notebook below and upload it to try it yourself!

You should do this for tomorrow :

1. Is there a conductor for the sets :math:`\{2, 5\}`, :math:`\{3, 5\}`,
   and :math:`\{3, 6\}` ? If so, what is it? Explore!
2. Is there a Frobenius number for the sets :math:`\{2, 5\}`,
   :math:`\{3, 5\}`, and :math:`\{3, 6\}`? If so, what is it? Explore!
3. Can you list all the set of numbers that can not be written as
   :math:`2 x + 5 y`, :math:`3 x + 5 y` and :math:`3 x + 6 y`? When is
   this set finite?
4. Write down a step-by-step explanation of the conductor of
   :math:`\{3, 4\}`. You should both compute it and explain how you know
   that all N greater than the conductor really can be written as
   :math:`3x+4y` for some :math:`x,y`.
5. Try leaning how to load a sage jupyter notebook and how to start your
   own sage notebook.

We’ve already been trying out some experimental mathematics; now it’s
time to see how one might start doing this with the computer.  Again,
**please feel free to try things**.  Don’t just try to follow the
lecture.

 By the way, the acronym “SAGE” stands for “Software for Algebra and
Geometry Experimentation”.

We saw some proof earlier.  Let’s start our exploration of programming
by recalling a theorem some of you may not have known needs to be a
theorem.

--------------

**Theorem (Division Algorithm)**:

.. math:: \text{If }a\text{ and }b\text{ are integers, }b>0\text{, then we can always write }a=qb+r\text{ with }0\leq r<b\text{ and }q\text{ an integer.}

--------------

In words, if you divide an integer :math:`a` by a positive integer
:math:`b`, you will always get an integer remainder :math:`r` that is
nonnegative, but less than :math:`b`.  Further, the remainder
is \ *unique*.

For example, if :math:`a=13,b=3` then

.. math:: 13=4\cdot 3+1\text{, so }q=4\text{ and }r=1.

We can prove this using the definition of divisibility and the
“Well-Ordering Principle”, but I do not assume you are too familiar with
induction and this principle.

By the way, notice that this relates directly to something we proved
earlier; in this event, :math:`\gcd(a,b)=\gcd(b,r)`.  This leads
directly to the Euclidean algorithm for computing the GCD. We can check
this works using Sage.

The code ``a % b`` will return the remainder upon division of :math:`a`
by :math:`b`. In other words, the result is the unique integer :math:`r`
such that

1. $ 0 :raw-latex:`\leq `r < b$, and
2. :math:`a = bq + r` for some integer :math:`q` (the quotient), as
   guaranteed by the Division Algorithm. Then :math:`(a − r)/b` will
   equal :math:`q`. For example,

.. code:: ipython2

    r = 14 % 3
    r

.. code:: ipython2

    q = (14 - r ) /3
    q

It is also possible to get both the quotient and remainder at the same
time with the ``.quo_rem()`` method (quotient and remainder).

.. code:: ipython2

    a = 14
    b = 3
    a.quo_rem ( b )

A remainder of zero indicates divisibility. So ``(a % b) == 0`` will
return ``True`` if :math:`b` divides :math:`a`, and will otherwise
return False.

.. code:: ipython2

    (20 % 5) == 0

.. code:: ipython2

    (17 % 4) == 0

The ``.divides()`` method is another option:

.. code:: ipython2

    c = 5
    c.divides (20)

.. code:: ipython2

    d = 4
    d . divides (17)

The greatest common divisor of a and b is obtained with the command
``gcd(a, b)``, where in our first uses, a and b are integers.

.. code:: ipython2

    gcd (2776 , 2452)

We can use the gcd command to determine if a pair of integers are
relatively prime

.. code:: ipython2

    a = 31049
    b = 2105
    gcd (a , b ) == 1

.. code:: ipython2

    a = 3563
    b = 2947
    gcd (a , b ) == 1

The command ``xgcd(a,b)`` (“eXtended GCD”) returns a triple where the
first element is the greatest common divisor of a and b (as with the
gcd(a,b) command above), but the next two elements are values of
:math:`r` and :math:`s` such that :math:`ra + sb = gcd(a, b)`.

.. code:: ipython2

    xgcd (633 ,331)

Portions of the triple can be extracted using [ ] (“indexing”) to access
the entries of the triple, starting with the first as number :math:`0`.
For example, the following should always return the result True, even if
you change the values of :math:`a` and :math:`b`. Try changing the
values of :math:`a` and :math:`b` below, to see that the result is
always True.

.. code:: ipython2

    a = 633
    b = 331
    extended = xgcd (a , b )
    g = extended [0]
    r = extended [1]
    s = extended [2]
    g == r * a + s * b

The method ``.is_prime()`` will determine if an integer is prime or not.

.. code:: ipython2

    a = 117371
    a . is_prime()

.. code:: ipython2

    b = 14547073
    b . is_prime()

The command ``random_prime(a, proof=True)`` will generate a random prime
number between :math:`2` and :math:`a`. Experiment by executing the
following two compute cells several times. (Replacing ``proof=True`` by
``proof=False`` will speed up the search, but there will be a very,
very, very small probability the result will not be prime.)

.. code:: ipython2

    a = random_prime (10^21 , proof = True )
    a

.. code:: ipython2

    a.is_prime ()

The commands ``next_prime(a)`` and ``previous_prime(a)`` are other ways
to get a single prime number of a desired size. Try it Yourself.

In addition to checking if integers are prime or not, or generating
prime numbers, Sage can also decompose any integer into its prime
factors, as described by the Fundamental Theorem of Arithmetic

.. code:: ipython2

    a = 2600
    a.factor()

So $2600 = 2^3 × 5^2 × 13 $ and this is the unique way to write 2600 as
a product of prime numbers (other than rearranging the order of the
primes themselves in the product).

Lists And Loops 
----------------

Now let’s start talking the basics of how one might go about exploring
questions with *SageMath*.

Lists
~~~~~

The command ``prime_range(a, b)`` returns an ordered list of all the
primes from :math:`a` to :math:`b- 1`, inclusive. For example,

.. code:: ipython2

    prime_range (500 , 550)

While Sage will print a factorization nicely, it is carried internally
as a list of pairs of integers, with each pair being a base (a prime
number) and an exponent (a positive integer). Study the following
carefully, as it is another good exercise in working with Sage output in
the form of lists

.. code:: ipython2

    a = 2600
    a.factor ()

.. code:: ipython2

    a = 2600
    factored = a . factor ()
    first_term = factored [0]
    first_term

.. code:: ipython2

    second_term = factored [1]
    second_term

.. code:: ipython2

    third_term = factored [2]
    third_term

.. code:: ipython2

    first_prime = first_term [0]
    first_prime

.. code:: ipython2

    first_exponent = first_term [1]
    first_exponent

The next compute cell reveals the internal version of the factorization
by asking for the actual list. And we show how you could determine
exactly how many terms the factorization has by using the length
command, ``len()``.

.. code:: ipython2

    list(factored)

.. code:: ipython2

    len(factored)

Can you extract the next two primes, and their exponents, from a?

In *SageMath*, it turns out that making a list of lists is a good way to
organize our data, by using the table command.

.. code:: ipython2

    L = [[(3,4),6,3],[(4,5),12,6],[(3,6),'none',oo]]
    table(L,header_row=["a,b",'conductor','number not possible'])

.. code:: ipython2

    L = [["$(3,4)$","$6$","$3$"],["$(4,5)$","$12$","$6$"],["$(3,6)$",'none',"$\infty$"]]
    table(L,header_row=["$a,b$",'conductor','number not possible'])

This is a good use for lists - making things easy to communicate to
others!

Loops
~~~~~

A related concept you should have also seen in your previous Python
class is a *loop*. Loops make doing tedious things many times easier.
Here is an example getting the determinant of powers of a matrix. I
*could* do it “by hand”.

.. code:: ipython2

    A = matrix([[1,2],[3,4]])
    print det(A^0)
    print det(A^1)
    print det(A^2)
    print det(A^3)
    print det(A^4)

This is not terrible, but it’s not exactly nice either.  Instead, let’s
use a loop.

.. code:: ipython2

    for i in [0,1,2,3,4]:
        print det(A^i)

What did we do?  

For (each) :math:`i` in the list (ordered set) [0,1,2,3,4], (return)
:math:`A^i`. The square brackets created a list, and the powers of the
original matrix come in the same order as the list. Remember, the colon
in the first line and the indentation in the second line are
**extremely** important; they are the basic syntactical structure of
Python.

.. code:: ipython2

    for i in [0,1,2,3,4]
        A^i

.. code:: ipython2

    for i in [0,1,2,3,4]:
    A^i

Streamlining
~~~~~~~~~~~~

For the curious, it’s worth knowing there are quicker ways to make the
possible values for :math:`i` quicker to write.  Here are two possible
options.

.. code:: ipython2

    for i in [0..4]:
        print det(A^i)

.. code:: ipython2

    for i in range(5):
        print det(A^i)

Notice that the ``range(5)`` starts counting at zero and ends *before*
reaching five; this is standard behavior. You can get other behavior,
too.

.. code:: ipython2

    range(3, 23, 2); [3,5..21]

List Comprehensions
~~~~~~~~~~~~~~~~~~~

There is a very powerful way to create such lists in a way that looks
like a loop. This way also looks like mathematical notation. It is
called a *list comprehension*.  We start with a relatively easy example:

.. math:: \{n^2\mid n\in\ZZ, 3 \leq n \leq 12\}

 Who hasn’t written something like this at some point? Here’s how to
translate this into Sage:

.. code:: ipython2

    [ n^2 for n in [3..12] ]

The notation is easiest if you think of it mathematically; “The set of
:math:`n^2`, for (any) :math:`n` in the range between 3 and 12.”

This sort of construction works for lots of things of mathematical
interest, of course - such as our first example.

.. code:: ipython2

    [ det(A^i) for i in [0..4] ]

How could we use lists and loops to explore the conductor?

1. Can I list *all* numbers of the form :math:`3x+4y`?
2. How would I list some numbers of this form?  (Hint: fix one of the
   variables).
3. How would I list more numbers of this form? (Hint: make a loop
   inside of a loop).
4. Can you turn this into a list and not just a loop?

I know you have done dictionaries in your python course. Can you try to
use dictionaries for this? 

Defining Functions (Extending Sage)
-----------------------------------

Hopefully you will experience success with your exploration! Still, it
is kind of annoying that there isn’t just a simple command for it. In
fact, it is often the case that Sage can do something, but doesn’t have
a simple command for it. For instance, you might want to take a matrix
and output the square of that matrix minus the original matrix.

.. code:: ipython2

    A = matrix([[1,2],[3,4]])
    A^2 - A

How might one do this for other matrices?  Of course, you could just
always do :math:`A^2-A` again and again.  But this would be tedious and
hard to follow, as with so many things that motivate a little
programming.  Here is how Python and Sage solve this problem.

.. code:: ipython2

    def square_and_subtract(mymatrix):
        return mymatrix^2 - mymatrix

.. code:: ipython2

    square_and_subtract(A)

.. code:: ipython2

    square_and_subtract(matrix([[1.5,0],[0,2]]))

As a technical matter, this is only a \ *Python function*, not a Sage
symbolic function (callable expression).  That means it does not have
access to all the same things as we discussed earlier today.  That is
okay, as it is often better for a function not to live in the symbolic
algebra world.

What if we are worried about forgetting what this function does?  To be
fair, we chose a great name for it.  Still, just in case, we can provide
a documentation string, putting it in triple quotes ``"""``.

.. code:: ipython2

    def square_and_subtract(mymatrix):
        """
        Return `A^2-A`
        """
        return mymatrix^2-mymatrix

.. code:: ipython2

    square_and_subtract?

That’s pretty cool!  And potentially quite helpful, especially if the
function is complicated.  The :math:`A` typesets properly because we put
it in backticks ``A``. (For the \ *real* experts, one can use “raw
strings” to include backslashes (say, for LaTeX) in these documentation
strings, like r""“:raw-latex:`\frac{a}{b}`”"".  Look at the
documentation for Bessel functions for some great examples). A very
careful reader *may* have noticed that there is nothing that requires
the input ‘mymatrix’ to be a matrix.  Sage will just try to square
whatever you give it and subtract the original thing.

.. code:: ipython2

    square_and_subtract(sqrt(5))

Functions are very flexible in what input they can allow or require, as
well as what output they give.  Below are three examples that show some
of this flexibility.   There are many other good resources out there.

.. code:: ipython2

    def func1( y, z ):
        return y+z

.. code:: ipython2

    def func2( y, z=3):
        return y+z

.. code:: ipython2

    def func3( y ):
        return y, y+3

Could we use functions to explore our question? 1. Could you write a
function that lists all numbers of the form :math:`3x+4`, given some
range of :math:`x`? 2. What about of the form :math:`3x+8=3x+4\cdot 2`?
3. What about of the form :math:`3x+4y`?  Note that here you will have
TWO ranges to consider, that of :math:`x` and that of :math:`y`. 4. Are
your lists sorted or not?  Can you think of a way to fix that? 5. What
about for the form :math:`4x+5y`? 6. What about for :math:`ax+by`?

Exercises 
----------

These exercises are about investigating basic properties of the
integers, something we will frequently do when investigating groups.

1. Use the ``next_prime()`` command to construct two different 8-digit
   prime numbers and save them in variables named a and b.

2. Use the ``.is_prime()`` method to verify that your primes a and b are
   really prime.

3. Verify that 1 is the greatest common divisor of your two primes from
   the previous exercises.

4. Find two integers that make a “linear combination” of your two primes
   equal to 1. Include a verification of your result.

5. Determine a factorization into powers of primes for c = 4598037234

6. Write a computer program that will implement the Euclidean algorithm.
   The program should accept two positive integers a and b as input and
   should output gcd(a, b) as well as integers :math:`r` and :math:`s`
   such that

--------------

.. raw:: html

   <!--NAVIGATION-->

< `2. Review of Python
programming <ems_2020_day_02_short_review_on_python_programming.ipynb>`__\ \|\ `Table
of contents <ems_2020_table_of_contents.ipynb>`__ \| `4. Starting To
Program <ems_2020_day_04_sum_of_squares_and_plotting>`__ >

--------------