Reading, Writing and Arithmetic

BSDNT - interlude

2011-10-23T08:23:00.000-07:00

You will notice that I have not worked on BSDNT for about a year now. Well, I'm thinking of restarting the project soon. I did complete two new revisions v0.25 and v0.26 since I stopped blogging. The first of these added random functions for a single word. They generate single words which have an increased probability of triggering corner cases, e.g. by a sparse binary representation. The second of these updates is a bsdnt_printf function. This is like printf but adds a %w format specifier for printing single words. There is also a %m for printing a len_t and %b for printing a bits_t.

I likely won't get much more done on BSDNT until early next year, but here is what I am planning:

1) I am tremendously grateful to Dr. Brian Gladman for his work on an MSVC version of the library. However, I started to struggle to keep up with this side of things more than I thought. Microsoft's MSVC doesn't support inline assembler in 64 bit x86. This means the entire plan of the MSVC version has to be different. It seems like far too much effort to combine both sets of code into a single library. I've therefore decided (sorry Brian) to ditch the MSVC code from my copy of BSDNT.

2) I wasn't happy with the interface of the division code. There are a few issues to consider.

The first issue is chaining. Obviously carry-in occurs at the left rather than the right. But for general division functions should the carry-in be a single limb or multiple limbs. It seems like the remainder after division is going to be m limbs and so the carry-in should be also. It is not clear what is better here. Internally, the algorithms deal with just a single carry-in limb because they use 2x1 divisions to compute the quotient digits. Perhaps chaining just means that we consider the first digit of the remainder to be carry-in for the next division.

Another problem associated with this is that when reading the carry-in from the array, if the carry-in happens to be zero then the array entry may not exist in memory. This means the code has to always check if the carry-in should be zero or not before proceeding.

The other issue to consider is shifting for normalisation. One assumes that the precomputed inverse is computed from a normalised value (the first couple of limbs of the divisor). Now, it is not necessary to shift either dividend or divisor ahead of time. One can still perform the subtractions that occur in division, on unshifted values. One does need to shift limbs of the dividend in turn however, as the algorithm proceeds, in order to compute the limbs of the quotient. But this shifting can occur in registers and need not be written out anywhere. This is implemented, but currently every limb gets shifted twice. This can be cut down to a single shift.

3) The PRNGs are currently quite hard to read. They have numerous macros to access their context objects. They are extremely flexible, but possibly overengineered. I'd like to simplify their implementations somewhat.

4) The configure script is a little overengineered. The idea of supporting lots of compilers is nice. But in reality GCC should exist almost everywhere. The original concept of BSDNT was to use inline assembly for architecture support. This gets around issues with global symbol prefixes and wotnot. It also makes the library really simple to read. Even on Windows 64 there is MinGW64 and this is the only setup I aim to target in that direction.

I hope to deal with all of these issues before proceeding with development of BSDNT. Give me some time as I am busy until about the end of the year. However, I do plan to continue development of BSDNT after sorting out these issues, because I think that fundamentally what we have is very solid.

Previous article: v0.24 = nn_bitset/clear/test and nn_test_random

BSDNT - v0.24 nn_bitset/clear/test and nn_test_random

2010-11-20T15:01:00.000-08:00

In today's update we make a long overdue change to bsdnt, again to improve our testing

of the library.

We are going to add a function for generating random bignums with sparse binary

representation. We'll also add some other random functions based on this primitive.

Using test integers with sparse binary representation in our test code will push our

functions harder because it will test lots of corner cases such as words that are all

zero, in the middle of routines, and so on. As it is currently, we'd be extremely

lucky for the random word generator we've been using to generate an all zero word, or

a word with all bits set to one for that matter.

The first step to generating such test randoms is to write a function for setting a

given bit in an integer. This will be an nn_linear function despite it not actually

taking linear time. In fact, it will take essentially constant time. However, it is an

nn type function, so it belongs in an nn module.

The routine is very straightforward. Given a bit index b, starting from 0 for the least

significant bit of a bignum, it will simply use a logical OR to set bit b of the bignum.

Rather than construct a bignum 2^b and OR that with our number, we simply determine

which word of the bignum needs altering and create an OR-mask for that word.

Computing which word to adjust is trivial, but depends on the number of bits in a word.

On a 64 bit machine we shift b to the right by 6 bits (as 2^6 = 64), and on a 32 bit

machine we shift b to the right by 5 bits (2^5 = 32). This has the effect of dividing

b by 64 or 32 respectively (discarding the remainder). This gives us the index of the

word we need to adjust.

Now we need to determine which bit of the word needs setting. This is given by the

remainder after dividing b by 64 or 32 respectively, and this can be determined by

logical AND'ing b with 2^6-1 or 2^5-1 respectively. This yields a value c between 0 and

63 (or 31) inclusive, which is a bit index. To turn that into our OR-mask, we simply

compute 2^c (by shifting 1 to the left by c bits).

Now that we have our OR-mask and the index of the word to OR it with, we can update the

required bit. We call this function nn_bit_set.

While we are at it we create two other functions, nn_bit_clear and nn_bit_test.

It's now straightforward to write test functions which randomly set, clear and test

bits in a random bignum.

Next we create a random bignum generator which sets random bits of a bignum. In order

to do this, we simply choose a random number of bits to set, from 0 to the number of words

in the bignum, then we set that many bits at random.

We call this function nn_test_random1.

We now add a second random bignum generator. It works by creating two bignums using the

function nn_test_random1 and subtracting one from the other. This results in test randoms

with long strings of 1's and 0's in its representation.

We call this function nn_test_random2.

Finally, we create a function nn_test_random which randomly switches between these two

algorithms and our original nn_random algorithm to generate random bignums. We switch all

our test code to use nn_test_random by changing the function randoms_of_len to use it.

At this point we can have greater confidence that our functions are all working as they

are supposed to be, as our test code has been suitably fortified at last! (Well, they are

working now, after I spent a day hunting down the bugs that these new randoms found - no,

I am not kidding. That's how good at finding bugs this trick is!)

Today's code is found here: v0.24

Previous article: v0.23 - sha1 and PRNG tests

Next article: BSDNT - interlude

BSDNT - v0.23 sha1 and PRNG tests

2010-11-12T07:30:00.000-08:00

In a recent update we added three PRNGs (pseudo random number

generators). What we are going to do today is add test code for

these.

But how do you test a pseudo random generator? It's producing

basically random values after all. So what does it matter if the

output is screwed up!?

Well, it does matter, as shown by the problems on 32 bit machines

which I wrote about in the PRNG blog post. It would also matter if

the PRNGs were broken on some platform such that they always output

0 every time!

There's a few ways of testing PRNGs. One is simply to run them for a

given number of iterations, write down the last value it produces and

check that it always does this.

The method we are going to use is slightly more sophisticated. We are

going to hash a long series of outputs from the PRNGs, using a hash

function, and check that the hash of the output is always the same.

Basically, our test code will take a long string of words from the

PRNGs, write them to an array of bytes, then compute the sha1 hash of

that array of bytes. It'll then compare the computed hash with a hash

we've computed previously to ensure it has the same value as always.

Moreover, we'll set it up so that the hash is the same regardless of

whether the machine is big or little endian.

The hash function we are going to use is called sha1. Specifically,

we'll be using an implementation of the same written by Brian Gladman

(he also supplied the new test code for the PRNGs for today's update).

The first step is to identify whether the machine is big or little

endian. This refers to the order of bytes within a word in physical

memory. On little endian machines (such as x86 machines) the least

significant byte of a word comes first. On big endian machines the

order is the other way around. Thus the number 0x0A0B0C0D would have

the byte 0x0D stored first on a little endian machine, but 0x0A stored

first on a big endian machine.

Endianness can be identified by architecture, or it can be identified

with a short program. We choose to use the latter method as it should

be hard to fool. At configure time a short C program will run that will

place bytes into a four byte array, then read that array as a single

32 bit number. We then compare the value to a 32 bit value that would

be stored in the given way on a little endian machine. If it compares

equal, then the machine must be little endian. If not we compare with

a number that would be stored in the given way on a big endian machine.

If the machine doesn't match either order, then it must be a very rare

machine with mixed endianness, which we don't support in bsdnt.

The configure script will write some defines out to config.h which

then allow bsdnt modules to identify whether the machine is little or

big endian at compile time, i.e. at zero runtime cost.

Now to discuss the sha1 hashing scheme.

A hashing scheme simply takes a piece of data and computes from it a

series of bits which serve to "identify" that piece of data. If

someone else has access to the same hashing algorithm and a piece of

data which purports to be an exact copy of the original, then they

can verify its identity by computing its hash and comparing.

Of course a hash is only valuable in this sense if it is much shorter

than the piece of data itself (otherwise you'd just compare the

actual data itself).

A very simple hashing scheme might simply add all the words in the

input to compute a hash consisting of a single word.

A secure hashing scheme has at least two other properties.

i) It shouldn't be possible to determine the original data from its

hash. (The data may be secret and one may wish to provide for the

independent verification of its authenticity by having the recipient

compare the hash of the secret data with a publicly published value.

Or, as is sometimes the case, the hash of secret data, such as a

password, might be transmitted publicly, to compare it with a

pre-recorded hash of the data.)

ii) It must be computationally infeasible to substitute fake data

for the original such that the computed hash of the fake data is the

same as that of the original data.

Of course, if the hash is short compared to the data being hashed,

then by definition many other pieces of data will have the same hash.

The only requirement is that it should be computationally infeasible

to find or construct such a piece of data.

The first step in the SHA1 algorithm is some bookkeeping. The

algorithm, as originally described, works with an input message which

is a multiple of 16 words in length. Moreover, the last 64 bits are

reserved for a value which gives the length of the original message in

bits.

In order to facilitate this, the original message is padded to a

multiple of 16 words in length, with at least enough padding to allow

the final 64 bits to be part of the padding, and to not overlap part

of the message.

The padding is done by first appending a single binary 1 bit, then

binary zeroes are appended until the message is the right length.

Then of course the length in bits of the original message is placed

in the final 64 bits of the message.

The hashing algorithm itself performs a whole load of prescribed

twists and massages of the padded message.

For this purpose some functions and constants are used.

Given 32 bit words B, C and D there are four functions:

1) f0(B, C, D) = (B AND C) OR ((NOT B) AND D)

2) f1(B, C, D) = B XOR C XOR D

3) f2(B, C, D) = (B AND C) OR (B AND D) OR (C AND D)

4) f3(B, C, D) = B XOR C XOR D

and four corresponding 32 bit constants:

1) C0 = 0x5A827999

2) C1 = 0x6ED9EBA1

3) C2 = 0x8F1BBCDC

4) C3 = 0xCA62C1D6

To begin the algorithm, we break the padded message up into 16 word

blocks M1, M2, M3, i.e. each Mi is 16 words of the padded message.

Each block is processed via a set of steps, and an "accumulated" hash

of 160 bits, consisting of five 32 bit words (the final hash we are

after) is computed: H0, H1, H2, H3, H4.

The algorithm starts by initialising the five "hash words" to the

following values:

H0 = 0x67452301, H1 = 0xEFCDAB89, H2 = 0x98BADCFE, H3 = 0x10325476

and H4 = 0xC3D2E1F0.

Each block of 16 words, Mi, of the padded message is then used in

turn to successively transform these five words, according to the

following steps:

a) Break Mi into 16 words W0, W1, ..., W15.

b) For j = 16 to 79, let Wj be the word given by

Wj = S^(-1)(W{j-3}) XOR W{j-8} XOR W{j-14} XOR W{j-16}),

where S^n(X) means to rotate the word X to the left through n bits

(a negative n means right rotation).

c) Make a copy of the hashing words:

A = H0, B = H1, C = H2, D = H3, E = H4

d) For j = 0 to 79 repeat the following set of transformations in

the order given:

TEMP = S^5(A) + f{j/20}(B, C, D) + E + Wj + C{j/20},

E = D, D = C, C = S^30(B), B = A, A = TEMP,

where j/20 signifies "floor division" by 20, and where f and C are

the above-defined functions and constants.

e) Update the hashing words according to the following:

H0 = H0 + A, H1 = H1 + B, H2 = H2 + C, H3 = H3 + D, H4 = H4 + E.

Note that steps a-e are repeated for each block of 16 words, Mi in

the padded message, further manipulating the five words with each run.

The resulting five words H0, H1, H2, H3, H4 after all the words of the

padded message have been consumed, constitutes the sha1 hash of the

original message.

The function to compute the sha1 hash of a message is given in the

files sha1.c and sha1.h in the top level directory of the source

tree.

A new test file t-rand.c is added in the test directory. It contains

the sha1 hash of a large number of words as output by our three

PRNGs. If a user of bsdnt has the same hash for the PRNGs when run

on their machine, then they can have a very high level of confidence

that they are working as expected.

Note that the sha1 algorithm is known as a secure hashing algorithm,

which means that in theory it can be used to hash very important

data so that the recipient can independently confirm the data hasn't

been tampered with (by computing the hash of the value and making

sure it matches some published value).

We don't explain how sha1 actually works. The mysterious constants

are not so mysterious. C0 is the square root of 2 in hexadecimal, C1 is

the square root of 3, C2 the square root of 5, C3 the square root of 10.

I don't know the meaning of the functions f0-f3.

What is worth noting is that in recent years, people have figured out

how to produce sha1 hash collisions (two messages with the same hash).

I don't pretend to be an expert in such things.

All we care about here is that a broken PRNG really can't pretend to

be working, and for that, sha1 works a treat.

Disclaimer: use the information in this post at your OWN RISK!! We

make no representations as to its correctness. The same goes for

bsdnt itself. Read the license agreement for details.

Warning: cryptography is restricted by law in many countries including

many of those where the citizens believe it couldn't possibly be so.

Please check your local laws before making assumptions about what you

may do with crypto.

The code for today's article is here: v0.23

Previous article: v0.22 - Windows support

Next article: v0.24 - nn_bitset/clear/test and nn_test_random

BSDNT - v0.22 Windows support

2010-11-11T05:58:00.000-08:00

Today's update is a massive one, and comes courtesy of Brian Gladman. At last we add

support for MSVC 2010 on Windows.

In order to support different architectures we add architecture specific files in the arch

directory. There are three different ways that architectures might be supported:

* Inline assembly code

* Standalone assembly code (using an external assembler, e.g. nasm)

* Architecture specific C code

Windows requires both of the last two of these. On Windows 64 bit, MSVC does not support

inline assembly code, and thus it is necessary to supply standalone assembly code for this

architecture. This new assembly code now lives in the arch/asm directory.

On both Windows 32 and 64 bit there is also a need to override some of the C code from the base

bsdnt library with Windows specific code. This lives in the arch directory.

Finally, the inline assembly used by the base bsdnt library on most *nix platforms is now in the

arch/inline directory.

In each case, the os/abi combination is specified in the filenames of the relevant files. For

example on Windows 32, the files overriding code in nn_linear.c/h are in arch/nn_linear_win32.c/h.

(Note win32 and x64 are standard Windows names for 32 and 64 bit x86 architectures, respectively.)

If the code contains architecture specific code (e.g. assembly code) then the name of the file

contains an architecture specifier too, e.g. arch/inline/nn_linear_x86_64_k8.h for code specific

to the AMD k8 and above.

It's incomprehensible that Microsoft doesn't support inline assembly in their 64 bit compiler

making standalone assembly code necessary. It would be possible to use the Intel C compiler on

Windows 64, as this does support inline assembly. But this is very expensive for our volunteer

developers, so we are not currently supporting this. Thus, on Windows 64, the standalone

assembly is provided in the arch/asm directory as just mentioned.

Brian has also provided MSVC build solution files for Windows. These are in the top level source

directory as one might expect.

There are lots of differences on Windows that requires functions in our standard nn_linear.c,

nn_quadratic.c and helper.c files to be overridden on Windows.

The first difference is that on 64 bit Windows, the 64 bit type is a long long, not a long. This

is handled by #including a types_arch.h file in helper.h. On most platforms this file is empty.

However, on Windows it links to an architecture specific types.h file which contains the

requisite type definitions. So a word_t is a long long on Windows.

Also, when dealing with integer constants, we'd use constants like 123456789L when the word type

is a long, but it has to become 123456789LL when it is a long long, as on Windows 64. To get

around this, an architecture specific version of the macro WORD(x) can be defined. Thus, when

using a constant in the code, one merely writes WORD(123456789) and the macro adds the correct

ending to the number depending on what a word_t actually is.

Some other things are different on Windows too. The intrinsic for counting leading zeroes is

different to that used by gcc on other platforms. The same goes for the function for counting

trailing zeroes. We've made these into macros and given them the names high_zero_bits and

low_zero_bits respectively. The default definitions are overridden on Windows in the architecture

specific versions of helper.h in the arch directory.

Finally, on Windows 64, there is no suitable native type for a dword_t. The maximum sized

native type is 64 bits. Much of the nn_linear, and some of the nn_quadratic C code needs to

be overridden to get around this on Windows. We'll only be using dword_t in basecase algorithms

in bsdnt, so this won't propagate throughout the entire library. But it is necessary to

override functions which use dword_t on Windows.

Actually, if C++ is used, one can define a class called dword_t and much of the code can

remain unchanged. Brian has a C++ branch of bsdnt which does this. But for now we have C code

only on Windows (otherwise handling of name mangling in interfacing C++ and assembly code

becomes complex to deal with).

Brian has worked around this problem by defining special mul_64_by_64 and div_128_by_64

functions on 64 bit Windows. These are again defined in the architecture specific version of

helper.h for Windows 64.

Obviously some of the precomputed inverse macros need to be overridden to accomodate these

changes, and so these too have architecture specific versions in the Windows 64 specific version

of the helper.h file.

In today's release we also have a brand new configure file for *nix. This is modified to handle

all the changes we've made to make Windows support easy. But Antony Vennard has also done

some really extensive work on this in preparation for handling standalone assembly on arches

that won't handle our inline assembly (and for people who prefer to write standalone assembly

instead of inline assembly).

The new configure file has an interactive mode which searches for available C compilers (e.g.

gcc, clang, icc, nvcc) and assemblers (nasm, yasm) and allows the user to specify which to use.

This interactive feature is off by default and is only a skeleton at present (it doesn't actually

do anything). It will be the subject of a blog post later on when the configure file is finished.

The code for today is found at: v0.22

Previous article: v0.21 - PRNGs

Next article: v0.23 - sha1 and prng tests

BSDNT - v0.21 PRNGs

2010-10-31T16:02:00.000-07:00

In this update we are going to replace the el cheapo random generator in
bsdnt with something of higher quality.

Some time ago, Brian Gladman brought to my attention the fact that on 32
bit machines, the test code for bsdnt actually caused Windows to hang!

The issue required some sleuthing work on Brian's part to track down.
He eventually discovered the cause of the problem, and it was, oddly
enough, the pseudo-random number generator (PRNG) I had used.

Brian suspected the PRNG immediately because of his past experience as a
professional cryptographer. In fact, it turns out that PRNG's of the
kind that I had used, aren't particularly good even if they don't have
bugs!

The particular kind of PRNG I had used is called a linear congruential
PRNG. They start with the PRNG initialised to some random seed value,
n = seed. Then each time they are called, they apply the transformation
n = (n*c1 + c2) % p for some explicit constants c1, c2 and a large
enough "prime" p.

One can take c2 = 0 in the transformation and it is also common to see
p = 2^b for some b (e.g. b = WORD_BITS, and yes, I know p = 2^b is not
usually prime). When p = 2^b it is usually the case that the top half of
the bits output have reasonable random properties, but the bottom half
do not. In this case it is acceptable to stitch two LC PRNG's together to
give the full number of bits.

When p is an actual prime, these PRNG's are called prime modulus linear
congruential PRNG's and they aren't too bad when implemented properly.
They still fail to meet some standards of random quality.

To return a whole word of random bits, one either needs to use a prime
p that is larger than a word, which is usually impractical, or again
stitch the output of two prime modulus LC PRNG's together.

However, one needs to be careful, as the period of the generator is p-1,
so if one is on a 32 bit machine, it doesn't do to use a prime p just
over half the size of a word (the first mistake I made), otherwise the
period is just over 65536. That isn't too good for generating random
values for test code!

But how was my LC PRNG causing Windows to crash!? The reason was that
in some of the test functions we required bignums that didn't overflow
when summed together. This of course depends almost entirely on the top
few bits of the bignums being added together.

The problem was that in the expression n = (n*c1 + c2) % p, I was using
values of c1 and c2 which were not reduced mod p. It turns out that this
is crucial to correct operation. It might seem that the result ends up
being reduced mod p anyway, and indeed it would be if n*c1 fit in a word.
However, because it doesn't you actually get ((n*c1 + c2) % 2^32) % p
which causes the binary representation of the value generated to be quite
regular.

Anyhow, on Windows (and probably on other 32 bit machines) the test code
generates length 90 bignums over and over at some point, looking in vain
to find pairs of such bignums which when added do not overflow. As these
are garbage collected at the end of the test function, memory starts
filling up with the orphaned bignums that are discarded by the test code
as it looks for appropriate values. This eventually overwhelms the heap
allocator on Windows and causes the entire OS to crash!

The problem of writing decent PRNG's has been studied extensively, and one
expert in the subject is George Marsaglia. He famously turned up on a
usenet forum in January of 1999 and dumped not one, but piles of fast, high
quality PRNG's which do not suffer from the problems that other PRNG's do.

Amazingly, many of the PRNG's in common usage today are either written by
George, or based on ones he wrote. So he's some kind of legend!

Anyhow, we will make use of two of his PRNG's, Keep It Simple Stupid (KISS)
and Super KISS (SKISS) and a third PRNG called Mersenne Twister, due to
Makoto Matsumoto and Takuji Nishimura in 1997.

George's generators are in turn based on some simpler PRNG's. He begins by
defining a linear congruential generator, with c1 = 69069 and c2 = 1234567.
This is taken p = mod 2^32 (yes, it's not prime). This has good properties
on its top 16 bits, but not on its bottom 16 bits, and for this reason had
been widely used before George came along. This generator has period 2^32.

Next he defines a pair of multiply with carry (MWC) generators. These are
of the form n = c1*lo(n) + hi(n) where lo(n) is the low 16 bits of n, hi(n)
is the high 16 bits and c1 is an appropriately chosen constant.

He stitches together a pair of these MWC PRNG's mod 2^16 to give 32 random
bits. For simplicity we'll refer to this combined random generator as MWC.
This has a period of about 2^60.

Thirdly he defines a (3-)shift-register generator (SHR3). This views the
value n as a binary vector of 32 bits and applies linear transformations
generated from 32 x 32 matrices L and R = L^T according to
n = n(I + L^17)(I + R^13)(I + L^5) where I is the 32 x 32 identity matrix.

In order to speed things up, special transformations are chosen that can
be efficiently implemented in terms of XOR and shifts. This is called an
Xorshift PRNG. We'll just refer to it as SHR3.

Now given appropriate seed values for each of these PRNG's Marsaglia's
KISS PRNG is defined as MWC ^ CONG + SHR3. This generator passes a whole
slew of tests and has a period of 2^123. In this update we make it the
default random generator for bsdnt.

Super KISS is a random generator defined by George later in 2009. It gives
immense periods by adding together the output of three PRNG's, one with a
massive order. It is basically defined by SKISS = SUPR + CONG + SHR3.

Here, the new generator SUPR is based on a prime p = a*b^r + 1 such that
the order of b mod p has magnitude quite near to p - 1.

It starts with a seeded vector z of length r, all of whose entries are
less than b and an additional value c which is less than a.

One then updates the pair (z, c) by shifting the vector z to the left by
one place and setting the right-most entry to (b - 1) - ((a*z1 + c) mod b)
where z1 is the entry shifted out at the left of z. Then c is set to t/b.

Naturally in practice one uses b = 2^32 so that all the intermediate
reductions mod b are trivial.

As with most generators which have massive periods the "state" held by this
generator is large. It requires data mod p for a multiprecision p.

Note the similarity with the MWC generator except for the "complement" mod
b that occurs. This is called a CMWC (Complemented-Multiply-With-Carry)
generator.

George proposed using the prime p = 640*b^41265+1, where the order of b is
5*2^1320481. The period of the CMWC generator is then greater than
2^1300000.

Of course, at each iteration of the algorithm, 41265 random words are
generated in the vector. Once these are exhausted, the next iteration of
the algorithm is made.

The algorithm SUPR in the definition of SKISS is thus just a simple
array lookup to return one of the words of the vector z. Each time SKISS
is run, the index into the array is increased until all words of the array
are exhausted, at which point the CMWC algorithm is iterated to refill the
array.

We now come to describing the Mersenne twister.

It is based on the concept of a feedback shift register (FSR). An FSR shifts
its value left by 1 bit, feeding at the right some linear combination of the
bits in its original value. The Mersenne twister is conceptually a
generalisation of this.

The difference with the Mersenne twister is that the "feedback" is effected
by a certain "twist". This is effected by applying a "linear transformation"
A of a certain specific form, with multiplication by A having addition
replaced by xor in the matrix multiplication. The twist can be described
more straightforwardly, and we give the more straightforward description
below.

One sets up a Mersenne twister by picking a recurrence degree n, a "middle
word" 1 <= m <= n and a number of bits for a bitmask, 0 <= r <= 32. One

picks these values so that p = 2^(n*w - r) - 1 is a Mersenne prime (hence

the name of this PRNG). Given a vector of bits a = [a0, a1, ..., a{w-1}] of length

w and a sequence x of words of w bits, the Mersenne twister is defined by a

recurrence relation x[k+n] = x[k+m] ^ ((upper(x[k]) | lower(x[k+1])) A)

where upper and lower return the upper w - r and lower r bits of their

operands, and where A is the "twist" spoken of and defined below, in terms of

a. Of course ^ here is the xor operator, not exponentiation. For a vector X of w

bits, XA is given by X>>1 if X[0] == 0 otherwise it is given by (X>>1) ^ a.

Some theory is required to find an A such that the Mersenne twister will have
maximum theoretical period 2^(n*w - r) - 1.

To finish off, the Mersenne twister is usually "tempered". This tempering
simply mangles the bits in a well understood way to iron out some of the
known wrinkles in the MT algorithm.

Only a couple of sets of parameters are in common use for Mersenne twisters.
These are referred to as MT19937 for 32 bit words and MT19937-64 for 64 bit
words.

As with all PRNG's, there is a whole industry around "cracking" these things.
This involves starting with a short sequence of values from a PRNG and
attempting to find the starting constants and seed values.

Obviously, in crytographic applications, there is not much point generating
"secure" keys with a PRNG with a single source of entropy. Even if your key
is generated by multiplying primes of many words in length, if those words
were generated from a PRNG seeded from the current time, it may only take
a few iterations and a good guess as to which PRNG you used, to determine
the constants used in the PRNG and thus your entire key. And that's
irrespective of which constants you chose in your PRNG!

So if you are doing crypto, you need to take additional precautions to
generate secure keys. Just seeding a PRNG from the time probably won't cut
it!

Some PRNG's are more "secure" than others, meaning that knowing a
few output values in a row doesn't give terribly much information about
which values may follow. But if you rely on a PRNG to be secure, you
are essentially betting that because you don't know how to get the
next few values and nor does anyone else that has written about the
subject, then no one at all knows. Of course one needs to ask oneself
if they would tell you if they did.

Another assumption one should never make is that no one has the computing
power to brute force your PRNG.

Some PRNG's are designed for cryptographic applications, and maybe one can
believe that these are "safe" to use, for some definition of safe.

Anyhow, we only care about random testing at this point. In today's update
32 and 64 bit KISS, SKISS and MT PRNG's are added in the directory rand.
Our randword, randinit, and randclear functions are all replaced with
appropriate calls to KISS functions.

There is also an option to change the default PRNG used by bsdnt. Is it my
imagination or does the test code now run faster, even on a 64 bit machine!

At some point we will add some tests of the new PRNG's. These will compare
the outputs with known or published values to check that they are working as
designed for a large number of iterations.

Brian Gladman contributed to this article and also did most of the work
in implementing Marsaglia's PRNG's in bsdnt. The Mersenne twisters were
originally written by Takuji Nishimura and Makoto Matsumoto and made available
under a BSD license. Brian did most of the work in adapting these for bsdnt.

The code for today's update is here: v0.21

Previous article: v0.20 - redzones

Next article: v0.22 - Windows support

BSDNT - v0.20 redzones

2010-10-30T12:12:00.000-07:00

In this update we implement another improvement to the test code in bsdnt. I don't know

what the correct name is, but I call them redzones.

The basic idea is this: suppose you have a function nn_blah say, and it writes to an nn_t

b say. If it writes well beyond the allocated space for b, then almost certainly a

segfault will occur. But what if it only writes a word or two before the beginning or

after the end of the allocated space? Very likely this will cause a segfault only on

some systems, depending on the granularity of the heap allocator and depending

on what other bsdnt data might be in the overwritten space!

So what if we could detect this kind of error? Well, that is what redzones hope to do.

Essentially if an nn_t b is allocated with m words of space, when redzones are turned on

it allocates m + 2C words of space for some small constant C. It then fills the first

and last C words of b with known words of data (usually some recognisable pattern of bits).

When the garbage collector cleans up, it examines the redzones to ensure that they have

not been altered. If they have, they raise an error.

The nn_t b is set to point just after the first C words, which contain the redzone, and in

every other respect act like a normal nn_t. The user needn't know that an extra C words

of data were allocated immediately before and after the length m nn_t they requested.

Nor do they need to be aware of the checking that goes on when the nn_t is finally cleaned

up, that the redzones haven't been touched.

Of course it's nice to be able to turn redzones off sometimes, when testing the library.

Therefore I've added a configure option -noredzones which turns off redzones if they are

not required. This works by setting a #define WANT_REDZONES 0 in config.h. The

memory allocator for nn_t's and the garbage collector both operate differently if redzones

are turned on.

At present, the only way to allocate memory for nn_t's in test code is to use

randoms_of_len, so it is convenient to rewrite this to call a function alloc_redzoned_nn

instead of malloc, and for the garbage collector to call free_redzoned_nn. These new

functions are defined in test.c.

The only difference when WANT_REDZONES is set in config.h is that REDZONE_WORDS, which is defined in test.h is changed from 0 to 4 words (meaning 4 redzone words are to be

allocated at each end of a redzoned nn_t). Having redzones of length 0 is the same as not

having them at all. So this makes the functions easy to write.

Also in test.h REDZONE_BYTE is defined to the hexaecimal byte 0xA which has binary bit

pattern 1010, i.e. alternating one's and zeroes. This is the value that is placed into the

redzones byte-by-byte before the nn_t is used. At the end, when they are cleaned up, the

garbage collector examines the redzones to ensure they are still filled with these bytes.

Fortunately checking redzones does not dramatically slow down our test code, and no new

test failures result. This means it is highly likely that our nn_t functions do not overwrite

their bounds. To check that the new redzones code works, it is a simple matter of mocking

up a broken function which overwrites its bounds. The new code complains loudly as it

should, unless redzones are switched off at configure time.

The code for today's update is here: v0.20

Previous article: v0.19 - asserts

Next article: v0.21 - prngs

BSDNT - v0.19 asserts

2010-10-25T19:38:00.000-07:00

About a week ago I got enthused to work on another coding project I've been

wanting to experiment with for a long time. I discovered that it was highly

addictive and I just couldn't put it down. It's also given me some interesting

ideas for a higher level interface to bsdnt. But more on that later when we

start working on it.

Unfortunately in that week of time there have been no bsdnt updates.

Moreover, on the weekend my main computer died (physical hard drive failure).

I pulled out my backup machine and Windows wanted to install 3 months of

"important updates". Naturally this caused the machine to crash, I was unable

to recover to the restore point it set, the startup repair didn't work and

the only solution was a format and reload.

*Fifteen and a half hours* later I had reinstalled Windows and it had finally

finished putting the approximately 165 "important updates" back on!!

Unfortunately all the bsdnt blog articles I had meticulously prepared in

advance were lost. Thus I am regenerating what I can from the diff between

revisions of bsdnt. Sorry if they end up being shorter than past updates.

Fortunately I did not lose any of the code I wrote, as that was backed up in

a git repo on an external server!

Anyhow, in this update we make a very simple change to bsdnt, again in an

attempt to improve the test quality of the library. We add asserts to the

code.

An assert is a check that is made at runtime in live code to test if some

predefined condition holds. If the assert fails, an error message is printed

specifying the line of code where the assert is located and what the

condition was that failed.

Now, I am not personally a great believer in asserts. As they are runtime

checks, they require computing cycles, which is just a no-no for a bignum

library. The other option is to turn them off when not testing code. However,

this simply leads to the asserts rarely being run when they are needed.

The other problem with asserts is that they pollute the code, making the

source files longer and appear more complex.

However, there is one situation where I believe they can be very helpful,

and that is in checking the interface of functions within a library and that

it is being respected both in intra-library calls and by the test code for

the library.

Specifically, assert are useful for checking that valid inputs have been

passed to the functions, e.g. you might have a restriction that a Hensel

modulus be odd. Adding an assert allows you to test that all the moduli

you pass to the function in your test runs are in fact odd.

The main advantage in putting asserts into the code is that it forces you

to think through what all the conditions should be that you assert. In

adding asserts to the code in bsdnt I discovered one function in which the

test code was pushing the code to do things I didn't write it to cover.

This forced me to either rewrite the test, or drop that as a condition (I

think I chose the former for consistency with other related functions in

bsdnt).

Of course we do not want to consume cycles when the library is run by the

end user, and so we make asserts optional. This is done using a configure

switch. By default the macro WANT_ASSERT is set to 0 in a file config.h by

configure. However, if the user passes the option -assert to configure, it

sets the value of this define to 1.

A macro ASSERT is then defined in helper.h which is either an empty macro

in the default case or is an alias for the C assert function if WANT_ASSERT

is set to 1 in config.h.

Of course we have to remember to turn asserts on to run the test code, and

this really highlights their main weakness. As I mentioned, the asserts I

added did clarify the interface, but I don't believe they showed up any

bugs in bsdnt. With this expectation, asserts can be a useful tool.

The code for today's update is here: v0.19

Previous article: v0.18 - printx_word, nn_printx

Next article: v0.20 - redzones

BSDNT - v0.18 printx_word, nn_printx

2010-10-16T14:51:00.000-07:00

It is time we improved our test code again. We'll spend a few days updating

things to make improvements in the way we test.

Today's update is quite straightforward. We currently have no way of printing

nn_t's. This is quite inconvenient when it comes to the test code, where

little to no diagnostic information is printed at all. In particular, we

aren't printing out any of the multiple precision integers for examination

when a test fails.

Now, it is actually quite a difficult job to print bignum integers in decimal.

In fact, as far as I can see, one requires a function which allocates temporary

space to efficiently print integers. This is an interesting challenge:

is there an algorithm to convert from binary to decimal and print the result,

with just O(1) temporary space, with any complexity.

I think it may be possible if one allows the input to be destroyed. If so, a

subsidiary question would be to do the same thing without destroying the

input. I doubt that is possible, but I do not have a proof. Of course, to be

practical, we'd require an algorithm which doesn't destroy the input.

To get around this issue, we'll start with a simple nn_printx algorithm,

which will print a bignum in hexadecimal. We also add an nn_printx_short

function which prints the first couple of words of a bignum, an ellipsis and

then the final couple of words. This is useful for large bignums that would

print for screens and screens due to their size. We'll use this in our test

code to prevent printing too much output upon test failure.

Another function we add is an nn_printx_diff function. It accepts two nn_t's

and prints information about the range of words where they differ and prints

the first and last differing word in each case.

There is one tricky aspect to our printing functions however. A word is often

an unsigned long, but on some platforms it will be an unsigned long long. For

this reason, when printing a word, we need to use %lx as the format specifier

on some platforms and %llx on others.

So we need to add a routine which will print a word and abstract away the

format specifier so the caller doesn't have to think about it. The function

we include to do this is caled printx_word. It prints a word without needing

to specify a format specifier.

We add files helper.c/h to bsdnt which will contain routines like this one

which aren't specific to our nn module. A few existing functions and macros

also get moved there. The configure system will automatically look for

architecture specific versions of helper.c, allowing us to override the

definition of the functions in that file on a per architecture basis.

We add the printx_word function to helper.c which can be overridden with

an architecture specific version. On a platform where %llx is required, an

architecture specific version will simply replace the generic version which

uses %lx.

In test.h we add some macros, print_debug and print_debug_diff which use the

stringizing operator to print the names of the variables and then print their

values. The stringizing operator (#) is a preprocessor macro which turns a

macro parameter into a string. In our case, we pass the variable name to the

macro and turn it into a string so that we can print the variable name.

A few modifications to the TEST_START and TEST_END macros in test.h also

allow us to give a unique name to each test which is then printed along with

the iteration at which the test failed. This also uses the stringizing

operator so that the caller of TEST_START can specify the unique name for the

test. It seems difficult to come up with an automatic way of generating

unique test names, so this will have to do.

It would also be a useful thing to have it print the value of the random seed

at the start of a failing iteration too. After we have improved the random

generation code in bsdnt v0.21, perhaps someone would like to try adding this

feature.

We roll out our new diagnostic printing routines to our test code. Of course,

to see any of this new code in action, one has to introduce a bug in one of

the tests so that the new diagnostic code is actually run. I leave it to you

to fiddle around introducing bugs to see that the new test code does actually

print useful diagnostic information.

Later on we'll add a bsdnt_printf function which will be variadic and accept

a format specifier like the C printf function and which will have a

consistent %w for printing a word. This will also make things easier on

Windows, where currently the format specifier will be wrong in many places.

We'll fix this problem in a later update.

The code for today's post is here: v0.18

Previous article: v0.17 - div_hensel

Next article: v0.19 - asserts

BSDNT - v0.17 div_hensel

2010-10-15T17:43:00.000-07:00

Now that we have nn_mullow_classical, we can add nn_div_hensel. As explained,

this will take an integer a of length n and divide by an integer d of length

m modulo B^n, returning a quotient q and an "overflow".

The overflow will be two words which agree with the overflow from mullow(q*d).

As per the euclidean division, the dividend a will be destroyed.

The algorithm is somewhat simpler than the euclidean algorithm. If d1 is the

least significant word of d then we use an inverse mod B of d1 (dinv say) and

multiply it by the current word of the dividend being considered (working from

right to left) to get a quotient word q1.

We then subtract q1*d (appropriately shifted) from the dividend. There is no

adjustment to do as the inverse mod B is unique (so long as d is odd).

Any borrows and overflows from the subtractions are accumulated in the two

overflow and returned.

In our test code, we check a few things. Firstly, for an exact division, we

want that the quotient is really the exact quotient of a by d. As the

quotient returned is m words, which may be larger than the actual quotient,

we check that any additional words of q are zero. We do this by normalising q.

The second test we do is for an inexact division. We check that the the

overflow words turn out to be the same as the overflow from mullow(q*d).

Note that if we wish to make Hensel division efficient for doing an exact

division, say of a 2m - 1 by m division, we merely pass m words of the

dividend in instead of all 2m - 1 words, so that the quotient is also m words.

Once again we don't allow an overflow-in to Hensel division. This wouldn't

give us any kind of chaining property anyway.

Instead, we'd have to do a mulhigh(q, d) and subtract that from the high

part of the chain before continuing, and the mulhigh will accept our overflow

words from the low Hensel division.

In fact, we add a chaining test that does precisely this. We do a Hensel

division on the low n words of a chain, subtract a mulhigh from the high

m words of the chain, then compute the high m words of the quotient using

a second Hensel division.

To make our test code work, we add an ODD flag to randoms_of_len so that only

odd divisors are used with Hensel division.

It isn't clear if it makes sense to allow a carry-in at the top of div_hensel

or not. On the one hand, it might seem logical to allow a carry-in on account

of the way mul_classical works. On the other hand, divrem_hensel1 took a

carry-in, but at the bottom. This was for chaining rather than a read-in of

an extra word. We choose the latter convention, as it seems to make more

sense here and stops the code from becoming horribly complex.

The code for today's post is here: v0.17

Previous article: v0.16 - mullow_classical

Next article: v0.18 - printx_word, nn_printx

BSDNT - v0.16 mullow_classical, mulhigh_classical

2010-10-14T12:17:00.000-07:00

The most logical routine to implement next would be Hensel division. Its
main application is in doing exact division.

For that reason, we might want to focus on Hensel division without
remainder first.

This would take an m word integer and divide it by an n word integer,
division proceeding from right to left. Essentially it gives us division
modulo B^m.

However, before we implement this, it makes sense to think about what its
"inverse" operation might be.

If q is the Hensel quotient of a by d mod B^m, then the inverse operation
can be thought of as multiplication of q by d mod B^m.

Hensel division gives an m word quotient q, so this would imply that its
inverse should be a multiplication of {d, n} by {q, m} mod B^m. We might call
this inverse operation mullow, as it returns the low m words of the product
d*q.

However, we need to be careful with this kind of multiplication. We'd also
like to have a mulhigh which returns the high part of the multiplication,
and we'd like the sum of mullow and mulhigh to be the same as a full mul.

However, there is a problem if mullow merely returns the product mod B^m. Any
carries out of mullow will have been lost. Also, all the word by word

multiplications that contribute to the high word of the product mod B^m
will be thrown away.

To rectify the problem we accumulate an "overflow" out of the mullow
corresponding to the sum of all these high words and carries. As this
overflow is essentially the sum of arbitrary words it may take up two words.

Thus, instead of mullow yielding m words it will yield m + 2 words. We'd
like to pass the extra two words as an "overflow-in" to mulhigh, thus the
logical thing is to return these two words from mullow separately from the
product mod B^m itself.

Hensel division will also return two overflow words. After all, what it
essentially does to the dividend is subtract a mullow of the quotient by the
divisor. So, the overflow from Hensel division will be defined as precisely
the overflow from mullow(q*d).

We manage the overflow by accumulating it in an dword_t. However, as we don't
wish the user to have to deal with dword_t's (these are used in our internal
implementations only), we split this dword_t into two separate words at the
end and return them as an array of two words representing the "overflow".

Today we shall only implement mullow and mulhigh. The first of these is a lot
like a full multiplication except that the addmul1's become shorter as the
algorithm proceeds and the carry-out's have to be accumulated in two words, as
explained.

At the same time we implement mulhigh. This takes two "overflow-in" words and
computes the rest of the product, again in a fashion similar to a full
multiplication.

Our test code simply stitches a mullow and mulhigh together to see that the
chain is the same as a full multiplication.

we have to be careful in that if one does an n by n mullow, the mulhigh that
we wish to chain with this must start with an n-1 by 1 multiplication,
not an n by 1, otherwise the sum of the mullow and mulhigh would contain the
cross product of certain terms twice.

We also have to be careful in the case where the full multiplication is only
a single word by a single word. Here the overflow out of the mullow part is
only a single word and there is no mulhigh to speak of. It merely passes
the overflow from the mullow straight through.

Both mullow and mulhigh must accept all nonzero lengths, as per full
multiplication. This causes a few cases to deal with in mulhigh. This doesn't
seem particularly efficient or elegant, but there seems to be little we
can do about that.

An interesting question is what the inverse of mulhigh is. Essentially,
this is our euclidean divapprox.

There's something slightly unsatisfying here though. Recall that the divapprox
algorithm proceeds by subtracting values q1*d' where q1 is the current quotient
word and d' is what is currently left of the divisor. We throw away one word of
the divisor at each iteration until finally we are left with just two words.

It would be wholly more satisfying if we didn't require this extra word of the
divisor throughout. We'd then be working right down to a single word in the
divisor so that we will really have subtracted a mulhigh by the time the
algorithm completes.

Any modification I can think of making to the euclidean division to make this
seem more natural also makes it much less efficient.

Perhaps some further thought will lead to a more satisfying way of thinking about
these things, which isn't also less efficient in practice.

The code for today is here: v0.16

Previous article: v0.15 - divapprox_classical

Next article: v0.17 - div_hensel

BSDNT - v0.15 divapprox_classical

2010-10-13T15:00:00.000-07:00

During the past few weeks, Brian Gladman has been doing

some tremendous updates, including some random number generators and

making bsdnt work on Windows (MSVC).

We discuss all matters related to bsdnt on our development list:

http://groups.google.co.uk/group/bsdnt-devel?hl=en

Anyhow, now for today's update.

We'd now like to implement a variant of our divrem_classical algorithm. This
time we'd like to just return a quotient, with no remainder. The question is,
can this be done in less time than a full divrem?

At first sight, the answer seems to be no. As we saw in the post on the divrem
algorithm, every word of both the dividend and divisor counts, and we need
to keep the dividend completely updated as the algorithm proceeds, otherwise
we will get the wrong quotient.

So, with the exception of perhaps the final update (which is only needed to
determine the remainder), there doesn't seem to be much we can save.

But what if we allowed the quotient to be approximate, say within 1 of the actual

quotient? In fact, let's demand that if q' is our approximate quotient, that |a -

d*q'| < d. In other words, we allow the quotient to be too large by 1, but not

too small by 1.

Ideally, what we would like is to be doing about half the updating work.
Specifically, we'd like to be truncating both the dividend and divisor as we
go.

Ordinarily we start with something like

AAAAAAAAAAAAAAA /
DDDDDDDDD

To get the first quotient word we shift the divisor and pad with zeroes, thus

AAAAAAAAAAAAAAA /
DDDDDDDD000000

After one iteration, a word of the dividend has been removed, and we then
shift the divisor again.

AAAAAAAAAAAAAA /
DDDDDDDD00000

We continue until we have

AAAAAAAAA /
DDDDDDDD

Now there is only one quotient word left. But we notice that we don't use
most of the remaining dividend or the divisor to determine this quotient
word. In fact, we could almost truncate to

AA /
D

What if we had truncated at this same point all along?

In fact, if we truncate so that the final division is a two word by one
word division (here we have to be careful, in that we are talking about the
number of words *after* normalisation), then clearly our quotient could be
out by as much as two on that final division, by what we have said in an
earlier post. That is of course ignoring any accumulated error along the way.

As we don't wish to multiply the entire thing out to see what we have to do
to correct it, it is clear that this amount of truncation is too great.

So let's truncate one further word to the right in both the dividend and
divisor, so that the final division (to get the final quotient word) is a
three word by two word division.

In fact, in the example above, as there will be five quotient words, there
will be five iterations of the algorithm, after which we want two words of
the divisor remaining. So, we will start with

AAAAAAA.A /
DDDDDD.D

(The decimal points I have inserted are arbitrary, and only for notational
purposes in what follows.)

After five iterations, throwing away one more word of the divisor each time,
we'll be left with

AA.A /
D.D

The first thing to notice is that our previous divrem algorithms, with the
adjustments they made as they went, gave the precise quotient given the
data they started with.

The second thing to notice is that truncating both the dividend and the
divisor at the same point, as above, will not yield a quotient that is too
small. In fact, the quotient we end up with will be the same as what we would
have obtained if we had not truncated the dividend at all, and only truncated
the divisor. Additional places in the dividend can't affect the algorithm.

Truncating the divisor, on the other hand, may result in a different quotient
than we would have obtained without truncation. In fact, as we end up
subtracting less at each update than we would if all those words were still
there, we may end up with a quotient which is too large. The divisor may also
divide more times, because of the truncation, than it would have if it had not
been truncated.

However, it is not enough to merely consider how the quotient changes with
truncation in order to see how far we can be out. We'll likely end up with a
very pessimistic estimate if we do this, because we may suppose that the
quotient can be one too large at each iteration, which is not true.

Instead, the quantity to keep track of is the original dividend minus the
product of the full divisor d and the computed quotient q. At the end of the
algorithm, this is the actual remainder we'll end up with, and we'd like to
keep track of how much our *computed* remainder (what's left of the dividend
after the algorithm completes) differs from this actual remainder.

Essentially, we accumulate an error in the computed remainder due to the
truncation.

Clearly, at each iteration, the word after the decimal point in what remains
of the dividend may be completely incorrect. And we may miss a borrow out of
this place into the place before the decimal point. So after n iterations of
the algorithm, the dividend may become too large by n. Of course n.0 is much
smaller than our original (normalised) divisor d (also considered as a decimal
D.DD...).

At the final step of the algorithm, we will have a dividend which is too large
by at most this amount, and we'll be using a divisor which is truncated to
just two words. However, the latter affects the computed remainder by an
amount much less than the original d (if Q is the final quotient word, it is
as though we added q*0.0DDDDDD to our divisor, so that the full divisor would
go Q times where it otherwise would only go Q-1 times).

So these two sources of error only increase the computed value q'*d + r (where
q' is the computed quotient and r is the computed remainder) by an amount less
than d. Thus, the computed quotient q' can be at most one larger than the
actual quotient q.

This is equivalent to the required |a - d*q'| < d.

So it seems that if we truncate out dividend at the start of the algorithm,
and our divisor after each iteration, we can get an approximate quotient q'
within the required bounds.

We'll leave it to another time to describe the usefulness of an algorithm
which computes a quotient which may be out by one. What we will note is that
we've done computations with much smaller integers. It therefore costs us
significantly less time than a full divrem.

In this week's branch we implement this algorithm. In the test code, we check
the required quotient is the same as the one returned by divrem, or at most
one too large.

The trickiest part is ensuring we truncate at the right point. We want to
finish on the last iteration with two words *after* normalisation of the
divisor.

Actually, if we are really smart, we realise that if d does not need
to be shifted by much to normalise it, we can get away with finishing with
just two *unnormalised* words in the divisor. The error will still be much
less than d.

To be safe, if the number of quotient words is to be n, I check if the
leading word of the unnormalised divisor is more than 2*n. If not, too much
normalisation may be required, and I set up the algorithm to finish with
three unnormalised words instead of two. Otherwise it is safe to finish
with two words in the unnormalised divisor.

The algorithm begins by computing the number of words s the divisor needs to
be to start. This is two more than the number of iterations required to
get all but the final quotient word, since we should have two words in the
divisor at this point. If normalisation is going to be a problem, we add one
to this so that we compute the final quotient word with three unnormalised
words in the divisor.

Now the number of words required to start depends only on the size of the
quotient, and thus it may be more than the number of words d actually has.
Thus we begin with the ordinary divrem algorithm until the number of words
required is less than the number of words d actually has.

Now we truncate d to the required number of words and the dividend to one
more than that. The remaining part of the algorithm proceeds in the same
manner, throwing away a divisor word every iteration.

Only one thing can go wrong, and that is the following: because we are

truncating the divisor at each point, we may end up subtracting too little from

the dividend. In fact, what can happen is that the top word of the dividend

does not become zero after we subtract q*d' (where d' is the truncated divisor).

When this happens, the top word of the dividend may be 1 after the subtraction.

We know that the least significant word of the dividend could be completely

wrong, and the next word may be too large by about as many iterations as

we've completed so far.

Thus, in order to fix our overflow, we subtract from the dividend as much as we

need to in order for the overflow to disappear. We don't mind our dividend

being too large, as we adjust for that in the algorithm. But we cannot allow it to

become too small. Thus we must only subtract from the dividend precisely the

amount required to make the overflow vanish.

We can safely subtract away whatever is in the bottom two words of the

dividend, as this is not even enough to remove the overflow. And then we can

subtract 1 from the whole dividend. This must remove the overflow and is

clearly the least we can get away with subtracting to do so.

The code for today's post is here: v0.15

Previous article: v0.14 - divrem_classical

Next article: v0.16 - mullow_classical

BSDNT - v0.14 divrem_classical

2010-10-06T03:49:00.000-07:00

It's time to implement our schoolboy division routine. I prefer the name

divrem_classical, in line with the multiplication routines.

This function will take a dividend a, divisor d, carry-in ci and returns

a quotient q and remainder r. We'll also need to pass in a precomputed

inverse to speed up the dword_t by word_t divisions that need to occur.

We are going to implement the first algorithm spoken of in the previous

post, namely the one which uses a single "normalised" word of the divisor.

We need to be careful not to just shift the top word of the divisor so

that it is normalised, but if it has more than one word, shift any high

bits of the second word as well. We want the top WORD_BITS bits of the

divisor, starting with the first nonzero bit.

We can use our macro for doing a dword_t by word_t division to get each

new word of the quotient. We start with the carry-in and the most

significant word of a. The macro will automatically shift these by the

same amount as we shifted the leading words of the divisor.

As per the divrem1 function, we require that the carry-in be such that

the quotient won't overflow. In other words, we assume that if the

divisor d is m words, then the top m words of the dividend including the

carry-in, are reduced mod d.

First we need to check that we are not in the special case where truncating

the dividend and divisor to two and one words respectively would cause

an overflow of the quotient word to be computed. This only happens

when the top word of the dividend equals the top word of the divisor, as

explained in the previous post.

If the truncation would cause an overflow in the quotient, we collect a

quotient word of ~0, as discussed in the previous post. If not, we compute

the quotient using our macro.

After this point, the remainder is computed. We allow the function to

destroy the input a for this purpose. We leave it up to the caller to make

a copy of a and pass it to the function, if this is not desired

behaviour.

We do our adjustment so that the quotient is correct. We need a while loop

for this, as mentioned in the previous article.

Finally we write out the quotient word and read in the next word of the

dividend that remains.

In the test code we use our mul_classical and muladd_classical functions to

check that divrem_classical is indeed the inverse of these functions, with

zero remainder and nonzero remainder respectively.

The code for today's post is here: v0.14

Previous article: divrem discussion

Next article: v0.15 - divapprox_classical

BSDNT - divrem discussion

2010-10-02T09:45:00.000-07:00

Today's post is a discussion only, without accompanying code. This is because

the topic is divrem using the "school boy" algorithm, and this is not as

straightforward as one might imagine. The discussion below is informal, and

may contain errors. Please let me know if so and I can make some adjustments

before releasing the code for the next article, where we implement this.

As a consolation, I have released a failed attempt at using C89 macros to make

even more simplifications to our test code. In particular, I tried to make the

chain tests more "automatic". Unfortunately, this didn't work out. It ended up

making the test code too complicated to use and it was going to be way too much

work to convert all of it over to the new test system.

This test code "simplification" was completely abandoned and will not be merged.

But if you are curious you can see it here: gentest

In particular, take a look at generic.c/h and the test code in t-nn.c.

Anyhow, back to today's discussion:

Firstly, what exactly do we mean by the school boy division algorithm?

Certainly I leaned how to do 39 divide 7 at school, and even 339 divide 7.

But how about 10615 divide 1769?

The first step is to divide 1769 into 1061. It doesn't go, so we grab another

digit. So it becomes 1769 into 10615. I have no idea how many times it goes!

Either I missed something at school when I studied "long division", or I only

think I learned an algorithm suitable for problems like this.

In fact, the only way I know to do this division is trial and error, i.e.

to go through all the likely quotient candidates from one to nine.

Using a calculator I find it goes 6 times.

I think I intuitively know how to proceed from here. We place the digit 6

into our quotient, multiply 1769 by 6 and subtract to form a remainder.

But now we have a problem. What if we are dividing 17699949 into 106150000.

Does it still go 6 times? In fact it does not. But how would I know this

without being able to divide 17699949 into 106150000 in the first place!

So I have two problems. Firstly, if I simply truncate the numerator and

denominator to a certain number of digits, even the first digit of my

quotient may be wrong, and secondly, if I use more digits in the first place

my intiuitive "algorithm" doesn't help. It basically says, guess and refine.

That's all very well, but when my digits are 64 bit words, I may not be so

lucky with the guessing thing.

Now an interesting thing to note in the example above is that no matter how

many digits we add to the dividend, 1769 will always go into the first 5

digits 10615 of the dividend, 6 times (with some remainder), no more, no

less.

To see this, we note that 6*1769 is less than 10615 and 7*1765 is greater.

So the only way we could get 1765 to go in 7 times would be to increase

those first five digits of the dividend. Any following digits are irrelevant

to that first digit, 6, of the quotient.

This is a general feature of integer division, and also applies for the

word based (base B) integer division problem.

However, adding digits to the divisor is not the same, as the example above

shows. In fact 17699949 only goes 5 times into 106150000.

Now the question is, how far off can it be if I truncate the dividend and

divisor? How bad can it get?

We can answer this as follows. Note that 17699949 < 17700000. So 10615xxxx

divided by 1769yyyy is greater than or equal to 10615 divided by 1770, which

happens to be 5. So the answer has to be either 5 or 6 times.

What I have done here is simply add 1 to the first four digits of the

divisor, to get a lower bound on the first digit of the quotient. In this

way, I can get a small range of possible values for the first digit of the

quotient. Then I can simply search through the possibilities to find the

correct quotient digit. This matches my intuition of the long division

algorithm at least.

It seems that *perhaps* a quotient digit obtained in this way will be at

most out by 1. In other words, roughly speaking, if X is the first few digits

of the dividend and Y the appropriate number of digits of the divisor (in the

example above, X = 10615 and Y = 1769), then perhaps X / Y >= X / (Y + 1)

>= (X / Y) - 1, where by division here, I mean integer division. Let's call

this condition R.

Let's suppose now that we are working on integers written as binary words

instead of decimal digits. Let's also suppose that X / Y < B. Let's call

this condition S.

Well, it is easy to generate a counterexample to condition R. Let Y = 2,

X = B. Clearly R is not satisfied, even though S is.

But this seems to only be a problem because Y is so small. So, what if we

constrain Y to be at least B/2?

It turns out that condition R still doesn't hold. A counterexample is

Y=B/2, X = B*Y - 2.

However, I claim that X / (Y + 1) >= (X / Y) - 2 does hold, so long as

condition S holds. Let us call this inequality T.

Clearly there does not exist a counterexample to T with X <= Y.

Let us suppose that for a certain Y >= B/2, there is a counterexample to T.

Clearly if X0 is the first such counterexample, then X0 is a multiple of Y,

but not of Y + 1. Nor is X0 - 1 a multiple of Y + 1 (else X0 - 2 would have

been a counterexample).

It is also clear that every multiple of Y from X0 on is a counterexample, as

the left hand side can never "catch up" again.

Thus we have the following result: if there exists a counterexample to T for

some Y >= B/2 and X < BY, then X = (B - 1)*Y is a counterexample.

Substituting this value of X into T yields that T holds if and only if

(B - 1)*Y / (Y + 1) >= B - 3 for all Y >= B/2. Thus, if we can show that this

last inequality holds for all Y >= B/2 then we have proven T.

Note that this inequality is equivalent to (B - 1)*Y >= (Y + 1)*(B - 3), i.e.

2*Y >= B - 3. However, as Y >= B/2, this holds under our hypotheses.

So putting all of the above together, suppose we are dividing A by D and that

Y >= B/2 is the leading word of the divisor D and B*Y > X >= Y is the first

word or two words of the numerator A (whichever is appropriate). Then if Q0

is the leading word of the quotient Q = A / D, then we have shown that

X / Y >= Q0 >= (X / Y) - 2.

In other words, X / Y can be no more than 2 away from the first word of the

quotient Q = A / D that we are after.

This leads us immediately to an algorithm for computing a quotient of two

multiprecision integers.

(i) Let Y be the first WORD_BITS bits of the divisor D, so that

B > Y >= B/2.

(ii) Let X be the first word or two words (appropriately shifted as per Y) of

A such that B*Y > X >= Y.

(iii) Let Q0 = X/Y.

(iv) Set R = A - D*Q0 * B^m (for the appropriate m), to get a "remainder".

While R < 0, subtract 1 from Q0 and add D * B^m to R.

(v) Write down Q0 as the first word of the quotient A / D and continue on by

replacing A by R and returning to step (i) to compute the next word of the

quotient, etc.

Another algorithm can be derived as follows. Instead of using Y >= B/2 in the

above algorithm, let's choose Y >= B^2/2. A similar argument to the above

shows that X / (Y + 1) >= (X / Y) - 1 for Y >= B^2/2 and B*Y > X >= Y. It boils

down to showing that (B - 1)*Y / (Y + 1) >= B - 2 for Y >= B^2/2, i.e. that

Y >= B - 2, which is clearly satisfied under our hypotheses.

The algorithm is precisely the same, but at step (iv) we can replace the while

statement with an if statement and perform at most one adjustment to our

quotient Q0 and to R.

Now we return to step (v), in which we said that we could just continue

on from step (i) to compute the next word of the quotient Q. If we do this

and set A = R then compute the new X, what we would like is something like

the divrem1 algorithm we implemented, where, (after possibly some kind of

initial iteration that we handled specially), it is always true that the new

X is two words and has its first word reduced mod Y.

However, this does not follow from 0 <= R < D*B^m, and it may be that the

first word of the remainder is equal to Y! This is due to the truncation of

R to get Y.

It is clear from 0 <= R < D*B^m that if X is the first two words of the

remainder that X <= B*Y. So to make the algorithm continue in the way we'd

like, we only have to deal with the special case where X = B*Y.

We know that in the first algorithm above, the next word of the quotient may

be B - 1 or B - 2, since we know already that it is not B. We must multiply

out and check which it is.

In the second algorithm above, where Y >= B^2/2, the only possibility for

the next quotient word is B - 1, as we know it is not B.

In the next article we will look at implementing the first of the two

algorithms above, leveraging the code we already have for computing and using

a precomputed inverse.

Previous article: v0.13 - muladd1, muladd

Next article: v0.14 - divrem_classical

BSDNT - v0.13 muladd1, muladd

2010-10-01T05:13:00.000-07:00

We noted in the last article that multiplication need not take a carry-in

and that it doesn't seem related to the linear functions we have been

implementing.

Instead, we think of multiplication in a different way, as an inverse of

division. We'll soon implement division with remainder, i.e. given a and

d, find q and r such that a = d*q + r, where 0 <= r < d

If d and q are of lengths m and n respectively, then r is of length at most

m and a is either of length m + n or m + n - 1.

Multiplication corresponds to the case where r = 0.

As we will certainly not wish to have to zero pad our division out to

m + n limbs if a is in fact only m + n - 1 limbs, it makes sense that our

division will take a carry-in. For this reason, it made sense for our

multiplication to yield a carry-out, i.e. it will write m + n - 1 limbs

and return an m + n - th limb, which may or may not be 0.

When r is not zero in our division, the inverse operation would take a

value r of n limbs and add d*q to it, returning the result as a value a of

m + n - 1 limbs (and a carry which may or may not be zero). We call this

routine muladd. In the case where r and a happen to be aliased, the result

will be written out, overwriting a.

We first implement nn_muladd1 which is identical to addmul1 except that

muladd1 writes its result out to a location not necessarily coincident with

any of the inputs. In other words, nn_addmul1 is an in-place operation

whereas nn_muladd1 is not.

Next we implement nn_muladd_classical. It takes an argument a of length m

to which it adds the product of b of length m and c of length n. The result

may or may not alias with a.

We also implement a version of the linear and quadratic muladds which writes

out the carry, naming the functions with the nn_s_ prefix, as per the

convention initiated in our nn_linear module.

As for multiplication, we don't allow zero lengths. This dramatically

simplifies the code.

An important application of the muladd function will be in chaining full

multiplication. We'll discuss this when we get to dealing with multiplication

routines with better than quadratic complexity.

Our test code simply checks that muladd is the same as a multiplication and

and an addition.

The github repo for this post is here: v0.13

Previous article: configure and assembly

Next article: divrem discussion

BSDNT - configure and assembly improvements

2010-09-25T12:59:00.000-07:00

It's time we added architecture dependent assembly support to bsdnt.

Here is how we are going to do it. For each of the implementation files

we have (nn_linear.c, nn_quadratic.c), we are going to add a _arch file,

e.g. nn_linear_arch.h, nn_quadratic_arch.h, etc.

This file will be #included in the relevant implementation file, e.g.

nn_linear_arch.h will be #included in nn_linear.c.

These _arch.h files will be generated by a configure script, based on

the CPU architecture and operating system kernel. They will merely

include a list of architecture specific .h files in an arch directory.

For example we might have nn_linear_x86_64_core2.h in the arch directory,

which provides routines specific to core2 processors running in 64 bit

mode.

In these architecture specific files, we'll have inline assembly routines

designed to replace various routines in the pure C implementation files

that we have already written. They'll do this by defining flags, e.g.

HAVE_ARCH_nn_mul1_c, which will specify that an architecture specific

version of nn_mul1_c is available. We'll then wrap the implementation of

nn_mul1_c in nn_linear.c with a test for this flag. If the flag is defined,

the C version will not be compiled.

In order to make this work, the configure script has to work out whether

the machine is 32 or 64 bit and what the CPU type is. It will then link

in the correct architecture specific files.

At the present moment, we are only interested in x86 machines running on

*nix (or Windows, but the architecture will be determined in a different

way on Windows).

A standard way of determining whether the kernel is 64 bit or not is to

search for the string x86_64 in the output of uname -m. If something else

pops out then it is probably a 32 bit machine.

Once we know whether we have a 32 or 64 bit machine, we can determine the

exact processor by using the cpuid instruction. This is an assembly

instruction supported by x86 cpus which tells you the manufacturer, family

and model of the CPU.

We include a small C program cpuid.c with some inline assembly to call the

cpuid instruction. As this program will only ever be run on *nix, we can

make use of gcc's inline assembly feature.

When the parameter to the cpuid instruction is 0 we get the vendor ID,

which is a 12 character string. We are only interested in "AuthenticAMD"

and "GenuineIntel" at this point.

When we pass the parameter 1 to the cpuid instruction, we get the processor

model and family.

For Intel processors, table 2-3 in the following document gives information

about what the processor is:

http://www.intel.com/Assets/PDF/appnote/241618.pdf

However the information is out of date. Simply googling for Intel Family 6

Model XX reveals other models that are not listed in the Intel documentation.

The information for AMD processors is a little harder to come by. However,

one can essentially extract the information from the revision guides, though

it isn't spelled out clearly:

http://developer.amd.com/documentation/guides/Pages/default.aspx#revision_Guides

It seems AMD only list recent processors here, and they are all 64 bit.

Information on 32 bit processors can be found here:

http://www.sandpile.org/ia32/cpuid.htm

At this point we'd like to identify numerous different architectures. We

aren't interested in 32 bit architectures, such as the now ancient Pentium

4 or AMD's K7. Instead, we are interested when 32 bit operating system

kernels are running on 64 bit machines. Thus all 32 bit CPUs simply identify

as x86.

There are numerous 64 bit processors we are interested in:

64 bit Pentium 4 CPUs were released until August 2008. We identify them as p4.

All the 64 bit ones support SSE2 and SSE3 and are based on the netburst

technology.

The Intel Core Solo and Core Duo processors were 32 bit and do not interest us.

They were an enhanced version of the p6 architecture. They get branded as x86

for which only generic 32 bit assembly code will be available, if any.

Core 2's are very common. They will identify as core2. They all support SSE2,

SSE3 and SSSE3 (the Penryn and following 45nm processors support SSE4.1 - we

don't distinguish these at this stage).

Atoms are a low voltage processor from Intel which support SSE2, SSE3 and are

mostly 64 bit. We identify them as atom.

More recently Intel has released Core i3, i5, i7 processors, based on the

Nehalem architecture. We identify these as nehalem. They support SSE2, SSE3,

SSSE3, SSE4.1 and SSE4.2.

AMD K8's are still available today. They support SSE2 and SSE3. We identify

them as k8.

AMD K10's are a more recent core from AMD. They support SSE2, SSE3 and SSE4a.

We identify these as k10. There are three streams of AMD K10 processors,

Phenom, Phenom-II and Athlon-II. We don't distinguish these at this point.

So in summary, our configure script first identifies whether we have a 32 or

64 bit *nix kernel. Then the CPU is identified as either x86, p4, core2,

nehalem, k8 or k10, where x86 simply means that it is some kind of 32 bit CPU.

Our configure script then links in architecture specific files as appropriate.

The only assembly code we include so far are the nn_add_mc functions we wrote

for 64 bit core2 and k10. As these are better than nothing on other 64 bit

processors from the same manufacturers, we include this code in the k8

specific files until we write versions for each processor. We also add an

nn_sub_mc assembly file for Intel and AMD 64 bit processors.

The configure script includes the architecture specific .h files starting from

the most recent processors so that code for earlier processors is not used

when something more recent is available.

The github branch for this revision is here: asm2

Previous article: v0.12 mul_classical

Next article: v0.13 - muladd1, muladd

BSDNT - v0.12 mul_classical

2010-09-24T12:25:00.000-07:00

I've been on holidays for a couple of days, hence the break in transmission.

(For academics: the definition of a holiday is a day where you actually

stop working and do nothing work related for that whole day. I realise the

definition is not widely known, nor is it well-defined.)

Note, due to accidentally including some files in today's release that I

intended to hold over till next time, you have to do ./configure before

make check. This detects your CPU and links in any assembly code

that is relevant for it. More on this when I do the actual blog post about it.

It's time to start implementing quadratic algorithms (i.e. those that

take time O(n^2) to run, such as classical multiplication and division).

Before we do this, we are going to reorganise slightly. The file which

is currently called nn.c we will rename nn_linear.c to indicate that it

contains our linear functions. We'll also rename t-nn.c to t-nn_linear.c.

The macros in that file will be moved into test.h.

We also modify the makefile to build all .c files in the current directory

and to run all tests when we do make check. A new directory "test" will

hold our test files.

Finally, we add an nn_quadratic.c and test/t-nn_quadratic.c to hold the

new quadratic routines and test code that we are about to write.

The first routine we want is nn_mul_classical. This is simply a mul1

followed by addmul1's. Of course this will be horrendously slow, but once

again we defer speeding it up until we start adding assembly routines,

which will start happening shortly.

In a departure from our linear functions, we do not allow zero lengths.

This dramatically simplifies the code and means we do not have to check

for special cases.

There does not appear to be any reason to allow a carry-in to our

multiplication routine. An argument can be made for it on the basis of

consistency with mul1. However, the main use for carry-in's and carry-out's

thus far has been for chaining.

For example, we chained mul1's s*c and t*c for s and t of length m and

c a single word in order to compute (s*B^m + t)*c.

But the analogue in the case of full multiplication would seem to be

chaining to compute (s*B^m + t)*c where s, t and c all have length m. But

it probably doesn't make sense to chain full multiplications in this way

as it would involve splitting the full product t*c say, into two separate

parts of length m, which amongst other things, would be inefficient.

It might actually make more sense to pass a whole array of carries to our

multiplication routine, one for every word of the multiplier. However it is

not clear what use this would be. So, for now at least, we pass no carry-ins

to our multiplication routine.

We settle on allowing a single word of carry out. This may be zero if the

leading words of the multiplicands are small enough. Rather than write this

extra word out, we simply return it as a carry so that the caller can decide

whether to write it.

The classical multiplication routine is quite straightforward, but that plus the

rearrangement we've done is more than enough for one day.

The github branch for this release is here: v0.12

Previous article: v0.11 - generic test code

Next article: configure and assembly

BSDNT - v0.11 generic test code

2010-09-20T18:37:00.000-07:00

It's time we improved our test framework for all the functions we are
adding, before things get out-of-control.

The best strategy is to actually write an interpreter for a small
expression parser which allows us to test arbitrary expressions.
However, that would be a distraction at this point, so we opt to
simply add some convenience functions.

In particular, we'll add some new random functions which generate
lists of random words and random garbage collected nn's. Our only aim
here is to simplify our test code and cut down the number of lines of
it.

We don't want to make it *more* complicated for people to use, and
C89 places some hard limits on us in that regard.

The first step is to add two new files, test.c and test.h which
will contain our garbage collection and random functions.

We begin by adding an enum type_t, which parameterises all the different
types our garbage collection can clean up. At the moment we only need
one type, NN.

Next we add a variadic function (variable number of arguments) which
creates a whole bunch of random words. Unfortunately C89 doesn't support
variadic macros, so we need to pass pointers to the words, so that the
actual words can be modified.

This random function takes a flag which specifies what type of random
number we want. Specifically we'll have flags for ODD, NONZERO, ANY.

C is also so stupid that we actually need to tell it explicitly, or mark
somehow, the number of arguments we are giving to the variadic function.
In C99 one can work around this by wrapping a variadic function with a
variadic macro. But we don't have that facility here. We choose to simply
pass NULL as the final argument to variadic functions.

We also introduce a global garbage stack which will keep a pointer to all
the objects allocated so far, so that a garbage collector can clean them
up later on, when called. The fact that this garbage is global makes it not
threadsafe, but we are only using the test support for test code at this
point and so we don't care about thread safety. Later we could pass a
context around as a parameter if we wanted to make it threadsafe.

We also add a gc _cleanup function, which cleans away all the garbage, so
that memory leaks don't occur. Later we could expand this function to do
redzone checking and various other automatic tests on the garbage it is
disposing of.

Now we can add some convenience functions for generating random
multiprecision integers.

The randoms_of_len function again takes a NULL terminated list of
*pointers* to nn_t's and both initialises them to the given length and
sets them to random limbs. We can specify various flags here, such as
ANY, FULL, the latter returning a multiprecision integer with exactly
the given number of limbs, with the top limb being nonzero (unless
the number of limbs requested is zero).

Now that we have all this, in our test file test.h, we define a macro
TEST_START which takes a number of iterations and simply sets up a loop
with that many iterations. A TEST_END function then calls garbage collection
to clean up at the end.

Finally, we add a test_generics function, which generates some random values
using our new generics and then cleans up.

By running a program called valgrind over our test code, we can see if all
the memory used by our generics actually got cleaned up after use.

We now roll out these changes to the rest of our test code. We note that the

entire test file is now about two-thirds of the length it was before we rewrote

it!

The new branch is on github: v0.11

Previous article: v0.10 - mod1_preinv

Next article: v0.12 - mul_classical

BSDNT - v0.10 mod1_preinv

2010-09-19T07:49:00.000-07:00

We have one more linear function to implement before we move on to some

assembly language optimisation and then quadratic functions, namely mod1.

It returns the remainder after division by a single word.

One might think that this cannot be computed any faster than divrem1,

however the following algorithm allows us to perform the computation

slightly faster.

The algorithm is credited to Peter Montgomery. It is based on the following

observation. Suppose we have computed, in advance, the values

b2 = B^2 mod d and b3 = B^3 mod d.

Then a number of the form a0 + a1 * B + a2 * B^2 + a3 * B^3 can be

reduced mod d by computing s = a0 + a1 * B + a2 * b2 + a3 * b3, then

later reducing s mod d.

Note that as d can be at most B - 1, the values bi are at most B - 2.

The values ai are at most B - 1. Thus ai * bi is at most B^2 - 3*B + 2.

This means that when summing up s, we are adding at most B^2 - 3*B + 2

to B^2 - 1 at each step. The total is at most 2*B^2 - 3*B + 1. In

other words, there might be a carry of *1* into a third limb. But we

already know that b2 = B^2 mod d. Thus we can get rid of this 1 from

the third limb by adding b2 to our total in its place. The total will

then be at most B^2 - 2*B - 1, which fits easily into two limbs.

We may need to do this adjustment twice when summing up. It is an

unpredicted branch to do this correction, but with a deep enough pipeline

it is possible the processor will compute both paths, minimising the cost

of a misprediction.

Thus, with some precomputation, we can reduce four limbs of our dividend

to two limbs with 2 multiplications. Another advantage is that these

multiplications are independent. This gives the processor lots of

opportunity to pipeline them and compute them quickly.

One tricky implementation detail is in testing when s overflows two limbs.

We can accumulate into three limbs, but this is not quite efficient.

One way of testing for overflow of a double word addition is to check if

the result is smaller than one of the summands. If so, an overflow

occurred and we can make our adjustment.

In the case where d is at most B/3, we can make another optimisation. Also

compute b1 = B mod d and perform a third multiplication a1 * b1. We have

that bi is less than B/3. Thus the three products ai*bi are at most (B/3 - 1)*(B - 1)

= B^2/3 - 4*B/3 + 1. Clearly s is now at most B^2 - 3*B + 2 and no overflow

occurs. Therefore no tests or adjustments are required to compute s.

I will leave this optimised case as an exercise and just implement the

generic case for now.

At the end of the algorithm we must reduce a double word mod d. This could

be optimised with a precomputed inverse a la divrem1, but again I skip this

optimisation and leave it as an exercise.

At the start of the algorithm we have to set up differently depending on

whether there are an even or odd number of words (including the carry-in).

If there is an odd number in total, we need to do a single reduction of

the extra word so that an even number remain.

We add a new mod_preinv1_t type and a function precompute_mod_inverse1

to compute it. This computes B, B^2 and B^3 mod d. For now this precomputation

is not optimised to use a precomputed inverse. This would actually save

significant time, but again I leave it as an exercise to optimise this.

The end result looks a little messy compared to the algorithms implemented

so far. I wonder if anyone can find any simplifications of the code.

Here v0.10 is the github branch for this post.

Previous article: v0.9 - divrem1_hensel

Next article: v0.11 - generic test code

BSDNT - v0.9 divrem1_hensel

2010-09-18T05:11:00.000-07:00

The next division function we'll introduce is so-called Hensel division, or

right-to-left division. This starts at the least significant word and applies

the usual division algorithm back towards the most significant word.

Of course, now any "remainder" will be at the most significant word end and

have a completely different meaning.

Hensel division is actually just division mod B^m. Thus if q is the Hensel

quotient, you have a = q * d mod B^m. In fact, if r is the Hensel remainder,

you have a = q*d + r*B^m.

If a fits in m words and the division is exact (i.e. a = q*d), the Hensel

division and ordinary (euclidean) division return the same quotient. If the

division is not exact, this is no longer the case.

However, it may be that chaining Hensel divisions together does produce an

exact division, even when division over the bottom part of the chain

wouldn't produce an exact division. Thus, the ability to chain Hensel

divisions, and thus the ability to return the Hensel remainder, is

important.

Well, now we come to the first issue if we want to implement this. What is

the C operator for Hensel division? Actually, there doesn't seem to be one.

We have to implement it ourselves.

If a0 is a word from our dividend, and d is our single word divisor, then

we require q0 such that d * q0 = a0 mod B. In other words, we want q0 and

r0 such that d * q0 = a0 + r0*B.

Suppose we can find a value dinv such that dinv * d = 1 mod B. Then

we have that dinv * d * q0 = dinv * a0 mod B, i.e. q0 = dinv * a0 mod B,

which is a single word-by-word mullow.

Once we have q0 we can compute d * q0, the high word of which is the Hensel

"remainder" r0. We need to subtract this from the next limb of our divisor

before continuing.

We'll make the convention that carry-ins and carry-outs from our Hensel

division are positive rather than negative. We'll just remember to subtract

them. This way, if q is our Hensel quotient and r our Hensel remainder,

then d * q = a + r * B^n.

So how do we compute dinv? Firstly, we better require that d is odd, else

dinv * d = 1 mod B is an impossibility. With this restriction, we have

gcd(d, B) = 1 and therefore the extended euclidean algorithm lets us find

s, t such that s*d + t*B = 1. Then modulo B we have s*d = 1 mod B, i.e.

dinv = s is the precomputed inverse that we are after.

Another way to solve dinv * d = 1 mod B is to use Hensel lifting. This

works by first solving the equation mod 2, then use that solution to

solve it mod 4 and so on. This can all be done with nothing more than

multiplications.

Firstly, we know a solution mod 2, namely 1 (as d is odd), i.e. if v1 = 1

then we have v1 * d = 1 mod 2. Now suppose we have a solution vk mod 2^k,

i.e. vk * d = 1 mod 2^k. We'd like to find a value wk such that

(vk + 2^k * wk) * d = 1 mod 2^2k. This is always possible -- I omit the

easy proof -- so we only need to solve this equation for wk, i.e.

(d * vk - 1) = -2^k * wk * d. We can get rid of the d on the right hand

side by multiplying by vk, i.e. 2^k * w = ((1 - d * vk) * vk) (mod 2^2k).

In other words, two multiplications will suffice to compute

2^k * wk (mod 2^2k) from vk. Then v{k+1} = vk + 2^k * wk is the inverse

of d mod 2^2k.

As a further optimisation of all this, one can just work mod 2^64 all

the way through the computation. In other words, start with v1 = 1

and compute v{k+1} = vk + ((1 - d * vk) * vk) (mod 2^64), six times.

To get from a solution mod 2 to a solution mod 2^64 clearly only takes

six Hensel lifting steps. Of course, with a lookup table mod 2^8 to start

from, instead of starting from a solution mod 2, one can do it in three

steps, requiring just six word-by-word multiplications (mullow only).

However, today I am feeling very lazy and I leave it as an exercise for

someone to implement the Hensel lifting lookup table method.

We add a simple type to hold our precomputed Hensel inverse, and a function

for computing it.

The actual nn_divrem_hensel_1_preinv function is again a straightforward

loop. We ensure our remainder is positive at each step and ensure that

we propagate any borrows from subtractions at each step.

The test code for the Hensel division is very similar to that for the

ordinary euclidean division except that d must be odd and the chaining

test will proceed right-to-left.

OK, that is enough brain strain for one day!

Here v0.9 is the github repository for today's post.

Previous post: v0.8 - divrem1_preinv

Next post: v0.10 - mod1_preinv

BSDNT - v0.8 divrem1_preinv

2010-09-16T18:18:00.000-07:00

Today we'll implement a not-stupid division function. Since the earliest

days of computers, people realised that where possible, one should

replace division with multiplication. This is done by computing a

precomputed inverse v of the divisor d and multiplying by v instead of

dividing by d.

One particularly ingenious example of this, which was developed recently,

is the algorithm of Möller and Granlund. It can be found in this paper:

http://www.lysator.liu.se/~nisse/archive/draft-division-paper.pdf

If one has a precomputed inverse v as defined at the start of section 3

of that paper, then one can compute a quotient and remainder using

algorithm 4 of the paper using one multiplication and one mullow.

There is one caveat however. The divisor d must be "normalised". This

means that the divisor is a word whose most significant bit is a 1.

This is no problem for us as we can simply shift n to the left so that it is

normalised, and also shift the two input limbs to the left by the same

amount. The quotient will then be the one we want and the remainder

will be shifted to the left by the same amount as n. We can shift it back

to get the real remainder we are after.

We can use this algorithm in our divrem1 function. However, rather than

shift the remainder back each iteration, as it is the high limb of the

next input to the same algorithm we can simply pass it unmodified to the

next iteration of the divrem1 function.

The name of our function is nn_divrem1_preinv. We also introduce a type,

preinv1_t, to nn.h which contains the precomputed inverse (dinv) and the

number of bits n has been shifted by (norm). We add a function,

precompute_inverse1, which computes this.

In order to compute the number of leading zeroes of n, we simply use a gcc

intrinsic __builtin_clzl which does precisely what we are after. This gives

the number of bits we have to shift by.

We implement the Möller-Granlund algorithm as a macro to save as many cycles

as we can. This time we really have to be careful to escape the variables we

introduce inside the macro, to prevent clashes with parameters that might be

passed in. The standard way of doing this is to prepend some underscores to

the variables inside the macro. We call the macro divrem21_preinv1 to signify

that it is doing a division with dividend of 2 words and divisor of 1 word.

Once this is done, the nn_divrem1_preinv algorithm is very simple to implement.

It looks even simpler than our divrem1_simple function. The test code is very

similar as well. The new function runs much faster than the old one, as can

easily be seen by comparing the pauses for the different functions when

running the test code.

Later we'll have to introduce a timing and profiling function so that we can

be more explicit about how much of a speedup we are getting.

I haven't experimented with the divrem21_preinv macro to see if I can knock

any cycles off the computation. For example, careful use of const on one or

more of the variables inside the macro, or perhaps fewer word_t/dword_t

casts might speed things up by a cycle or two. Let me know if you find a

better combination.

Of course this is really just doing optimisations that the C compiler should

already know about. The best way to get real performance out of it is to add

some assembly language optimisations, which we'll eventually do.

The github branch for today's post is here: v0.8

Previous article: bsdnt assembly language

Next article: v0.9 - divrem1_hensel

BSDNT - assembly language

2010-09-14T18:20:00.000-07:00

Today a few of us (Antony Vennard, Brian Gladman, Gonzalo Tornaria and myself) had a go at writing some assembly language.

The results of our labour, for AMD64 can be found on the asm branch on github:

http://github.com/wbhart/bsdnt/tree/asm

This is not a proper revision of bsdnt, but a heavily hacked version, just to demonstrate the principles.

In the nn.c file is some inline assembly code which implements the nn_add_mc function in x64 assembler. The format used is Intel format assembly code, instead of the usual AT&T syntax used by inline gcc normally. We made a change to the makefile so that gcc would use the Intel format assembly code pervasively. Note, this will cause gcc to miscompile bsdnt for many revisions of gcc. This is due to bugs in gcc. Specifically gcc versions 4.4.1 and 4.4.4 work, however.

At this point we haven't unrolled the loop (repeated the contents of the loop over and over to save some time on loop arithmetic). But the results are still pretty nice.

The straight C version we had before took about 12 cycles per word (with gcc 4.4.1 on an AMD Opteron K10 machine). With this new assembly code it takes about 3 cycles per word.

With loop unrolling we might hope to halve that again, but this messy optimisation will have to wait, otherwise work on the functionality of the library will slow right down.

We also wrote some assembly code for Core 2 machines. This runs in about 4.5 cycles per word. At the moment this is not committed anywhere, but you can find a copy here:

http://groups.google.co.uk/group/bsdnt-devel/msg/15e2883d94014196?hl=en

The important thing to realise about assembly language is it is highly architecture dependent. Our AMD64 code is actually slower than C on an Intel Core 2 machine! The assembly code needs to be prepared for the machine in question. Moreover, this code won't work at all on Windows 64 or on any 32 bit machine!

At some point we'll introduce a configuration script which will select the correct assembly code for the architecture being used.

To actually time this code, we simply hacked t-nn.c to only run one of the nn_add_m tests, but with the nn_add_m executed 1000 times and the actual test that the code works, turned off.

We set the number of words being added in each nn_add_m to 24, as the test machine we used was 2.4GHz. This means the number of seconds the test takes to run is approximately the number of cycles per word. This is of course a temporary hack, just to illustrate working assembly code (the test will pass if you remove the loop that calls nn_add_m a thousand times instead of once).

Previous article: v0.7 - divrem1simple

Next article: v0.8 - divrem1_preinv

BSDNT - v0.7 divrem1_simple

2010-09-13T10:24:00.000-07:00

It's time for our first division function. This will again be a linear function

in that it divides a multiprecision integer by a single word and returns the

(single word) remainder.

This is the first time we can do substantially better than a simple for loop

(other than when we coded nn_cmp and other functions, which allowed early exit).

To get us going, we'll write a really dumb division function first. It'll just

use C's "/" and "%" operators in a for loop.

C is supposed to optimise simultaneous occurrences of / and % together in the

same piece of code to use a single processor instruction where available.

However, even with that optimisation, our function will be exceedingly slow.

We'll see why that is when we implement a not-stupid divrem1, next time around.

To emphasise its stupidity, we'll call our exceedingly dumb division function

nn_divrem1_simple. Perhaps we can use it in our test code to compare against

when we write the more sophisticated version.

Hopefully our naming conventions won't be defeated by this function, which

proceeds left to right (or most significant word first).

One doesn't usually need to write the remainder down immediately following

the quotient (though one could imagine a function operating that way).

However, it is convenient to be able to accept a "remainder-in", essentially

thought of as an extra limb of the dividend, so that divrem can be chained.

In this way, divrem1 is most similar to right shift in semantics (the latter

is just division by a power of two, after all).

Now, if we blindly try to implement this function, we quickly discover a

problem. Even if we work with double words, a division of a two limb quantity

by a single limb can be problematic.

We can illustrate this using decimal digits, instead of machine words. Suppose

we divided 94 by 7. The resulting quotient is 13 with remainder 3 (I hope!).

The point is, the quotient takes up two "words" 1 and 3, not just a single

word.

We can get around it by reducing the top word first. We have 9 divide 7 is 1

remainder 2. Thus the first word of our quotient is 1. Now we are left with

24 divide 7. This time the quotient is 3 with remainder 3 and everything is

fine. Note the remainder, 3, is less than 7.

The critical thing to note is that once we get things started, everything

works fine from then on. The top word in our double words are always reduced

modulo the divisor after the first iteration.

So the almost trivial algorithm almost works, except for the lead-in, where

we have to get things right. After that, we can divide our dividend word on

word until we reach the bottom, where we'll have a single limb remainder.

Now as mentioned, it is also convenient to supply a "carry-in". This acts like

an additional limb of the dividend. Assuming this carry-in was reduced mod

the divisor, we find that we end up with as many words in our quotient as we

had in our original dividend.

Clearly however, if the carry-in is not reduced, we will end up with even one

more quotient limb! However, there is an easy way to avoid this problem.

We simply decide to restrict the carry-in to be reduced modulo the divisor!

This restriction is not a problematic restriction because we will still be

able to chain our division functions together with this restriction (so long

as we use the same divisor throughout).

Anyhow, with this restriction, the quotient will have precisely the same number

of limbs as the dividend, and since we'll start with the carry-in, not the top

limb of the dividend, and that carry-in will already be reduced, there is no

bootstrapping iteration required, i.e. the problem of a non-reduced first limb

as explained above, simply doesn't exist.

In this way we avoid any complex feed-in and wind-down code, and as for all

the functions so far, the length m can be zero. The result is satisfyingly

symmetric, even if it is terribly slow.

In the test code, if q = a / c with remainder r, we can check that a = q*c + r.

We can use the functions we have already implemented, including our mul1_c code

to verify this.

We also check that passing in a reduced carry-in results in a correct result,

by chaining two divisions together.

This will do for today, as a more serious implementation of divrem1 will

require much more serious code. We'll take a look at that next time.

The github branch here v0.7 is for this post.

Previous article: bsdnt v0.6 - addmul1, submul1

Next article: bsdnt assembly language

BSDNT - v0.6 addmul1, submul1

2010-09-12T17:04:00.000-07:00

At last we come to addmul1 and submul1. After implementing these, we would

already be able to implement a full classical multiplication routine, using the

naive O(n^2) multiplication algorithm.

It turns out that addmul1 and submul1 are quite similar to mul1. Instead of

just writing out the result, we have to add it to, or subtract it from the

first operand.

These are really combined operations, i.e. combined mul1 and add_m or sub_m

respectively.

The reason we introduce these is that the processor has multiple pipelines

and merging two operations like this gives us the chance to push more

through those pipelines. We'll add more combined operations later on.

The fact that we can add b[i]*c, a[i] and ci and not overflow a double word

needs some justification.

Clearly b[i] and c are at most B-1. Thus b[i]*c is at most B^2-2B+1. And

clearly a[i] and ci are at most B-1. Thus the total is at most

B^2-2B+1 + (B-1) + (B-1) = B^2-1, which just fits into a double word.

The dword_t type allows us to get away with adding all these in C as though

we don't care about any overflow (in fact we know there is none). Of course

that doesn't make it efficient. This function is still a candidate for

assembly optimisation and loop unrolling later on.

The submul function is essentially the same, so long as we have taken care to

perform operations in the right order.

The naming of our functions causes us slight difficulties here. We'd ideally

like versions of addmul1 and submul1 which write the carry out to the high

limb and versions which add/sub it from the high limb. We opt for the latter

for now. After all, if one were accumulating addmuls, this is what one would

require most often.

Maybe someone reading this has a better idea how to handle these.

We add a test which checks that doing two addmuls in a row is the same as

doing a single addmul with the multiply constant equal to the sum of the

original two. We also add a test to check that chaining addmuls works. We

add similar tests for submul.

We delay coding up a quadratic multiplication basecase as there are a few more

linear functions to work on, most notably the various kinds of division and

remainder functions. These are all interesting functions to write.

The github repo for this version is here: v0.6

Previous article: v0.5 add, sub, cmp

Next article: v0.7 divrem1_simple

BSDNT - v0.5 add, sub, cmp

2010-09-11T12:44:00.000-07:00

It's time we added a full add and sub routine. These perform addition and
subtraction with possibly different length operands. For simplicity, we
assume the first operand is at least as long as the second. This saves us
a comparison and a swap or function call.

Add is simply an add_m followed by an add1 to propagate any carries right
through to the end of the longest operand. A similar thing holds for sub.

We'd like to code these as macros, however we'd also like to get the return
value. Thus, we use a static inline function, which the compiler has the
option of inlining if it wants, saving the function call overhead.

We add tests for (a + b) + c = (a + c) + b where b and c are at most the same
length as a. We also check that chaining an add with equal length operands,
followed by one with non-equal operands, works as expected. We add similar
tests for subtraction too.

The next function we add is comparison. It should return a positive value
if its first operand is greater than its second, a negative value if it is
the other way, and zero if the operands are equal.

As for equal_m and equal, we introduce different versions of the function
for equal length operands and possibly distinct operands.

The laziest way to do comparison is to subtract one value from the other
and see what sign the result is. However, this is extremely inefficient
in the case that the two operands are different. It's very likely that
they already differ in the most significant limb, so we start by checking
limb by limb, from the top, until we find a difference.

Various tests are added. We test that equal things are equal, that operands
with different *lengths* compare in the right order, and operands with
different *values* but the same length compare in the right order.

We also take the opportunity to do a copy-and-paste job from our cmp test
code and quickly generate some test code for the nn_equal function, which
didn't have test code up to this point.

From the next section onwards, we will only deal with one or two functions
at a time. All the really simple functions are now done, and we can move on
to more interesting (and useful) things.

The github branch for this post is here: v0.5

Previous article: 0.4 - add1, sub1, neg, not

Next article: v0.6 addmul1, submul1

BSDNT - v0.4 add1, sub1, neg, not

2010-09-10T07:30:00.000-07:00

Today we introduce some more convenience functions, neg and not, and we
also add the functions add1 and sub1.

All of these functions operate slightly differently to the ones we have
introduced so far.

The functions add1 and sub1 simply add a single word to a multiprecision
integer, propagating any carries/borrows all the way along.

The main loops of add1 and sub1 need to stop if the carry becomes zero. This
is for efficiency reasons. In most cases when adding a constant limb to a
multiprecision integer, only the first limb or two are affected. One doesn't
want to loop over the whole input and output if that is the case.

However, we must be careful, as in the case where the input and output are
not aliased (at the same location), we still need to copy the remaining
limbs of the input to the output location.

We add tests that (a + c1) + c2 = (a + c2) + c1 and do the same thing for
subtraction.

We also check that chaining of add1's and chaining of sub1's works. Until
we can generate more interesting random test integers this test doesn't
give our functions much of a workout. We eventually want to be able to
generate "sparse" integers, i.e. integers with only a few binary 1's or a
few binary 0's in their binary representation. The latter case would be
interesting here as it would test the propagation of carries in our add1
and sub1 functions. We'd also eventually like to explicitly test corner
cases such as multiprecision 0, ~0, 1, etc.

A final test of add1/sub1 that we add is a + c1 - c1 = a.

The not function is logical not. It complements each limb of the input. It
is a simple for loop.

The neg function is twos complement negation, i.e. negation modulo B^m. In
fact, twos complement negation is the same as taking the logical not of the
integer, then adding 1 to the whole thing. The implementation is similar to
add1, except that we complement each limb after reading it, but before adding
the carry.

One difference is that we still need to complement the remaining limbs after
the carry becomes zero, regardless of whether the input and output are
aliased.

The carry out from (neg a) is notionally what you would get if you were
computing 0 - a. In other words, the carry is always 1 unless a is 0. In
order to allow chaining, neg must notionally subtract the carry-in from the
total.

We test that (not (not a)) = a and that neg is the same as a combination of
not and add1 with constant 1. We can also test that adding -b to a is the
same computing as a - b. And finally we can test chaining of neg, as always.

I wonder what the most interesting program is that we could implement on top
of what we have so far. Tomorrow we add a few more convenience functions
before we start heading into the more interesting stuff.

I think I may have solved the test framework problem. More on that when we
get to v0.11 and v0.12.

There is a github branch here v0.4 for this article.

Previous article: v0.3 - copy, zero, normalise, mul1

Next article: v0.5 - add, sub, cmp

BSDNT - v0.3 copy, zero, normalise, mul1

2010-09-09T13:54:00.000-07:00

In this section we add a few more trivial convenience functions, nn_copy,

nn_zero and nn_normalise.

The first two are simple for loops which we implement in nn.h.

A pair {c, m} will be considered normalised if either m is zero

(representing the bignum 0) or the limb c[m-1] is nonzero. In other words,

the most significant word of {c, m} will be non-zero.

The only tricky thing to do with nn_normalise is to make sure nothing goes

wrong if all the limbs are zero.

We add numerous tests of these new functions. We zero an integer and check

that it normalises to zero limbs. We copy an integer and make sure it

copies correctly. We also copy an integer, then modify a random limb and

then check that it is no longer equal. This provides a further check that

the nn_equal function is working correctly.

Finally, we do a more thorough test of nn_normalise by zeroing the top few

limbs of an integer then normalising it. We then copy just this many limbs

into a location which has been zeroed and check that the new integer is

still equal to the original *unnormalised* integer. This checks that

nn_normalise hasn't thrown away any nonzero limbs.

Next we add a mul_1 function. As for all of the linear functions we have

been adding, this will be very slow in C. This time we need to retrieve

the top limb of a word-by-word multiply. Again, ANSI C doesn't provide this

functionality, so we use the gcc extensions that allow us to upcast to a

double word before doing the multiply. Again the compiler is supposed to

"do the right thing" with this combination.

We add a test to check that a * (c1 + c2) = a * c1 + a * c2. We also add

a test for chaining of mul1's together, passing the carry-out of one to

the carry-in of the next.

The github branch for this post is here: v0.3

Previous article: v0.2 - subtraction, shifting

Next article: v0.4 - add1, sub1, neg, not

BSDNT - v0.2 subtraction, shifting

2010-09-08T08:00:00.000-07:00

In today's revision we add a subtraction function and left and right

shift functions.

Subtraction is almost trivial after addition. We copy what we did for

addition. The main change is that carry-ins become borrows. we have to

negate our borrows at each step, so that we are always dealing with

positive quantities. The returned borrow is zero or positive.

In the test code for subtraction, we can add a test which checks that

a + b - b = a. This is a further test of the addition function as well.

Left shift is straightforward. Again we make use of the dword_t type

which is twice as wide as a word_t. This allows us to do a single shift

(which the compiler has the option of converting to two shifts if required).

We shift each limb by the specified bits, then hang onto the bits that were

shifted out the top, to be added back in at the next iteration.

For right shift we only have three functions instead of four. It doesn't make

sense to write out the "carry-out" at the bottom of a right shift. Instead

we do a "read-in" nn_r_rsh version which is the exact opposite of nn_s_lsh. It

reads the "carry-in" from an (m + 1)-th limb (which it shifts appropriately).

The nn_r_rsh function, also returns the bits that are shifted out the bottom end.

We can provide an extra test for addition now, by checking that a + a = a << 1

We also check that (a << sh1) >> sh1 gets us back to a.

Finally, we test that we can chain our carry-in carry-out functions together. For

example, I should be able to shift the first few limbs of an integer, feed the

carry-out as a carry-in in to the next shift, which will then shift the rest of

the integer. This should give me the same result as if I had shifted the whole

integer in one go.

We also do similar carry-in/carry-out chaining with the addition and subtraction

functions.

We are going to need to set up a proper test framework at some point, but for

now, whilst we are still working on implementing "linear" functions, I'm simply

going to keep duplicating the test functions over and over and making the

minor modifications required to test the new functions we are adding.

Here is v0.2 on github for this post.

Previous article: v0.1 - basic types and addition

Next article: v0.3 - copy, zero, normalise, mul1

BSDNT v0.1 - basic types and addition

2010-09-07T11:56:00.000-07:00

We start with v0.1 of bsdnt. See the github branch:

http://github.com/wbhart/bsdnt/tree/v0.1

Firstly, we set up a bit of infrastructure before we add our first function,

the addition function.

We begin that with some basic types in nn.h:

word_t - represents a machine word (either 32 or 64 bit)

dword_t - twice the size of a machine word (to handle carries from addition)

nn_t - our basic multiprecision integer type, consisting of an array of words

nn_src_t - same as an nn_t, but declared const, so it can't be modified, used

for input/source parameters

len_t - the length in words of a multiprecision integer. We don't use a struct

as we'd have to pass it by reference, which would be less efficient due

to having to dereference the pointer all the time

rand_t - a random state, unused for now. This will ensure thread safety of our

random functions later on.

We also define WORD_BITS, the number of bits in a machine word. To simplify

things, we'll informally let B refer to the number 2^WORD_BITS in what follows.

Our basic bignum type is a pair {c, len} consisting of an nn_t c and a len_t len

counting the number of limbs of our bignum. At this stage we restrict len to

being non-negative, and we allow leading zero limbs in our representation of

our bignums. This allows for twos complement arithmetic with fixed length

bignums (arithmetic modulo B^len for some fixed len).

The first function we write is a simple addition function. This needs to be

pretty sophisticated. In particular, we want to be able to pass in a carry, so

that addition functions can be chained together, or on the end of other functions.

The first addition function will add two operands of the same length (number of

machine words), which we signify with an m at the end of the function name.

We signify that the function takes a carry-in by appending c to the function name.

We will also want the function to pass a carry-out. It'll return this carry out as a

word_t.

The function nn_add_mc is defined in nn.c. Notice how we cast each word to a

dword first, then do the addition, then cast back to a word for the low word of

the sum, and shift by WORD_BITS to get the high word of the result.

The compiler is *supposed* to optimise this combination. Actually, it makes poor

use of the processor carry flag, and screws up the loop, which is why C is not

the best language for a bignum library. But for now it is the best we can do.

Later on we'll add assembly optimisations to our library.

We could unroll the loop (write the contents four times in a row, say, to amortise

the cost of the loop counter arithmetic over four iterations). However, we are

after a cleanly coded library, so we resist this temptation until we write some

assembly code.

Now we introduce some extra macros in nn.h for our addition function. These

represent various permutations of allowing a carry in or not, and one more

interesting macro, nn_s_add_m. This function checks whether the result

and first input are aliased. If so (in other words, we have a = a + c), then the

carry-out gets *added* to the correct limb of a.

If a and b are not aliased (we have a = b + c), then the carry is *written*

and not added, to the correct limb of a. This macro will make coding cleaner

later on.

Notionally, the extra 's' in the function name stands for 'store'.

We use macros for these functions to stop code duplication, and because macros

prevent extra function call overhead, and sometimes offer an opportunity for the

compiler to optimise further.

Note that C macros are not hygienic. In fact they are just text replacement macros.

Thus we need to be careful about naming the macro variables so they are unlikely

to conflict with symbols passed in by the caller. In this case we don't introduce

any variables inside our macros, so this precaution is not necessary.

We also need to take care when using the macro parameters inside the macro. In

some cases it is necessary to put parentheses around the parameters in case

the caller passes in a complex expression which combines badly with expressions

inside our macro (e.g. due to operator precedence). Where there can be no

confusion when substituting a macro, we do not need to do this.

Next we will need to add some random functions, to generate random integers to add

in our test code. In nn.c, we add functions to generate a random word of data, a

random integer up to a limit and a random multiprecision integer.

The random functions take a random state as a parameter. This is essential to

ensure that we can make our random functions threadsafe at a later date. For now

the random state and associated init and clear functions do nothing.

The random word function generates two random integers slightly bigger than half

a word and stitches them together. This is done using an el cheapo pseudorandom

function which works modulo a prime slightly bigger than half a word. We take some

care not to always end up with even output, or output divisible by a given small

prime.

The test code is in t-nn.c. We add a comparison function nn_equal_m to nn.h, to

test equality of two mp integers, to be used in our test code, and some other

convenience functions. We generate random length mp integers, and check an

identity, in this case the associative law for addition. This is not a very

sophisticated test, and we'll have to improve out test code at a later date.

For now, it passes, and we move on to grander things.

Previous article: BSDNT-introduction

Next article: v0.2 - subtraction, shifting

BSDNT - Introduction

2010-09-07T11:38:00.001-07:00

Over the next few weeks to months, I'll be blogging about a new project of mine, called bsdnt.

The aim of this project is to develop a cleanly coded bignum library with (eventually) reasonable performance and a BSD license. I'll be developing this code for a while in a blog/tutorial style.

First some organisational matters.

For those who wish to follow along with the code as it is written, it can be found at:

http://github.com/wbhart/bsdnt/

If you have a github account, you can clone my project on github.

To make a clone of my repository on your local machine, you can do:

git clone git://github.com/wbhart/bsdnt.git bsdnt

The various branches, as I add them, will be called v0.1, v0.2, etc.

To switch branches, simply do, for example:

git checkout origin/v0.1

To make your own branch of this, to mess around in, do:

git checkout -b mybranch

Antony Vennard, set up a google groups discussion list. You can access that

here:

http://groups.google.com/group/bsdnt-devel?hl=en

A word on licensing: if you substantially copy my code, then I appreciate you

retaining the copyright, however if you write substantially your own code,

merely being "inspired" by my code, I am fine with you not adding my copyright to your project. But still let me know about your code, because I'd be interested in seeing it.

To compile the project, you will need a recent version of gcc. I'm using

gcc 4.4.1, and some of the features may not work with earlier than gcc 4.4.

You can proceed to the first blog about the code here:

Next article: BSDNT v0.1 - basic types and addition

Stacks and Elliptic Cohomology

2009-06-21T10:18:00.000-07:00

This week I became interested in two different topics due to conversations that I overheard. The first is the topic of stacks and the second is elliptic cohomology.

Stacks

Apparently there are numerous different kinds of stacks - Deligne-Mumford Stacks, Artin Stacks and for the die hard, apparently more general kinds of stacks.

The following is probably completely wrong, but is my understanding of what stacks are about.

Consider elliptic curves defined over the complex numbers K = C say. It is a classical result that up to isomorphism, these can be parameterised by points in the complex upper half plane, modulo the action of the modular group.

Now the upper half plane, modulo the modular group can be compactified by adding a point at the cusp, and made into a Riemann surface (of genus zero in this case). We can put coordinates on this Riemann surface (the complex j-function as it happens) and turn it into an algebraic curve of genus 0.

In fact, two elliptic curves are isomorphic iff their j-invariants are equal. In other words the Riemann sphere or j-line as it is often called, can be thought of as classifying all elliptic curves. In fact we call C the course moduli space for elliptic curves.

Now suppose we try to construct a universal space for families of elliptic curves over this base space. The problem is that an elliptic curve can have extra automorphisms and it is possible for a family of elliptic curves to contain isomorphic elliptic curves for this reason. That prevents us from having a univeral space for families of elliptic curves.

The way we get around this is using stacks. We define the stack of elliptic curves which is a category whose objects are families of elliptic curves over a base space (fixed for that family) and we define a morphism to be a map between families of elliptic curves along with a map between the corresponding base spaces such that the map between families is compatible with the map between base spaces.

Furthermore, for it to be a morphism (X'->B') -> (X->B), we require that if we pull the family of elliptic curves X back along the map B'->B of base spaces, we get a family of elliptic curves isomorphic to the family X'.

We can restrict to the subcategory of families of elliptic curves over a fixed base space B if we want. We call this subcategory the fibre over B.

Now note that the fibre over a base space B is a groupoid (a categories whose only morphisms are isomorphisms). We say that the stack of elliptic curves is a category fibred in groupoids.

Now it is clear that there is a universal family of elliptic curves with respect to this construction.

There's more to stacks than this (some of the critical components of the definition are omitted above).

A stack is of Deligne-Mumford type (formerly an algebraic stack) if it satisfies some additional conditions, in particular that there is an etale surjective morphism (called an atlas) from a scheme U to the stack F, amongst other things.

An Artin stack (nowadays what is referred to as an algebraic stack) simply replaces etale with smooth in the previous definition.

Anyhow, what was interesting to me is that nowadays stacks are replacing schemes as the ultimate objects of interest. A lot of work has been done to popularise them. Hey, if you want to know more theres only 1000 pages to read: The Stacks Project.

Elliptic cohomology

The second topic which piqued my interest this week was elliptic cohomology. I thought maybe this might be related to parabolic cohomology, which is defined in terms of parabolic cusps (fixed points of parabolic elements of SL_2(R)). But I don't know that this is the origin of the term.

Instead I found this enormous survey on the web, which is written helpfully:

http://www-math.mit.edu/~lurie/papers/survey.pdf

That's a long document to be called a summary, so I'll give a summary of the summary.

To a topological space X we can associate the singular cohomology groups A^n(X) = H^n(X; Z) which can be characterised by the Eilenberg-Steenrod axioms. Any collection of functors and connecting maps satisfying these axioms necessarily gives you the usual integral cohomology functors (X \subseteq Y) -> H^n(X, Y; Z). More generally we can replace Z with any abelian group M.

Now if we drop the last of the E-S axioms:

* If X is a point then A^n(X) = { 0 if n \neq 0 and Z is n = 0 }

then we get something more general, called a cohomology theory.

Interestingly, complex K theory is an example of such a cohomology theory!

Complex K-theory is a so-called multiplicative cohomology theory, because A^n(*) is a graded commutative ring. Another nice feature is it is periodic:

* According to the Bott periodicity theorem for complex K-theory, there is an element \beta in K^2(*) such that multiplication by \beta induces an isomorphisms: \beta : K^n(X) -> K^{n+2}(X)

and even:

* K^i(*) = 0 if i is odd.

Ordinary cohomology H*(X; A) for a commutative ring A is even, but to make it periodic we need to take products over every second cohomology group A^n(X; A) = \prod_k H^{n+2k}(X; A)

Now when A is an even periodic cohomology, it turns out that A(CP^\infty) of the infinite dimensional complex projective space, is isomorphic to a formal power series ring A(*)[[t]] over the commutative ring A(*).

We can view the element t as the first Chern class of the universal line bundle O(1). The space CP^\infty is a classifying space for complex line bundles, i.e. for any complex line bundle L on a space X there is a (classifying) map \phi : X -> CP^\infty and an isomorphism L <-> \phi*(O(1)).

We then define c1(L) = \phi*(t) \in A(X), the first Chern class of the cohomology A.

In ordinary cohomology the first Chern class of a tensor product of line bundles is simply the sum of the first Chern classes of the line bundles.

In the case of complex K-theory line bundles L can be thought of as representatives of elements of K(X) itself. We write such an element [L]. Then c1(L) = [L] - 1.

Now c1(L1 \tensor L2) = c1(L1) + c1(L2) + c1(L1)c1(L2).

In the general case we have:

c1(L1 \tensor L2) = f(c1(L1), c1(L2)) for some f \in A(*)[[t1, t2]].

It turns out that the following properties hold:

* f(0, t) = f(t, 0) = t

* f(u, v) = f(v, u)

* f(a, f(b, c)) = f(f(a, b), c)

A power series with these properties is called a commutative 1-dimensional formal group law over the commutative ring A(*). It gives a group structure on the formal scheme Spf A*(X)[[t]] = Spf A(CP^\infty).

We call the first of the group laws above f(a, b) = a + b the additive formal group law denoted \hat{G_a} and the other f(a, b) = a + b + a*b the multiplicative formal group law, denoted \hat{G_m}.

Now one wonders if there are other possible formal group laws. It turns out that the Lazard ring is a ring classifying formal group laws. Quillen proved that it comes from a cohomology called periodic complex cobordism denoted MP. There is a canonical isomorphism MP(CP^\infty) <-> MP(*)[[t]] and the coefficient ring MP(*) is the Lazard ring.

One may construct the moduli stack of all formal group laws M_{FGL} so that for a commutative ring R, the set of homomorphisms Hom(Spec R, M_{FGL}) can be identified with the power series f(u, v) \in R[[u, v]] satisfying the three conditions above. Then M_{FGL} is an affine scheme, in fact M_{FGL} = Spec MP(*), as we'd expect for a moduli stack.

Well it gets a bit more complicated than that. One must mod out by the action on M_{FGL} by the group of automorphisms of the formal affine line Spf Z[[x]]. This yields the stack of formal groups M_{FG}.

Now what is an elliptic cohomology?

Well we need to relax one thing slightly. Instead of demanding that we have a periodic cohomology, we'll just require our multiplicative cohomology A to be weakly-periodic:

* The natural map A^2(*) \tensor_{A(*)} A^n(*) -> A^{n+2}(*) is an isomorphism for all n \in Z.

Now we can define an elliptic cohomology A:

* R is a commutative ring

* E is an elliptic curve over R

* A is a multiplicative cohomology which is even and weakly-perioidic

* There are isomorphisms A(*) <-> R and \hat{E} <-> Spf A(CP^\infty) of formal groups, over R which is isomorphic to A(*)

Here \hat{E} represents the formal completion of E along its identity section. It is a commutative 1-dimensional formal group over R. It is classified by a map \phi : Spec R -> M_{FG}.

For finite cell complexes X we may interpret the complex cobordism groups MP^n(X) as quasi-coherent sheaves on the moduli stack M_{FG} and in the case of a formal group over R we can define A^n(X) to be the pullback of these sheaves along \phi.

When \phi is flat then we call the formal group Landweber-exact. Landweber gave a criterion for determining when a formal group is Landweber-exact.

Anyhow it turns out that in the elliptic cohomology case, when the formal group is Landweber-exact, the elliptic curve and the isomorphism of the definition of elliptic cohomology are uniquely given. In this case, giving the elliptic cohomology theory is exactly the same thing as giving an elliptic curve.

Now in the case of the formal multiplicative group, there is a universal formal group. But in the case of elliptic cohomology, there is no such thing as a universal elliptic curve over a commutative ring. The moduli stack of elliptic curves M_{1,1} is not an affine variety, and not even a scheme. However, it is a Deligne-Mumford stack.

For each etale morphism \phi : R \to M_{1,1} there is an elliptic curve E_\phi and happily, when \phi is etale, the associated formal group \hat{E_\phi} is Landweber-exact.

Well, this allows us to create on M_{1,1} a presheaf taking on values in the category of cohomology theories. But this is a nasty construction and so we try to represent our cohomology theories by representing spaces.

Roughly speaking this is where the theory of E-\infty rings or E-\infty spectra comes in. Basically a cohomology theory A_\phi can be represented by a spectrum. That eventually allows one to develop a universal elliptic cohomology theory.

We define a presheaf O_{M}^{Der} of E_\infty rings representing elliptic cohomology theories on the category {\phi : Spec R -> M_{1,1}} mentioned above.

Now we extract our universal cohomology theory by taking a "homotopy limit" of the functor O_{M}^{Der}. This gives us the E-\infty ring of topological modular forms tmf[\Delta^{-1}].

(Note tmf is what you get when you replace M_{1,1} by its Deligne-Mumford compactification. After inverting 2 and 3 [to do with the extra automorphisms on elliptic curves] there is an isomorphism from tmf to the ring of integral holomorphic modular forms.)

OK that summarises the first 9 pages of the long summary above. But probably it isn't much of an improvement on the original. But it helps the memory to write it all down somewhere.

MPIR - version 1.2

2009-06-07T11:25:00.000-07:00

Finally version 1.2 of MPIR (Multiple Precision Integers and Rationals), of which I am a lead developer, is released: http://www.mpir.org/

MPIR is an open source project based on the GNU Multi Precision (GMP, see http://www.gmplib.org/) library, but still licensed under version 2 of the LGPL.

About a month ago version GMP 4.3.0 was released, which they had been preparing for a LONG time. We expected some nice features, and found some, which we have been subsequently catching up with.

In particular we needed:

* Faster assembly code for multiplication basecase

* Faster unbalanced integer multiplication (where you are multiplying integers of different sizes)

* Improvements to the speed of multiplying medium sized integers (50-2000 words where 1 word = 2^64 on a 64 bit machine)

* Asymptotically fast extended GCD

* Faster integer multiplication for large integers (> 2000 limbs)

* Faster integer squaring

* Other assembly improvements

Multiplication Basecase

Jason Moxham, a brilliant MPIR developer from the UK decided to take on the mul_basecase challenge. He's been writing an assembly optimiser for quite a few months. It takes hand written assembly code and reorganises the instructions over and over, within permitted bounds, to try and find an optimal sequence.

The results over the past year are pretty impressive to see:

Multiplications per second
	128 x 128 bits	512 x 512 bits	8192 x 8192 bits
MPIR 1.2.0	53794646	12488043	117404
GMP 4.3.0	52766506	10879150	114927
MPIR 1.1.2	51802252	11802334	111772
MPIR 1.0.0	35856598	10928085	111641
MPIR 0.9.0	37299412	8122452	86301
GMP 4.2.1	25896136	6383542	60819

Note that the number of multiplications that can be done per second has more than doubled in most cases since last year, and all this just from assembly language improvements.

The timings for this table were made on an Opteron (AMD K8) server, running at 2.8 GHz.

Toom Multiplication

Toom multiplication can be described in terms of a polynomial multiplication problem. Firstly the large integers to be multiplied are split apart into pieces, which form the coefficients of polynomials. Then the original integer multiplication problem becomes one of polynomial multiplication, then a reconstruction phase at the end, where the polynomial coefficients of the product are stitched together to make the product integer.

This can be thought of in terms of writing the original integer in terms of some large base, 2^M, where M is usually a multiple of the machine word size (e.g. M = 64B on a 64 bit machine).

For example, Toom-3 breaks the two integers into 3 pieces each, corresponding to the digits in base 2^M:

A = a0 + a1 * 2^M + a2 * 2^2M

B = b0 + b1 * 2^M + b2 * 2^2M

So we construct two polynomials:

f(x) = a0 + a1 * x + a2 * x^2

g(x) = b0 + b1 * x + b2 * x^2

Then A * B = f(2^M) * g(2^M). So we first compute h(x) = f(x) * g(x), and then A * B = h(2^M).

To multiply the polynomials f(x) and g(x), we note that their product will be degree 4, and so we can determine it fully if we know its value at 5 independent points. We choose, for convenience, the points 0, 1, -1, 2, infinity. Finally we note that h(0) = f(0) * g(0), h(1) = f(1) * g(1), etc.

Thus to compute the product we only need to find the values of f(x) and g(x) at the chosen points, do five small multiplications to get the value of h(x) at those points, then interpolate to get the coefficients of h(x).

So we have replaced the original large multiplication with five small ones. Note that each of the small multiplications involves integers at most one third of the size. If we just used schoolboy multiplication to multiply the "digits" of our large integers, we'd do 3 x 3 = 9 small "digit multiplications". Instead, through the magic of evaluation/interpolation we only have 5 small "digit multiplications" to do (and a few additions and subtractions, etc., for the evaluation and interpolation phases).

The basecase multiplication code already handles unbalanced multiplication, so I decided to focus on unbalanced variants of Toom multiplication algorithms.

In MPIR we use Toom-2 (also known as Karatsuba multiplication - though we use a variant called Knuth multiplication where evaluation happens at 0, -1 and infinity), Toom-3, Toom-4 and Toom-7.

I spent a fair bit of time optimising Toom-3. As we have a lot of new assembly instructions available in MPIR, I was actually able to get the interpolation sequence used down from about 11 to 8 steps in MPIR 1.2. I also discovered that it was possible to use the output integer for temporary space. If one is careful, one can even set it up so that the temporaries stored in the output space don't have to be moved at the end, i.e. they are in precisely the right place as part of the output integer at the end of the algorithm.

I made similar optimisations for Toom-4 and Toom-7 and I also switched the interpolation phase of the algorithms over to twos complement.

Here are the results of this work in the Toom range:

Toom Multiplication
	Kara (3200 bits)	Toom3 (7680 bits)	Toom4 (25600 bits)	Toom7 (131072 bits)
MPIR 1.2.0	300414	74337	11428	1153
GMP 4.3.1	274738	68498	11263	1042
MPIR 1.1.2	260501	60039	9890	993
MPIR 1.0.0	261599	62275	9900	826
GMP 4.2.1	163273	33460	4980	408

Timings are on a 2.4GHz Core 2 (eno).

The differences in timing in the karatsuba region for MPIR give an indication of how much of a speedup is occurring on account of better assembly code. Anything beyond that indicates a speedup in the Toom algorithms themselves.

Note GMP 4.2.1 had Karatsuba and Toom-3 only, GMP 4.3.1 does not have Toom-7.

Unbalanced Multiplication

Now unbalanced multiplication proceeds in the same way, except that we have integers, and thus polynomials, of different length:

f(x) = a0 + a1 * x + a2 * x^2 + a3 * x^3

g(x) = b0 + b1 * x

Again the product is degree four and so we need to interpolate at five points, which we can choose to be precisely the same points that we used for Toom-3, even reusing the interpolation code, if we wish.

We call this Toom-42, denoting that we split our integers into 4 and 2 "digits" respectively, in our base 2^M.

The result is 5 small multiplications instead of 4 x 2 = 8; still a substantial saving.

Finally I implemented some of the unbalanced variants. In particular we now have Toom-42, an unbalanced version of Toom-33 (where the top coefficients are not necessarily exactly the same size) and Toom-32.

The results for unbalanced multiplications in the Toom range are now quite good:

Unbalanced Multiplication
	15000 x 10000 bits	20000 x 10000 bits	30000 x 10000 bits
MPIR 1.2.0	32975	25239	14995
GMP 4.3.1	31201	24099	14104
MPIR 1.1.2	23190	19970	13289
GMP 4.2.1	13235	10855	7235

Timings are on the 2.4GHz Core 2 (eno).

Toom Squaring

Squaring using the Toom algorithms is much the same, except that there is no need to evaluate the same polynomial twice. We also save because the pointwise multiplications are now also squares and we can recurse, right down to a basecase squaring case.

All of the same optimisations can be applied for Toom squaring as for ordinary Toom multiplication. Here are the comparisons:

Toom Squaring
	Kara (5120 bits)	Toom3 (12800 bits)	Toom4 (51200 bits)	Toom7 (131072 bits)
MPIR 1.2.0	358656	82274	10514	2762
GMP 4.3.0	357031	83676	10430	2564
MPIR 1.1.2	349686	78510	9917	2347
GMP 4.2.1	185201	42256	5112	1237

The performance is now comparable to GMP until the Toom7 region, where we pull ahead considerably. These timings were done on the 2.8GHz Opteron (K8) server.

FFT Multiplication

For multiplication of large integers, MPIR uses a Fast Fourier Transform (FFT) method. Instead of working over the complex numbers, one can work over a ring Z/pZ where p is a special prime, e.g. p = 2^M + 1 where M is some power of 2. This trick is due to Schonhage and Strassen.

For MPIR 1.2 we decided to switch over to the new FFT of Pierrick Gaudry, Alexander Kruppa and Paul Zimmermann. An implementation of ideas they presented at ISAAC 2007 was available to download from Paul and Alex's websites.

Most of the work in merging this FFT was removing bugs, making the code work efficiently on Windows and writing and running tuning code for all the platforms we support.

I wrote the tuning code, Brian Gladman discovered some primitives to use on Windows which replace inline assembler available on Linux and Jason Moxham, Brian Gladman, Jeff Gilchrist and myself tuned the FFT on a range of systems.

Here are some examples of the speedup:

FFT Multiplication
	2000000 bits	6000000 bits	20000000 bits	64000000 bits
MPIR 1.2.0	74.9	12.8	5.20	1.47
GMP 4.3.0	52.4	13.4	3.66	0.813
MPIR 1.1.2	47.2	12.5	3.09	0.742
GMP 4.2.1	32.8	8.66	2.11	0.528

Timings are again on the 2.8GHz Opteron server.

Extended GCD

A longstanding issue with MPIR has been the lack of fast extended GCD. We had merged Niels Mollers asymptotically fast GCD code, but unfortunately no suitable extended GCD implementation was available to us.

I began by reading Niels' paper to understand the half-GCD algorithm (ngcd) that he invented. Essentially half-GCD algorithms get their asymptotic improvement over the ordinary Euclidean algorithm as follows:

Begin by noting that the first few steps only depend on the uppermost bits of the original integers. Thus instead of working to full precision, one can split off the topmost bits and compute the first few steps of the Euclidean GCD on those bits. One keeps track of the steps taken in matrix form. One then has an exact representation for the steps taken so far in the algorithm.

Next one applies the matrix to the original integers. The integers that result are smaller than the original integers. One then finishes off the algorithm by applying the usual steps of the GCD algorithm to the smaller integers.

Of course to get an asymptotic improvement one needs to take care how to apply the matrix and to recurse the whole algorithm down to some fast basecase.

Once I understood that Niels' implementation explicitly computed the matrix at each step, it occurred to me that in order to turn it into an extended GCD algorithm, all I had to do was apply the matrix to the cofactors and keep track of any other changes that would affect the cofactors, throughout the algorithm.

As an example of the speedup, for 1048576 bit integers extended GCD benches as follows:

MPIR 1.1.2: 0.453

MPIR 1.2.0: 4.06

Again the benchmarks are for the 2.8GHz Opteron.

Assembly Improvements

Finally this mammoth MPIR release also contained numerous sped up assembly routines due to Jason Moxham for AMD Opterons (which usually also translate to improvements for other 64 bit x86 platforms).

Below I include timings (smaller is better) at each point where improvements have been made in the assembly routines in MPIR.

Assembly Improvements
	MPIR 0.9	MPIR 1.0	MPIR 1.1.2	MPIR 1.2
mpn_add_n	1598	1524
mpn_sub_n	1598	1524
mpn_xor_n	3015	2271
mpn_and_n	3014		1772	1523
mpn_ior_n	3014		1771	1525
mpn_nand_n	3013		2035	1775
mpn_nior_n	3013	1773
mpn_andn_n	3013	2273
mpn_iorn_n	3015	2271
mpn_addadd_n		2525		2190
mpn_addsub_n		2526		2189
mpn_subadd_n			2526	2190
mpn_addmul_1	3094	2524
mpn_submul_1	3092	2524
mpn_mul_1	3024	2522
mpn_sublsh1_n		2527	2400	2190
mpn_com_n	3014	1271
mpn_divexact_by3	12016	2278
mpn_lshift	2531	1701
mpn_rshift	2532	1617
mpn_hamdist	8281	1791
mpn_popcount	7281	1527
MPN_COPY	3012		2017	1021
mpn_divrem_1	23649		15119
MPN_ZERO	2013			783

This time timings are taken on a 2.6GHz AMD K10 box.

Numerous new assembly functions have also been added for 64 bit x86 machines since MPIR began, including: mpn_addmul_2, mpn_addadd_n, mpn_sublsh1_n, mpn_divexact_byff, mpn_rsh1add_n, mpn_lshift1, mpn_rshift1, mpn_rsh1sub_n, mpn_mul_2, mpn_lshift2, mpn_rshift2.