BSDNT - configure and assembly improvements

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

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

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

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

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

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

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