Tuesday, November 10, 2009

Mersenne Numbers - short sage example

I have used Sagemath and Pari/GP in the past and have found them both very capable Mathematics tools.

Pari/GP has been around a bit longer*, and has many recorded sample code snippets in The On-Line Encyclopaedia of Integer Sequences which AT&T kindly host.

(*Pari/GP has been GPL licensed for free use and modification for the last 15 years or more)

Sage is a bit different than a dedicated number theory program, and could perhaps be described in several ways:
  • A unifying front-end to many computerised algebra and mathematical programs
  • A means of accessing lots of computerised systems using a common language (python)
  • A front-end to Pari/GP (if you are doing Number Theory)
Anyway back to the maths...

See example output from below Sage script on sagemath cloud.

# Most Positive Integers of the form -1+2**n are Composite rather than Prime
#
# There is no significance to the range of numbers 7 through 37
# I selected that range so the primality test is_prime would not
# take too long.
#
# There are more than Mersenne Composites and Mersenne Primes generated by this
# loop. For example 8 gives -1+2**8 = 255 which is not a Mersenne Composite
# nor a Mersenne Prime.
#
# Expect many Falses and just a few Trues (Mersenne Primes) in the output.
#
for n in range(7,37):
mn=-1+2**n
print mn, is_prime(mn)
# True and False are intended to indicate if the number is Prime
# False indicates that the preceding number is a Composite

There is very little involved, and I could have written just the 3 lines and left it uncommented, however a few explanatory comments will hopefully illustrates things a bit better.

What inspired this example was some tabulation I did a while back; an extract is pasted below (click to enlarge):


If you have worked with Mersenne numbers before, then you will have spotted, I hope, that the example loop I used in sage is a bit too simple.

What makes a Mersenne number of any kind, is a prime power of two rather than just any old power of 2. So to be strictly useful, the sage example should only iterate through prime numbers for n, rather than all numbers in the range 7..37.

The great thing about an online version of sage is that anybody can register for an account and submit an example. Please feel free to alter the example I have given here to make it work better as just described, and try sage notebook yourself.
( Easy to take a copy of what I have done and paste it into your own notebook. )

Several universities already host sage online (convenient student access), although it is simple to install on any Ubuntu laptop if you want your own installation rather than relying on the Internet.

Downloads of Sage mathematics for other systems are available here.

Having read this article, and in particular my words on how to improve the example I gave, you might want to come up with your own enhancement to the opening line of Wikipedia's entry (2009 current):
In mathematics, a Mersenne number is a positive integer that is one less than a power of two

Finally a tip for any Pari/GP users who try Sage mathematics:
You may be using is_prime() rather than writing in pari isprime()
and here is quick extract to confirm what you are getting:
sage.rings.arith.is_prime(n, flag=0)
Returns True if prime, and False otherwise. The result is proven correct - this is NOT a pseudo-primality test!.

Friday, May 1, 2009

python and divisor summing

There are lots of areas of number theory concerned with exact divisors of integers - perfect numbers, amicable numbers, sociable numbers and much more.

A good place to begin understanding the behaviour of the divisor sum functions
[ divsum() in maxima, sigma() in pari/gp ] is to consider the integers 1 to 100.

A simple python program to show you what happens for
those integers is given below.

#!/usr/bin/python
## Aliquot sequence demonstrator
## This program was written as a very simple demonstrator of what
## happens for numbers in the range 1 through 100 as, for each,
## you sum the 'proper' divisors and repeat.
## In order to make viewing on a terminal output practical the
## routines work in batches of 10 (or 20) at a time as
## required.
## A crude approach to knowing when to stop printing relies in
## part on a list of known terminators.
## A more general program would instead maybe look at how many
## repetitions of the same value had occurred instead. Try that
## once you have decided how many repetitions you think are
## appropriate.
## Further reading on aliquot sequences might be had by searching:
##  The Lehmer Five are the first unterminated aliquot sequences
##
## (c) 2009 Gary Wright
##
##This program is free software: you can redistribute it and/or modify
##it under the terms of the GNU General Public License as published by
##the Free Software Foundation, either version 3 of the License, or
##(at your option) any later version.
##
##This program is distributed in the hope that it will be useful,
##but WITHOUT ANY WARRANTY; without even the implied warranty of
##MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##GNU General Public License for more details.
##
##You should have received a copy of the GNU General Public License
##along with this program.  If not, see http://www.gnu.org/licenses/
##
## http://www.wrightsolutions.co.uk/guide/code
## V1.0 20090501GW
##
#import divsum_out as ds
#ds.divsum_out20(range(20,40))# Readable when do divsum_out20(range(20)) ...
#...but not much higher
#ds.divsum_out10(0);divsum_out10(10)# More reabable and presentation suitable

known_terminators = ( 0, 6, 28 )     # tuple as immutable data

def divsum_out(n1):
print "%s%3i%s" % ('#',n1,'#')   # For first entries/line use
                                 # bang rather than pipe (|)
                                 # so easy to distinguish the first entries
                                 # which are just the numbers themselves.
#printf( ,'|',n1,'|')
while(n1 not in known_terminators):
    n1 = divsum(n1)
    print "%s%3i%s" % ('|',n1,'|')
return

def divsum(n):
# Sum all those values in the range 1..n-1 that are aliquot divisors of n
return sum(i for i in range(1, n-1) if n % i == 0)

def divsum_out20(nlist):
""" Run through divisor sums for 20 integers given as an input list
"""
delim = '#'
still_running_counter = 20
while(still_running_counter > 0):
    line = ''
    for i in nlist:
        line += (len(line) == 0)*delim
        line += str(i).rjust(3)+delim
    print line
    still_running_counter = sum(i for i in nlist if i not in known_terminators)
    delim = '|'
    nlist = map(divsum, nlist)
return

def divsum_out10(sval=0):
""" Run through divisor sums for 10 integers beginning with the supplied integer sval
sval is an abbreviation for Start value and defaults to 0.
You might try 10 or 20.
"""
delim = '#'
nlist = range(sval,sval+10)
still_running_counter = 10
while(still_running_counter > 0):
    line = ''
    for i in nlist:
        line += (len(line) == 0)*delim
        line += str(i).rjust(3)+delim
    print line
    still_running_counter = sum(i for i in nlist if i not in known_terminators)
    delim = '|'
    nlist = map(divsum, nlist)
return


I ran this as follows:
python -i divsum_out.py
>>> divsum_out10(0)

If you are new to python then you should refrain from typing the '>>>' which I included just to indicate (as python does) that you are then within an interactive python shell.

If you prefer to do your maths in maple then you might websearch for aliquotparts and aliquotsum functions.

Note: The web url and mail address given in the code are yet to be operational and will be usable when the server setup is complete in late summer 2009. I include them for future reference only.

Wednesday, April 29, 2009

maxima calg and eigenvalues

Early matrix learning was supported by some laptop work involving the open source computer algebra system (calg) maxima.

I find the system to be quite capable and for more advanced use it has some great plotting facilities also.

Determinants and Eigenvalues are useful things to know about any matrix and the maxima help is shown below followed by some simple examples I tried.

? eigenvalues;

-- Function: eigenvalues ()
-- Function: eivals ()
Returns a list of two lists containing the eigenvalues of the
matrix . The first sublist of the return value is the list of
eigenvalues of the matrix, and the second sublist is the list of
the multiplicities of the eigenvalues in the corresponding order.

`eivals' is a synonym for `eigenvalues'.

`eigenvalues' calls the function `solve' to find the roots of the
characteristic polynomial of the matrix. Sometimes `solve' may
not be able to find the roots of the polynomial; in that case some
other functions in this package (except `innerproduct',
`unitvector', `columnvector' and `gramschmidt') will not work.

In some cases the eigenvalues found by `solve' may be complicated
expressions. (This may happen when `solve' returns a
not-so-obviously real expression for an eigenvalue which is known
to be real.) It may be possible to simplify the eigenvalues using
some other functions.

The package `eigen.mac' is loaded automatically when `eigenvalues'
or `eigenvectors' is referenced. If `eigen.mac' is not already
loaded, `load ("eigen")' loads it. After loading, all functions
and variables in the package are available.

Here is the output for determinant and eigenvalues of some simple 2x2 matrices which
are deliberately simple in order to get familiar with the eigenvalue
output form described in the help above.

batching 0607MatDeterminantZeroAndEigenvalues-Simple2x2triangular.mac
(%i2) m0 : matrix([0, 0], [0, 0])
[ 0 0 ]
(%o2) [ ]
[ 0 0 ]
(%i3) determinant(m0)
(%o3) 0
(%i4) eigenvalues(m0)
(%o4) [[0], [2]]
(%i5) m1 : matrix([1, 0], [0, 0])
[ 1 0 ]
(%o5) [ ]
[ 0 0 ]
(%i6) determinant(m1)
(%o6) 0
(%i7) eigenvalues(m1)
(%o7) [[0, 1], [1, 1]]
(%i8) m2 : matrix([0, 1], [0, 0])
[ 0 1 ]
(%o8) [ ]
[ 0 0 ]
(%i9) determinant(m2)
(%o9) 0
(%i10) eigenvalues(m2)
(%o10) [[0], [2]]
(%i11) m3 : matrix([0, 0], [1, 0])
[ 0 0 ]
(%o11) [ ]
[ 1 0 ]
(%i12) determinant(m3)
(%o12) 0
(%i13) eigenvalues(m3)
(%o13) [[0], [2]]
(%i14) m4 : matrix([0, 0], [0, 1])
[ 0 0 ]
(%o14) [ ]
[ 0 1 ]
(%i15) determinant(m4)
(%o15) 0
(%i16) eigenvalues(m4)
(%o16) [[0, 1], [1, 1]]

As described in the help the 2nd of the two lists which eigenvalues() returns give you the multiplicities of those eigenvalues.

All of these matrices I have used are Triangular and so the determinant is the product of the entries on the main diagonal.
For all the examples I have used expect determinant zero.