午後なコード

float func1(float a, float *buf, int n){
    int i;
    float x = 0;
    for(i = 0; i < n; i++){
        x += buf[i] / a;
    }
    return x;
}

float func1a(float a, float *buf, int n){
    int i;
    float x = 0;
    for(i = 0; i < n; i++){
        x += buf[i];
    }
    return x / a;
}

float func2(float a, float *buf, int n){
    int i;
    float x = 1;
    for(i = 0; i < n; i++){
        x *= pow(buf[i], a);
    }
    return x;
}

float func2a(float a, float *buf, int n){
    int i;
    float x = 1;
    for(i = 0; i < n; i++){
        x *= buf[i];
    }
    return pow(x, a);
}

float func2b(float a, float *buf, int n){
    int i;
    float x = 1, y;
    int r, q;
    q = n >> 3; /* 8 で割った商 */
    r = n & 7;  /* 8 で割った余り */
    for(i = 0; i < q; i++){
        int j;
        y = 1;
        for(j = 0; j < 8; j++) y *= buf[i];
        x *= pow(y, a);
    }
    if(r){
        y = 1;
        for(i = 0; i < r; i++) y *= buf[i];
        x *= pow(y, a);
    }
    return x;
}

float func2c(float a, float *buf, int n){
    int i;
    float x = 0;
    for(i = 0; i < n; i++){
        x += log(buf[i]);
    }
    return exp(a * x);
}

#define ADJ (1.0e-10)
float func2d(float a, float *buf, int n){
    int i;
    double x, y, z;
    x = 1;
    for(i = 0; i < n; i++){
        x *= (buf[i] * ADJ);
    }
#if 0
    /* もしn, aが固定ならz(=pow(ADJ, -n*a)は先に計算しておける */
#else
    y = pow(ADJ, -a);
    z = 1;
    while(n){
        if(n & 1) z *= y;
        y *= y;
        n >>= 1;
    }
#endif
    return pow(x, a) * z;
}

/* func2cのlogを分解 */
float func2ca(float a, float *buf, int n){
    int i;
    int sum_s;
    float mul_t, x;
    sum_s = 0;
    mul_t = 1;
    for(i = 0; i < n; i++){
        int s, *p;
        volatile float t = buf[i];     /* buf[i]を破壊するのでコピー */
        p = (int *)&t;
        s = (*p >> 23) - 127;             /* 指数部 */
        *p = *p & ((1<<23)-1) | (127<<23);/* 仮数部 */
        sum_s += s;
        mul_t *= t;
    }
    x = sum_s * log(2) + log(mul_t); /* 底の変換 */
    return exp(a * x);
}

#define BITLEN 10
#define SIZE (1<<BITLEN)
const float a = 0.1;
static float idxTbl[256], powTbl[SIZE];

/* 第0次近似によるpow(x, a)の計算 */
float pow_fast(float x){
    static int initFlag = 1;

    float t;
    int s, *p;
    if(initFlag){
        int i;
        for(i = 0; i < 256; i++){
            idxTbl[i] = pow(2, (i - 127) * a);
        }
        for(i = 0; i < SIZE; i++){
            powTbl[i] = pow(1 + i / (double)SIZE, a);
        }
        initFlag = 0;
    }
    t = x;
    p = (int *)&t;
    s = (*p >> 23);                   /* 指数部 */
    *p = *p & ((1<<23)-1) | (127<<23);/* 仮数部 */
    return idxTbl[s] * powTbl[(int)((t - 1) * SIZE)];
}

    powTbl[(int)((t - 1) * SIZE)];

    powTbl[(int)(t * SIZE) - SIZE];

    *p = *p & ((1<<23)-1) | ((127+BITLEN)<<23);/* 仮数部 */
    return idxTbl[s] * powTbl[(int)(t) - SIZE)];

__ftol:
    push        ebp
    mov         ebp,esp
    add         esp,0FFFFFFF4h
    wait
    fnstcw      word ptr [ebp-2]
    wait
    mov         ax,word ptr [ebp-2]
    or          ah,0Ch
    mov         word ptr [ebp-4],ax
    fldcw       word ptr [ebp-4]
    fistp       qword ptr [ebp-0Ch]
    fldcw       word ptr [ebp-2]
    mov         eax,dword ptr [ebp-0Ch]
    mov         edx,dword ptr [ebp-8]
    leave
    ret

_f05:    0.5

    fld     dword [x]
    fsub    dword [_f05]     ;0.5を引いて
    fistp   dword [int_x]    ;四捨五入=>切捨て

;int loadCW(void);
proc    loadCW
        push    eax
        fstcw   word [esp]
        pop     eax
        ret
    
;void changeCW(int cw);
proc    changeCW
        mov     eax,[esp+4]
        push    eax
        fldcw   word [esp]
        pop     eax
        ret

    int cw;
    cw = loadCW();
    changeCW(cw | (3<<10)); /* 0方向へ切り捨て */

    /* intへのcastにfistを使った演算ループ */

    changeCW(cw);   /* 元のモードに戻す */

/* pow_fastの改良(テーブル初期化部分は除く) */
float pow_fast2(float x){
    int s, t;
    s = t = *(int *)&x;
    s >>= 23;                               /* 指数部 */
    t = (t & ((1<<23)-1)) >> (23 - BITLEN); /* 仮数部 */
    return idxTbl[s] * powTbl[t];
}

/* 一次補間ありの平方根 */
#define BITLEN 8
#define SIZE (1<<BITLEN)
static float idxTbl[256];
static float sqrtTbl1[SIZE];
static float sqrtTbl2[SIZE];

init_sqrt_fast(void){
	int i;
	for(i = 0; i < 256; i++){
		idxTbl[i] = (float)pow(2, (i - 127) * 0.5);
	}
	for(i = 0; i < SIZE; i++){
		sqrtTbl1[i] = (float)sqrt(1 + i / (double)SIZE);
	}
	for(i = 0; i < SIZE-1; i++){
		sqrtTbl2[i] = (sqrtTbl1[i+1] - sqrtTbl1[i]) / (float)SIZE;
	}
	sqrtTbl2[SIZE-1] = ((float)sqrt(2) - sqrtTbl1[SIZE-1]) / (float)SIZE;
}

static float sqrt_fast(float x){
#if 0	/* 補完無し */
	int s, t;
	s = t = *(int *)&x;
	s >>= 23;
	t &= ((1<<23)-1);
	t >>= (23-BITLEN);
	return idxTbl[s] * sqrtTbl1[t];
#else	/* 一次補完 */
	int s, t, u;
	s = t = u = *(int *)&x;
	s >>= 23;
	t &= ((1<<23)-1);
	u &= ((1<<(23-BITLEN))-1);
	t >>= (23-BITLEN);
	u >>= (23-BITLEN*2);
	return idxTbl[s] * (sqrtTbl1[t] + sqrtTbl2[t] * u);
#endif
}

static float sqrtTbl[SIZE*2]; /* Tbl1とTbl2をインタリーブ */
init_sqrt_fast(void){
	int i;
	for(i = 0; i < 256; i++){
		idxTbl[i] = (float)pow(2, (i - 127) * 0.5);
	}
	for(i = 0; i < SIZE; i++){
		double a, b, sa, sb, t;
		a = 1 + i / (double)SIZE;
		b = a + 1 / (double)SIZE;
		sa = sqrt(a);
		sb = sqrt(b);
		t = (sa + sb) * (sa + sb);
		sqrtTbl[i*2] = (float)((b*sb + 2 * b*sa + 8 * a*sb + 4 * a*sa) * 4 / (15 * t));
		t *= (sa + sb);
		sqrtTbl[i*2+1] = (float)((a + 3 * sa*sb + b) * 4 / (5 * t * SIZE * SIZE));
	}
}

double log_fast(double x, double a, double b, double c){
    double z, zz;
    z = (x - sqrt(2)) / (x + sqrt(2));
    zz = z*z;
    return 0.5 * log(2) + 2 *z * (a + zz * (b+c*zz));
}

void calc_diff(double a, double b, double c){
    double diff, x, s, t;
    const int N = 1000;
    int i;
    diff = 0;
    for(i = 0; i < N; i++){
        x = i / (double)(N+1) + 1;
        s = log(x);
        t = log_fast(x, a, b, c);
        diff += fabs(s - t);
    }
    printf("coef a=1%+e, b=1/3%+e, c=1/5%+e\n", a-1, b-1.0/3, c-1.0/5);
    printf("ave. diff=%e\n", diff / N);
}

      (1/3 1/5 1/7)
    A=(1/5 1/7 1/9)
      (1/7 1/9 1/11)

float logC( float x )
{
    float a, z, zz;
    unsigned int *p;
    int b;
    const float A = 1.0 + 3.104605e-7;      /* A*2 = 2.0000006209 */
    const float B = 1.0/3 - 9.440747e-5;    /* B*2 = 0.6664778517 */
    const float C = 1.0/5 + 6.987293e-3;    /* C*2 = 0.4139745860 */
    p = (unsigned int *)&x;
    b = (*p >> 23) - 127;   /* suppose x_org > 0 */
    *p= ( *p & 0x007FFFFF ) | 0x3F800000;
    a = *(float *)p;
    z = (a - SQRT2) / (a + SQRT2); 
    zz = z * z;
    z = (b+0.5) * LOG_E_2 + z * ((A*2) + zz * ((B*2) + zz * (C*2)));
    return z;
}

心理解析部分(){
    foreach(チャンネル[右, 左, まんなか, 両端]){
        ごちゃごちゃっと300〜400行;
    }
}

#define SIZE 128
static void test1(const int src[SIZE][SIZE], int dst[SIZE][SIZE])
{
    int i, j;
    for (i = 0; i < SIZE; i++) {
        for (j = 0; j < SIZE; j++) {
            dst[j][i] = abs(src[j][i]);
        }
    }
}

static void test2(const int src[SIZE][SIZE], int dst[SIZE][SIZE])
{
    int i, j;
    for (i = 0; i < SIZE; i++) {
        for (j = 0; j < SIZE; j++) {
            dst[i][j] = abs(src[i][j]);
        }
    }
}

int main(void)
{
    int src[SIZE][SIZE];
    int dst[SIZE][SIZE];
    test1(src, dst);
    test2(src, dst);
    return 0;
}

typedef unsigned char uint8;
typedef signed char sint8;
void add_us_s_and_clip(uint8 *p_us, sint8 *p_s, int count)
{
    int i;
    for (i = 0; i < count; i++) {
        int a, b;
        a = p_us[i];
        b = p_s[i];
        a += b;
        if (a < 0) {
            a = 0;
        } else if (a > 255) {
            a = 255;
        }
        p_us[i] = (uint8)a;
    }
}

• b ≧ 0 かつ a + b ≦ 255
• b ≧ 0 かつ a + b ＞ 255
• b ＜ 0 かつ a + b ≧ 0
• b ＜ 0 かつ a + b ＜ 0

1. 符号無し同士の飽和加算

2. 符号付き同士の飽和加算

3. 符号無しと符号付き同士の飽和加算

    add_us_s_and_clip(a, b) = paddsb(a^0x80, b) ^ 0x80

%rep 20
%if 0
	pxor		mm0, mm0
	paddq		mm0, [esi]
	movq		mm1, mm0
	pand		mm1, mm7
	movq		mm1, mm0		; [H:L]
	psrad		mm1, 31		; [sign H:*]
	punpckhdq	mm0, mm1		; [sign H:H]
%else
	movq		mm0, [esi]
	movq		mm1, mm0
	pand		mm1, mm7
	movq		mm1, mm0		; [H:L]
	psrad		mm1, 31		; [sign H:*]
	punpckhdq	mm0, mm1		; [sign H:H]
%endif
%endrep

[Back]
御意見は光成滋生<herumi@nifty.com>までお願いします