Article 39152 of comp.sys.dec:
johnl@Radix.Net (John A. Limpert) writes:

> On Thu, 16 Jan 1997 19:37:54 -0800, Jeffrey Potoff
> <jpotoff@latoso.cheme.cornell.edu> wrote:
> 
> >What hype would that be ?  Facts are facts.  Send whatever source 
> >code you're running that's so much faster on a PPro and I'll 
> >benchmark it on a 433mhz 21164 Alpha so we can end this nonsense once
> >and for all.
> 
> I benchmarked an Alpha with a Reed-Solomon encoder that I wrote and my
> 90 MHz Pentium beat it. The Alpha seemed to have alot of trouble with
> all of the byte sized operations in the encoder software. The Alpha is
> a neat chip but it stumbled on my application.

John was kind enough to send me his Reed Solomon code.  Sure enough, I
got the following results on a 200MHz P6 Linux machine and two 300MHz
ev5 machines (21164 without subword load/store).  One of the ev5
machines was running Linux, the other OSF/1.  Under Linux, I used gcc
-O2 for compilation, under OSF/1 I used cc -migrate -O4 -tune ev5:

original code:

686:200MHz:Linux:gcc -O2:			2.8MB/sec
ev5:300MHz:Linux:gcc -O2:			0.74MB/sec
ev5:300MHz:OSF/1:cc -migrate -O4 -tune ev5:	0.74MB/sec

The code had some obvious inefficiencies, so I spent half an hour
getting rid of the worse offender.  The results with that code are
below:

improved code:

686:200MHz:Linux:gcc -O2:			 3.4MB/sec
ev5:300MHz:Linux:gcc -O2:			10.0MB/sec
ev5:300MHz:OSF/1:cc -migrate -O4 -tune ev5:	11.4MB/sec

This is obviously much more reasonable.  Note that the 686 improved as
well, but only by 21% as compared to a factor of 15.4 for the Alpha.
I suspect there are two reasons the Alpha did so poorly with the old
code: (a) some of the arrays are exactly 8KB big (size of primary
cache), so there might be some thrashing going on and (b) the byte
accesses certainly hurt on the pre-ev56 CPUs with the original code.

Anyhow, spending half an hour to get a 21% improvement seems pretty
reasonable to me and on the Alpha it happens to win big.  The improved
code is appended below.  I included the original core-loop code as a
comments, so it's easier to see what changes I made.

If you look at the code, you'll notice that it takes advantage of the
64 bit registers on the Alpha.  Since some people might argue it is
not worth spending extra time optimizing code for the Alpha, I also
measured the Alphas when using exactly the same C source as for the
686.  The results for this case are below:

ev5:300MHz:Linux:gcc -O2:			6.0MB/sec
ev5:300MHz:OSF/1:cc -migrate -O4 -tune ev5:	5.9MB/sec

So the larger registers definitively help (though the improvement may
well be due to reducing the number of memory accesses as opposed to
doing more work per arithmetic operation).

Note that I literally spent just half an hour getting rid of the bad
core-loop code.  I'm sure the code could be tuned much more for both
CPUs, but that's not really the point of this exercise.

Lessons:

	(a) The absence of subword load/stores can really hurt
	    on POORLY written code.

	(b) Just making minor changes, performance increased
	    by over 20% for the 686 and 1400% on the Alpha.
	    The point here is that the same changes result
	    in benefits for both architectures (though the
	    benefit is smaller in the 686 case).

	(c) The use of 64 bit register on the Alpha really
	    helps in the above code.

  --david
--
/*
 * rse.cpp
 * Version 1.0
 * 05/09/95
 *
 * Reed-Solomon encoder for the standard NASA/CCSDS (255,223) code.
 *
 *     The Berlekamp Reed-Solomon encoder is implemented as
 *     a conventional Reed-Solomon encoder with dual/conventional
 *     basis tranforms on the input and output of the encoder.
 *
 *     Interleave (1..8) and virtual fill are supported.
 *
 * References:
 *
 *     CCSDS Blue Book 101.0-B-3
 *     Telemetry Channel Coding
 *     May 1993
 *
 *     Error Control Coding: Fundamentals and Applications
 *     Shu Lin & Daniel Costello
 *
 *     Reed-Solomon Encoders - Conventional vs. Berlekamp's Architecture
 *     M. Perlman & J. Lee
 *     JPL Publication 82-71
 *
 * Author:
 *
 *     John A. Limpert
 *     ATSC
 *     Code 530.2
 *     NASA/GSFC
 *     Greenbelt, Maryland 20771
 */

#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <assert.h>

#include <sys/time.h>

/*
  Sample encoder output:

  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
  00 00 00 00 00 00 00 00 00 00 00 00 00 00 FF 22
  07 1D 51 81 3F F6 86 E7 5D F5 64 1D 37 B8 D7 B8
  37 1D 64 F5 5D E7 86 F6 3F 81 51 1D 07 22 FF
*/

/*
  Conventional to Dual Basis Transform
  
  Tal
  
  1 0 0 0 1 1 0 1
  1 1 1 0 1 1 1 1
  1 1 1 0 1 1 0 0
  1 0 0 0 0 1 1 0
  1 1 1 1 1 0 1 0
  1 0 0 1 1 0 0 1
  1 0 1 0 1 1 1 1
  0 1 1 1 1 0 1 1
*/

static unsigned char tal[] = {
    0xFE, 0x69, 0x6B, 0x0D, 0xEF, 0xF2, 0x5B, 0xC7
};

/*
  Dual to Conventional Basis Transform

  Tal -1

  1 1 0 0 0 1 0 1
  0 1 0 0 0 0 1 0
  0 0 1 0 1 1 1 0
  1 1 1 1 1 1 0 1
  1 1 1 1 0 0 0 0
  0 1 1 1 1 0 0 1
  1 0 1 0 1 1 0 0
  1 1 0 0 1 1 0 0
*/

static unsigned char ital[] = {
    0x9B, 0xDD, 0x3E, 0x1C, 0x37, 0xB3, 0x60, 0x94
};

/* F(x) = x^8 + x^7 + x^2 + x + 1 over GF(2) */
const unsigned char poly = 0x87;	/* field generator polynomial */

#define N 255	/* number of symbols in codeword */
#define K 223	/* number of information symbols in codeword */
#define E 16	/* symbol error correction capability */

void rsInit();
void rsEncode(unsigned char* codeBlock, int interleave, int virtualFill);

typedef union {
    unsigned char u8[32];
    unsigned int  u32[8];
    unsigned long u64[4];
} big_word;

static unsigned char convToDual[256];		/* conventional to dual */
static unsigned char dualToConv[256];		/* dual to conventional */
static unsigned char alogAlpha[256];		/* GF(2^8) alog */
static unsigned char logAlpha[256];		/* GF(2^8) log */
static unsigned char mulAlpha[256][256];	/* GF(2^8) multiplication */
static unsigned char g[32];			/* RS G(x) polynomial */
static big_word products[256];			/* precomputed products */














#ifdef TEST
static unsigned char codeBlock[(N - 95)*8];

int
main()
{
    struct timeval tv_start, tv_stop;
    double delta;
    int i;

    /* test data */
    memset(codeBlock, 0, sizeof(codeBlock));
    codeBlock[1016] = 0xFF;
    codeBlock[1017] = 0xFF;
    codeBlock[1018] = 0xFF;
    codeBlock[1019] = 0xFF;
    codeBlock[1020] = 0xFF;
    codeBlock[1021] = 0xFF;
    codeBlock[1022] = 0xFF;
    codeBlock[1023] = 0xFF;
    
    /* initialize encoder tables */
    rsInit();

    /* test/benchmark encoder */
    gettimeofday(&tv_start, 0);
    for (i = 0; i < 10000; ++i)
        rsEncode(codeBlock, 8, 95);
    gettimeofday(&tv_stop, 0);
    delta =
	(tv_stop.tv_sec + tv_stop.tv_usec / 1000000.0)
	- (tv_start.tv_sec + tv_start.tv_usec / 1000000.0);
    printf("%g seconds/10000 iterations => %gMB/sec\n",
	   delta, 10000.0*sizeof(codeBlock)/delta/1024.0/1024.0);
    
#if 1
    /* display encoded codeblock */
    for (i = 0; i < sizeof(codeBlock); ++i) {
        printf("%.2X", codeBlock[i]);
        if (i%16 == 15)
            printf("\n");
        else
            printf(" ");
    }
    printf("\n");
#endif
    
    return 0;
}
#endif /* TEST */










/* rsInit -- initialize Reed-Solomon encoder */

void
rsInit()
{
    int i, j, k;
    unsigned char x, t, r;
    
    /* create tables for conventional/dual basis conversion */
    for (i = 0; i <= 255; ++i) {
	x = i;
	r = 0;
	/* multiply x by transform matrix */
	for (j = 0; j <= 7; ++j) {
	    t = x & tal[j];
	    t ^= (t >> 4);
	    t ^= (t >> 2);
	    t ^= (t >> 1);
	    t &= 0x01;
	    t <<= (7 - j);
	    r |= t;
	}
	convToDual[x] = r;
	dualToConv[r] = x;
    }

    /* create log and alog tables for GF(2^8) */
    for (x = 0x01, i = 0; i <= 254; ++i) {
	logAlpha[x] = i;
	alogAlpha[i] = x;
	if (x & 0x80)
	    x = (x << 1) ^ poly;
	else
	    x = (x << 1);
    }

    /* create multiplication table for GF(2^8) */
    for (i = 0; i <= 255; ++i) {
	for (j = 0; j <= 255; ++j) {
	    if (i == 0 || j == 0) {
		mulAlpha[i][j] = 0;
	    }
	    else {
		k = logAlpha[i] + logAlpha[j];
		if (k >= 255)
		    k -= 255;
		mulAlpha[i][j] = alogAlpha[k];
	    }
	}
    }

    /* store generator polynomial coefficients */
    g[0] = g[30] = alogAlpha[249];
    g[1] = g[29] = alogAlpha[59];
    g[2] = g[28] = alogAlpha[66];
    g[3] = g[27] = alogAlpha[4];
    g[4] = g[26] = alogAlpha[43];
    g[5] = g[25] = alogAlpha[126];
    g[6] = g[24] = alogAlpha[251];
    g[7] = g[23] = alogAlpha[97];
    g[8] = g[22] = alogAlpha[30];
    g[9] = g[21] = alogAlpha[3];
    g[10] = g[20] = alogAlpha[213];
    g[11] = g[19] = alogAlpha[50];
    g[12] = g[18] = alogAlpha[66];
    g[13] = g[17] = alogAlpha[170];
    g[14] = g[16] = alogAlpha[5];
    g[15] = alogAlpha[24];
    g[31] = alogAlpha[0];

    /* precompute products of feedback term and generator polynomial */
    for (i = 0; i < 256; ++i)
        for (j = 0; j < 32; ++j)
            products[i].u8[j] = mulAlpha[i][g[j]];
}















/* rsEncode -- Reed-Solomon encoder */

void
rsEncode(unsigned char* codeBlock, int interleave, int virtualFill)
{
    int encoder, i, j;
    unsigned char feedBack, *P;
    big_word checkSymbols[8];
    big_word *p;
    big_word *q;
    
    /* check parameters */
    assert(codeBlock != NULL);
    assert(interleave >= 1 && interleave <= 8);
    assert(virtualFill >= 0 && virtualFill <= 222);

    /* clear encoding register(s) */
    memset(checkSymbols, 0, 32*interleave);

    /* encode codeblock(s) */
    for (encoder = 0, i = 0; i < (K - virtualFill)*interleave; ++i) {
        /* convert to conventional basis and produce feedback term */
        feedBack = dualToConv[codeBlock[i]] ^ checkSymbols[encoder].u8[0];

        /* multiply by generator polynomial and add to register */
        p = &checkSymbols[encoder];
        q = &products[feedBack];

	/* the original code look like this (p and q were byte pointers:

        p[0] = p[1] ^ q[0];
        p[1] = p[2] ^ q[1];
        p[2] = p[3] ^ q[2];
        p[3] = p[4] ^ q[3];
        p[4] = p[5] ^ q[4];
        p[5] = p[6] ^ q[5];
        p[6] = p[7] ^ q[6];
        p[7] = p[8] ^ q[7];
        p[8] = p[9] ^ q[8];
        p[9] = p[10] ^ q[9];
        p[10] = p[11] ^ q[10];
        p[11] = p[12] ^ q[11];
        p[12] = p[13] ^ q[12];
        p[13] = p[14] ^ q[13];
        p[14] = p[15] ^ q[14];
        p[15] = p[16] ^ q[15];
        p[16] = p[17] ^ q[16];
        p[17] = p[18] ^ q[17];
        p[18] = p[19] ^ q[18];
        p[19] = p[20] ^ q[19];
        p[20] = p[21] ^ q[20];
        p[21] = p[22] ^ q[21];
        p[22] = p[23] ^ q[22];
        p[23] = p[24] ^ q[23];
        p[24] = p[25] ^ q[24];
        p[25] = p[26] ^ q[25];
        p[26] = p[27] ^ q[26];
        p[27] = p[28] ^ q[27];
        p[28] = p[29] ^ q[28];
        p[29] = p[30] ^ q[29];
        p[30] = p[31] ^ q[30];
        p[31] = q[31]; */
#ifdef __alpha__
	{
	    unsigned long p0, p1, p2, p3, q0, q1, q2, q3, r0, r1, r2, r3;

	    p0 = p->u64[0]; p1 = p->u64[1]; p2 = p->u64[2]; p3 = p->u64[3];
	    q0 = q->u64[0]; q1 = q->u64[1]; q2 = q->u64[2]; q3 = q->u64[3]; 

	    r0 = ((p1 << 56) | (p0 >> 8)) ^ q0;
	    r1 = ((p2 << 56) | (p1 >> 8)) ^ q1;
	    r2 = ((p3 << 56) | (p2 >> 8)) ^ q2;
	    r3 =               (p3 >> 8)  ^ q3;

	    p->u64[0] = r0; p->u64[1] = r1; p->u64[2] = r2; p->u64[3] = r3;
	}
#else
	{
	    unsigned int p0, p1, p2, p3, p4, p5, p6, p7;
	    unsigned int q0, q1, q2, q3, q4, q5, q6, q7;
	    unsigned int r0, r1, r2, r3, r4, r5, r6, r7;

	    p0 = p->u32[0]; p1 = p->u32[1]; p2 = p->u32[2]; p3 = p->u32[3];
	    p4 = p->u32[4]; p5 = p->u32[5]; p6 = p->u32[6]; p7 = p->u32[7];
	    q0 = q->u32[0]; q1 = q->u32[1]; q2 = q->u32[2]; q3 = q->u32[3]; 
	    q4 = q->u32[4]; q5 = q->u32[5]; q6 = q->u32[6]; q7 = q->u32[7]; 

	    r0 = ((p1 << 24) | (p0 >> 8)) ^ q0;
	    r1 = ((p2 << 24) | (p1 >> 8)) ^ q1;
	    r2 = ((p3 << 24) | (p2 >> 8)) ^ q2;
	    r3 = ((p4 << 24) | (p3 >> 8)) ^ q3;
	    r4 = ((p5 << 24) | (p4 >> 8)) ^ q4;
	    r5 = ((p6 << 24) | (p5 >> 8)) ^ q5;
	    r6 = ((p7 << 24) | (p6 >> 8)) ^ q6;
	    r7 =               (p7 >> 8)  ^ q7;

	    p->u32[0] = r0; p->u32[1] = r1; p->u32[2] = r2; p->u32[3] = r3;
	    p->u32[4] = r4; p->u32[5] = r5; p->u32[6] = r6; p->u32[7] = r7;
	}
#endif

        /* switch to next encoder register */
        if (++encoder >= interleave)
            encoder = 0;
    }

    /* store check symbols in code block, converting back to dual basis */
    P = &codeBlock[(K - virtualFill)*interleave];
    for (i = 0; i < 32; ++i)
        for (encoder = 0; encoder < interleave; ++encoder)
            *P++ = convToDual[checkSymbols[encoder].u8[i]];
}
-- 
David Mosberger-Tang; Dept. of CS; The University of Arizona; Tucson, AZ 85721
david.mosberger@acm.org            http://www.cs.arizona.edu/scout/suncam.html


