Wednesday, August 17, 2011

gmp mpz gotchas - part 1

The GNU mp (gmp) library is a fast and stable library for working with large integers.

As with all software, when you begin using it, you find out things that you should and shouldn't do.


NULL arguments to functions - check for each function:

mpz_get_str() is a function which seems to handle a NULL first argument just fine.

Example function call:

outM = mpz_get_str(NULL,10,M);

(outM above is a char pointer and gets the decimal representation of mpz M)

In the documentation, the function signature for mpz_get_str() is given as:

char * mpz_get_str (char *str, int base, mpz_t op )

Second example - be careful what you assume!

if (mpz_mod_ui(NULL,C,6) == 0) { /* do something */ }

The above throws a runtime error - not surprising when you consult the documentation for the function signature.

unsigned long int mpz_mod_ui (mpz_t r, mpz_t n, unsigned long int d ) 

mpz_mod_ui() is quite clear that the first argument should be an mpz type, so don't disappoint it :)

When working with gmp, you are effectively working in two type spaces - C and native gmp types. Functions such as mpz_mod_ui() are really flexible, in that they give you the mod() result in two different type spaces.

Libraries are not mind readers, and cannot know in which type space your next line of code will be working. By giving you the result as both an unsigned long int and the mpz_t labeled here as r, the function is doing it's damnedest to help you out!

The downside is that you may well see little point in retaining the result in the other type space.

My runtime error was caused by me attempting to throw away the mpz_t result as it felt redundant.

Here is my new preferred way - use an mpz_t you set up particularly for transient results.

if (mpz_mod_ui(MOD_NULL,C,6) == 0) {
 /* flag as six multiple */
}

Just create your own mpz_t MOD_NULL and use it where you are not interested in keeping the result. You will need to run mpz_init(MOD_NULL) before executing the if statement illustrated.



Transient variables - bad coding practices to avoid:

In the past, I have been tempted to reuse spare variables - example:
    if (mpz_mod_ui(BSQR,C,6) == 0) {

It is bad coding practice to use a seemingly proper variable, as a holder for transient data / as a dump.

Using dump variables generally is not without issue, however naming the thing correctly is at least a step in the right direction :)


The important thing to remember is that whatever you call your transient variable, that it is not used to hold values that really are used in your logic.

Avoid the temptation to reuse a label wherever possible - it rarely leads to good things.


NULL as the first argument to mpz_get_str - probably going wrong here:

Earlier I used mpz_get_str() as an example of a function that worked with NULL as a first argument.

outBSQRM4AC = mpz_get_str(NULL,10,BSQRM4AC);

Knowing the function signature, and looking at the above line of code, you should be thinking char pointer for outBSQRM4AC.

Assigning the result of mpz_get_str() using an equals symbol is okay, and you are asking the gmp library to do your allocation for you.

If however, you known the maximum size of the MPZ / character array that you require, then you may prefer to predefine your character array, and tell mpz_get_str() to use it. In that case the the first argument of mpz_get_str() will not be NULL.

If you have predefined your character array (including size, then you might instead be using a call of this form:

mpz_get_str(outBSQRM4AC,10,BSQRM4AC);


Notes and Further Reading:

The publisher SAMS have a book "Teach yourself C in 24 hours"

In the Chapter "Hour 17 - Allocating Memory" there is a discussion about malloc() and free() and how to use them in conjunction with character pointers.







A more indepth volume is titled "Teach yourself C in 21 days", but it really depends on how far you want to take your C skills, as to which is the better option.













Neither of these books are particularly 'new', however that is one of the plus points of C, it has a stable api and a body of literature that never goes out of date.

Neither of these books will teach you C++, however I have not and probably never will have, any interest in learning C++

If you do, then the 'in 24 hours', probably does not come near to describing, the time that you will have to devote to becoming skilled in C++.

Better C++ literature is available, if you decide ANSI C is not where you wish to spend your time.

Saturday, August 13, 2011

Cython - notes and gotchas - part 1

Experimenting with Cython, as bringing together many variants of the same .c file, has me yearning to find a better Python / C balance.

Working through tutorials and experiments, I will note anything tricky, or any easy to fall into traps.


TypeError: 'NoneType' object is not callable:

Because Python is flexible, and does not insist on you declaring variables ahead of time, you may be tempted to do the same in Cython.


Better to declare the variable 'primes' explicitly, or to avoid name clashes.


So in the second example the result of the function call ends up in 'primes_list_strings' and we work on it from there.

To clarify the original error: Calling the variable that holds the result 'primes', and making a function call to primes() in the same line, might lead to confusion.


Notes and Further Reading:

If you want completion features when working with Python, then one solution is to invoke iPython and use that.

There are several full IDEs for Python including Idle and Eric, but if you just need a little completion (and syntax highlighting), then iPython is okay also.


Screenshot of Eric from screenshots.debian.net:


Sunday, June 26, 2011

how many bits - are you 32 bit or 64 bit

My laptop usually has a 32 bit GNU / Linux install and a 64 bit install also.

It pays to remember which you are logged into when using openpfgw

-su: ./pfgw64: cannot execute binary file

The answer as to why you might get that response is simple, however it is easy to mistake the symptom.

Running off and websearching 'cannot execute binary file', was the wrong thing to do, and assuming that 'su' was somehow broken was a red herring.

Logging into a console as root and running ./pfgw64 gives a more meaningful:

-bash: ./pfgw64: cannot execute binary file

...which should be enough to tell you that the problem is not 'su', which in my case is just sitting in front of bash and relaying the message.

I hinted at the real cause earlier. Running ./pfgw32 will see a better result, because you are logged into a 32 bit system!

It is great that openpfgw provide optimised binaries for 32 bit and 64 bit, and it is up to you, as the user, to select the correct one!

Remember that openpfgw will also needs some libgmp libraries in order to run, otherwise expect a complaint about a missing .so file.

Wednesday, April 20, 2011

making primes - fun with niece and nephew

You cannot make a prime by multiplying other numbers together

Whichever naive definition* you choose to give, some examples always help.

Once the students have an idea of what is and is not a prime, a fun exercise can be to pick a significant year (Man on the moon), and work some addition and subtraction.

Using just the individual digits of 1969, challenge the student to see, how many primes they can construct, using addition and subtraction.

This exercise is a reinforcing mechanism in that instead of focusing on multiplication, you are working a secondary construction method.

  • 2 = 9 - 6 - 1
  • 3 = 9 - 6
  • 5 = 6 - 1
  • 7 = 6 + 1
  • 11 = 9 + 9 - 1 - 6
  • 13 = 9 + 9 + 1 - 6
  • 17 = 9 + 9 - 1
  • 19 = 9 + 9 + 1

Doing this with a group of students, it is more likely that at least one will spot that a +1 in the construction, can be switched to -1 in the construction, and often obtain a new prime.

Obviously this is not a general rule, however it does introduce the student to the fact that sometimes when you have a prime, another can be found, just two away.

Now use another significant year, say 1989 (Fall of the Berlin Wall), and repeat the exercise.
( Not so easy with the 1989 example )

And use the year of birth for some of the students, see how many constructions can be made.

A useful follow up discussion can be to examine why 1969 is a better specimen than 1989, for this particular exercise.

How about 1979? Better? Worse?

*Note: I used the phrase naive definition at the beginning. With younger learners it is sometimes a useful exercise, to refine a definition, rather than give a rigorous definition at the outset.
If the student asks but what about 1, I can make 7 from 1*7, then you have an interaction, and have stimulated some thought.

Sunday, February 20, 2011

pari GP terminal colours - readable darkbg

After 3 years or more of using Pari/GP, the default highlighting (in terminal), bothered me enough, to read up how to change it.

default(colors,darkbg)

Here is a comparison of the effects of the highlighting change:


The lines prefixed %7 and %8 and in between, show the output is now in grey, and input highlighting is now greenish.

The original "1, 6, 3, 4, 5, 2, 3" profile known as lightbg is reactivated briefly, so it can be compared again.

Making the change permanent on Debian by:

emacs /etc/gprc

and changing which lines are commented out, results in the top of my /etc/gprc file looking like this:



So now when I start up Pari/GP my highlighting is set okay from the outset:




Notes and further reading:
On Ubuntu you will want to use the keyword sudo at the front of the command for editing /etc/gprc

If you want a lighter alternative to emacs on your system, then zile is a lightweight alternative that supports basic editing. Alternatively your system default editor is probably already chosen and use that.

The pari/GP command ...

default(colors,d)

... can be used interactively if you do not wish to make the change permanent.
darkbg abbreviated to just 'd' should work.

Friday, January 7, 2011

r-project / r-cran and install.packages()

My version of the excellent statistics and mathematics package from r-project.org, is now 2.11, having updated to Debian Squeeze.

( The major work from Google Summer of code 2010 is probably in 2.12 released December 2010 )

I tend to install just the minimal r-project packages and then immediately add rgl for pretty plot3d.

But before rushing into install.packages("rgl"), perhaps I should check where exactly the install code will end up.



looks to me as if /usr/local/lib/R/site-library might be the destination (set in R_LIBS_SITE), so let us see...

install.packages("rgl") pops up a country ordered list, for you to select a mirror for the install


Interesting to see the number of mirrors, 5 for Germany, and 5 for Canada, and I should scroll down to below Taiwan to pick my country.


( Above confirms my guesswork that any install.packages on my system defaults to /usr/local/lib/R/site-library )

Now what happens next, really changes the course of this article...so read on.

The install.packages() command failed as it could not locate GL/gl.h or GL/glu.h, and here I consulted my laptop.
The laptop already has rgl installed and working, so I used find and dpkg, and discovered my silly mistake.

There is already a debian package in squeeze named 'r-cran-rgl' which will do all the hard work for me :)

apt-get install r-cran-rgl 

will replace the manual install.packages() method and get you into plot3d() in no time at all :)

( for those who are interested in where GL/glu.h lives, try libglu1-mesa-dev in support of manual install of rgl )


The plot3d() function in rgl allows you 'drag to spin' functionality so you can rotate, and examine the points from any angle you choose.

In the above I have dragged the back round to the front so zero is away from you, and hopefully, you can just make out the points curving up towards point (5,5,60) which is the end of the curve.

Some further illustrations of what rgl can do can be found on statmethods.net in section 'Spinning 3D Scatterplots'.


Further reading and links:
For technical folks you really want more packages and are interested in install.packages(), here is the function definition with arguments:

Sunday, November 28, 2010

Gnumeric isprime() trial division

Gnumeric is one of several spreadsheet choices available to a Debian and / or GNU / Linux user.

How useful is Gnumeric to a mathematician? Quite.

It does include several functions that a more general package LibreOffice / OpenOffice Calc / Excel might not include.

Statistics and Number Theory are areas where Gnumeric currently outpaces the competition, in terms of available functions.

Here is one example, the isprime() function:
/*
 * Returns -1 (out of bounds), 0 (non-prime), or 1 (prime).
 */
static int
isprime (guint64 n)
{
      int i = 1;
      guint64 p = 2;

      if (n <= 1)
            return 0;

      for (i = 1; p * p <= n; i++) {
            if (ithprime (i, &p))
                  return -1;
            if (n % p == 0)
                  return 0;
      }

      return 1;
}

Do not be alarmed if you are not comfortable looking through the above code, it is simply there so that I may elaborate next.

The method here is trial division, which is a reasonable way of going about things for integers say of a billion or so. It is implemented in C and so is likely pretty fast at getting it's answer.

In fact a little bit of trial and error can show (very roughly) where the threshold for this isprime() function lies:

Firstly a candidate number (fairly large by spreadsheet standards):

factorial(18)+1=6402373705728001

which isprime() answers with #LIMIT!

Now those handy zeros at the end of the number will serve to reduce the number in size a little. Removing a zero from ...001 give ...01 and a new number:  640237370572801

...which Gnumeric isprime() is happy to process and return FALSE.

If you are wondering, then the new number ending ...2801 has a factor of 2801, how interesting, and completely accidental.

Anyway, I have provided a 15 digit number that isprime() will process, and demonstrated that it does have a limit to what it will attempt.

So no, you are not going to find huge primes by using entries in a Gnumeric spreadsheet, however it does provide very convenient environment, for exploring integer construction.

Should you use Gnumeric, to test out your hypothesis for prime number patterns?
Well why not, there is value in knowing that constructive forms hold true up to 15 digits certainly.

Later you will maybe explore a dedicated number theory program, there are several great open source offerings. However Gnumeric seems like a great place to just dive in.

Links and further reading: