Meeting Pearls 3/Meeting_Pearls_III.iso/Pearls/bench/BIS/bis.c

C/C++ Source or Header | 1995-05-17 | 3KB | 231 lines

#include <stdio.h>
#include <time.h>
#include <math.h>

#define LUT_SIZE 10000
#define IMA_SIZE 256

#define IX(f,x,y) ((f)->data[(f)->width * (y) + (x)])

#define c_add(r,a,b)  (r).x = (a).x + (b).x, (r).y = (a).y + (b).y
#define c_sub(r,a,b)  (r).x = (a).x + (b).x, (r).y = (a).y + (b).y
#define c_mul(r,a,b)  (r).x = (a).x * (b).x - (a).y * (b).y, \
                      (r).y = (a).x * (b).y + (a).y * (b).x
#define c_cmul(r,a,b) (r).x = (a).x * (b).x + (a).y * (b).y, \
                      (r).y = (a).x * (b).y - (a).y * (b).x

typedef struct {
    double x;
    double y;
} Complex;

typedef struct {
    float x;
    float y;
} ShortComplex;

typedef struct {
    short x;
    short y;
} ShortVec;

typedef struct {
    ShortVec u;
    ShortVec v;
} ShortVecPair;

typedef struct {
    int width;
    int height;
    float *data;
    float vec[1];
} Ima;

typedef struct {
    int width;
    int height;
    Complex *data;
    Complex vec[1];
} Fos;

typedef struct {
    int length;
    ShortVecPair *list;
    ShortVecPair vec[1];
} Lut;

typedef struct {
    int length;
    Complex *list;
    Complex vec[1];
} Bis;

Lut *the_lut;
Fos *the_fos;
Bis *the_bis;

Lut *alloc_lut(n)
int n;
{
    Lut *lut;

    lut = (Lut *)malloc(sizeof(Lut) + (n-1) * sizeof(ShortVecPair));
    if (!lut) fatal("can't alloc_lut");

    lut->length = n;
    lut->list = &lut->vec[0];

    return lut;
}

Bis *alloc_bis(n)
int n;
{
    Bis *bis;

    bis = (Bis *)malloc(sizeof(Bis) + (n-1) * sizeof(Complex));
    if (!bis) fatal("can't alloc_bis");

    bis->length = n;
    bis->list = &bis->vec[0];

    return bis;
}

Fos *alloc_fos(n)
int n;
{
    Fos *fos;

    fos = (Fos *)malloc(sizeof(Fos) + (n*n-1) * sizeof(Complex));
    if (!fos) fatal("can't alloc_fos");

    fos->width = n;
    fos->height = n;
    fos->data = &fos->vec[0];

    return fos;
}

cleanup()
{
    if (the_bis) free(the_bis);
    if (the_fos) free(the_fos);
    if (the_lut) free(the_lut);
}

fatal(s)
char *s;
{
    if (s) fprintf(stderr,"FATAL: %s\n",s);
    cleanup();
    exit(99);
}

bispectrum(fos,lut,bis)
Fos *fos;
Lut *lut;
Bis *bis;
{
    int i;
    ShortVec u,v;
    Complex b,c;
    Complex one, two;

    for (i=0; i<lut->length; ++i) {
        u = lut->list[i].u;
        v = lut->list[i].v;
        one = IX(fos, u.x, u.y);
        two = IX(fos, v.x, v.y);
        c_mul(b, one, two);
        one = IX(fos, u.x+v.x, u.y+v.y);
        c_cmul(c, b, one);
        b = bis->list[i];
        c_add(b, b, c);
        bis->list[i] = b;
    }
}

clear_bis(bis)
Bis *bis;
{
    int i;

    for (i=0; i<bis->length; ++i) {
        bis->list[i].x = 0.0;
        bis->list[i].y = 0.0;
    }
}

fake_fos(fos)
Fos *fos;
{
    int i,j;

    for (i=0; i<IMA_SIZE; ++i)
        for (j=0; j<IMA_SIZE; ++j) {
            IX(fos,j,i).x = exp(-(double)(i*i+j*j));
            IX(fos,j,i).y = 1.0;
        }
}

void fake_lut(lut)
Lut *lut;
{
    int i,j;
    int n;
    ShortVec u,v;

    n = 0;
    for (i=0; i<IMA_SIZE; ++i)
        for (j=0; j<=i; ++j) {
            u.x = i;
            u.y = j;
            v.x = IMA_SIZE-1 - i;
            v.y = IMA_SIZE-1 - j;
            lut->list[n].u = u;
            lut->list[n].v = v;
            ++n;
            if (n >= lut->length) return;
        }
}

double sumall(bis)
Bis *bis;
{
    int i;
    double sum;

    sum = 0.0;
    for (i=0; i<bis->length; ++i)
        sum += bis->list[i].x + bis->list[i].y;
    return sum;
}

main()
{
    time_t start,stop;
    int i;

    the_lut = alloc_lut(LUT_SIZE);
    the_fos = alloc_fos(IMA_SIZE);
    the_bis = alloc_bis(LUT_SIZE);

    fake_fos(the_fos);
    fake_lut(the_lut);
    clear_bis(the_bis);

    time(&start);

    for (i=0; i<1000; ++i)
        bispectrum(the_fos,the_lut,the_bis);

    time(&stop);

    printf("sum = %f\n",sumall(the_bis));

    printf("%d seconds\n",stop-start);

    cleanup();
}