Hlavní navigace

Algoritmus CORDIC s hodnotami uloženými ve formátu FX

Pavel Tišnovský 26. 7. 2006

V předposlední části seriálu si ukážeme algoritmy realizované v programovacím jazyku C, ve kterých budou implementovány výpočty některých goniometrických funkcí pomocí algoritmu CORDIC. Veškeré výpočty přitom budou prováděny pouze s hodnotami uloženými ve formátu pevné řádové binární čárky (FX).

Obsah

1. Algoritmus CORDIC s hodnotami uloženými ve formátu pevné řádové binární čárky
2. Datové typy, makra a pomocné funkce pro práci s formátem FX
3. Funkce určená pro naplnění tabulky atans[]
4. Výpočet tangenty pomocí algoritmu CORDIC – neoptimalizovaná verze
5. Výpočet tangenty pomocí algoritmu CORDIC – optimalizovaná verze
6. Výpis celého demonstračního příkladu
7. Obsah poslední části tohoto seriálu

1. Algoritmus CORDIC s hodnotami uloženými ve formátu pevné řádové binární čárky

V předchozích částech tohoto seriálu jsme si uvedli, jakým způsobem je možné algoritmus CORDIC (COordinate ROtation DIgital Computer) použít pro výpočet goniometrických funkcí s hodnotami uloženými ve formátu plovoucí řádové binární čárky (FP). Tentýž algoritmus je však po mírné modifikaci možné použít i při práci s formátem pevné řádové binární čárky (FX) a dá se říci, že teprve zde se plně ukazuje jeho implementační jednoduchost a současně i velká vyjadřovací síla – pomocí jednoho iteračního postupu je možné vypočítat poměrně velké množství navzájem značně odlišných funkcí, přičemž paměťové nároky algoritmu jsou minimální a rovnají se tabulce s cca deseti až dvaceti prvky s rozsahem 32 (popř. pouze 16) bitů. To je výhodné především při implementaci na mikrořadičích (oblíbená řada 8051, PICy, řada 68HC11 atd.) a také při práci s programovatelnými obvody typu FPGA.

2. Datové typy, makra a pomocné funkce pro práci s formátem FX

Před vlastním uvedením implementace algoritmu CORDIC s čísly uloženými ve formátu pevné řádové binární čárky je zapotřebí provést náležitou přípravu, protože formát FX (většinou) není ani mikroprocesory ani překladači podporován. Z tohoto důvodu si v této kapitole pro připomenutí uvedeme všechny základní aritmetické funkce určené pro práci s formátem FX. Jedná se o funkce fx_add(), fx_sub(), fx_mul() a fx_div() (tyto funkce byly podrobněji vysvětleny v předchozím díle). Kromě toho jsou uvedena i těla dalších pomocných funkcí, zejména fx_print() i maker určených pro převod stupňů na radiány a naopak. Pomocí konstant A a B je určen rozsah a přesnost numerických hodnot uložených ve formátu fx – ten je v naší implementaci shodný s třicetidvou­bitovým datovým typem signed int, jelikož potřebujeme pracovat jak s kladnými, tak i se zápornými hodnotami.

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

/* počet míst před a za binární řádovou tečkou */
#define A 16
#define B 16

/* Ludolfovo číslo */
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* maximální počet iterací při běhu algoritmu */
#define MAX_ITER 16

/* "zesílení" při rotacích */
#define K 0.6073

/* převody mezi stupni a radiány */
#define rad2deg(rad) ((rad)*180.0/M_PI)
#define deg2rad(deg) ((deg)/180.0*M_PI)

/* datový typ, se kterým budeme pracovat */
typedef signed int fx;

/* hlavičky použitých funkcí */
void   fx_print(fx x);
fx     fp2fx(double x);
double fx2fp(fx x);

/* tabulka arkustangentu úhlů */
fx atans[MAX_ITER];

/* tabulka záporných celočíselných mocnin hodnoty 2 */
fx pows[MAX_ITER];

/*
 * Tisk numerické hodnoty uložené ve formátu pevné
 * řádové binární čárky (FX)
 */
void fx_print(fx x)
{
    int i;
    int val=x;                         /* pomocná proměnná pro převod do dvojkové soustavy */
    printf("bin: ");
    for (i=0; i<A+B; i++) {            /* převod na řetězec bitů (do dvojkové soustavy) */
        putchar(!!(val & (1<<(A+B-1)))+'0'); /* výpis hodnoty aktuálně nejvyššího bitu */
        if (i==B-1) putchar('.');      /* po řádové binární čárce vypsat značku */
        val=val<<1;                    /* posun na další (méně významný) bit */
    }

    printf("   hex: %08x   fp: %+11.5f\n", x, fx2fp(x));
}

/*
 * Převod z formátu plovoucí řádové binární čárky (FP)
 * do formátu pevné řádové binární čárky (FX)
 */
fx fp2fx(double x)
{
    return (fx)(x*(2<<(B-1)));
}

/*
 * Převod z celočíselného formátu (integer)
 * do formátu pevné řádové binární čárky (FX)
 */
fx int2fx(int x)
{
    return (fx)(x<<B);
}

/*
 * Převod z formátu pevné řádové binární čárky (FX)
 * do formátu plovoucí řádové binární čárky (FP)
 */
double fx2fp(fx x)
{
    return (double)x/(2<<(B-1));
}

/*
 * Součet dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_add(fx x, fx y)
{
    return x+y;
}

/*
 * Rozdíl dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_sub(fx x, fx y)
{
    return x-y;
}

/*
 * Součin dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_mul(fx x, fx y)
{
    fx result=(x>>(B/2))*(y>>(B/2));
    return result;
}

/*
 * Podíl dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_div(fx x, fx y)
{
    fx result=x/(y>>(B/2));
    return result<<(B/2);
} 

3. Funkce určená pro naplnění tabulky atans[]

Podobně jako při implementaci algoritmu CORDIC pro hodnoty ve formátu FP, i v případě práce s formátem FX je nutné některé hodnoty předpočítat. Opět se jedná o tabulku atans[] obsahující arkustangenty vybraných úhlů. Kromě toho naplníme i tabulku pows[], ve které jsou uloženy záporné mocniny čísla 2. Jak však uvidíme v dalším textu, není v případě FX výpočtů tato tabulka využita a algoritmus CORDIC se poněkud zjednoduší. Funkce pro výpočet obou zmiňovaných tabulek vypadá následovně:

/*
 * Vytvoření tabulky pro výpočet goniometrických
 * funkcí pomocí algoritmu CORDIC
 */
void fx_create_tables(void)
{
    int i;
    for (i=0; i<MAX_ITER; i++) {
        double p=pow(2.0, -i);
        atans[i]=fp2fx(atan(p));
        pows[i]=fp2fx(p);
    }
} 

To, zda je tabulka korektně naplněna, je možné otestovat velmi jednoduše:

    /* kontrolní výpis tabulky atans[] */
    for (i=0; i<MAX_ITER; i++)
        printf("%d\t%f\n", i, fx2fp(rad2deg(atans[i]))); 

Výsledkem běhu předchozího testu je tabulka hodnot arkustangent úhlů (zhruba) odpovídajících FP verzi:

Index Hodnota úhlu
0 44.999252
1 26.564514
2 14.035431
3 7.124374
4 3.575729
5 1.789612
6 0.894363
7 0.446747
8 0.222931
9 0.111023
10 0.055069
11 0.027100
12 0.013107
13 0.006119
14 0.002609
15 0.000870

S takto připravenou tabulkou je možné pokračovat v dalších výpočtech.

4. Výpočet tangenty pomocí algoritmu CORDIC – neoptimalizova­ná verze

Pravděpodobně nejjednodušším výpočtem, který je možné pomocí algoritmu CORDIC provést, je vyjádření tangenty zadaného úhlu. Nejprve si uvedeme neoptimalizovanou verzi, která vznikla v podstatě přímým převodem dříve uvedené implementace určené pro FP reprezentaci. Převod spočívá v náhradě aritmetických funkcí jejich FX ekvivalenty:

/* výpočet funkce tan() pro zadaný úhel delta */
/* (neoptimalizovaná verze) */
fx fx_tan_cordic(fx delta)
{
    int i;
    /* nastavení počátečních podmínek */
    fx x0=fp2fx(1.0);
    fx y0=fp2fx(0.0);
    fx xn;
    if (delta==0) return 0;                 /* ošetření nulového úhlu */
    for (i=0; i<MAX_ITER; i++) {            /* iterační smyčka */
        if (delta<0) {                      /* úhel je záporný => rotace doleva */
            xn=fx_add(x0, fx_mul(y0, pows[i]));
            y0=fx_sub(y0, fx_mul(x0, pows[i]));
            delta=fx_add(delta, atans[i]);
        }
        else {                              /* úhel je kladný => rotace doprava */
            xn=fx_sub(x0, fx_mul(y0, pows[i]));
            y0=fx_add(y0, fx_mul(x0, pows[i]));
            delta=fx_sub(delta, atans[i]);
        }
        x0=xn;
/*        printf("%i\t%+f\t%+f\t%+f\n", i, fx2fp(x0), fx2fp(y0), fx2fp(delta)); */
    }
    if (x0==0)                              /* ošetření tangenty pravého úhlu */
        if (y0<0) return 0;
        else return 0;
    else
        return fx_div(y0,x0);               /* vrátit výsledek operace */
} 

Tuto funkci si můžeme jednoduchým způsobem otestovat, postačí použít smyčku, ve které se počítají tangenty úhlů v rozsahu 0°..89° spolu s vyjádřením absolutních a relativních chyb vzniklých použitím algoritmu CORDIC a bitově omezené FP reprezentace:

    /* výpis tabulky tangent úhlů v rozsahu 0..89° */
    for (i=0; i<90; i++) {                  /* výpočetní smyčka */
        delta=deg2rad(i);                   /* převod úhlu na radiány */
        tanfx=fx_tan_cordic(fp2fx(delta));  /* aplikace algoritmu CORDIC */
        tanval=fx2fp(tanfx);                /* výpočet funkce tan */
        tanerr=fabs(tanval-tan(delta));     /* výpočet absolutních chyb */
                                            /* tisk výsledků */
        printf("%02d\t%14.10f\t%14.10f\t%12.10f\t%8.3f%%\n",
                i,
                tanval,
                tan(delta),
                tanerr,
                tanerr==0 ? 0:100.0*tanerr/tan(delta));
    } 

V tabulce s výsledky si všimněte, že pro některé úhly je vzniklá chyba rovna až deseti procentům. Je to způsobeno nízkým rozlišením použitého FX formátu, zejména hodnot arctat() těch nejmenších úhlů, a chybu již není možné snížit ani zvýšením počtu iterací (ten je nastaven na šestnáct, to odpovídá počtu bitů za binární řádovou čárkou):

Úhel tan podle CORDIC tan podle FPU absolutní chyba relativní chyba
00 0.0000000000 0.0000000000 0.0000000000 0.000%
01 0.0156250000 0.0174550649 0.0018300649 10.484%
02 0.0312500000 0.0349207695 0.0036707695 10.512%
03 0.0468750000 0.0524077793 0.0055327793 10.557%
04 0.0703125000 0.0699268119 0.0003856881 0.552%
05 0.0859375000 0.0874886635 0.0015511635 1.773%
06 0.1015625000 0.1051042353 0.0035417353 3.370%
07 0.1171875000 0.1227845609 0.0055970609 4.558%
08 0.1406250000 0.1405408347 0.0000841653 0.060%
09 0.1562500000 0.1583844403 0.0021344403 1.348%
10 0.1718750000 0.1763269807 0.0044519807 2.525%
11 0.1914062500 0.1943803091 0.0029740591 1.530%
12 0.2148437500 0.2125565617 0.0022871883 1.076%
13 0.2304687500 0.2308681911 0.0003994411 0.173%
14 0.2460937500 0.2493280028 0.0032342528 1.297%
15 0.2656250000 0.2679491924 0.0023241924 0.867%
16 0.2890625000 0.2867453858 0.0023171142 0.808%
17 0.3046875000 0.3057306815 0.0010431815 0.341%
18 0.3242187500 0.3249196962 0.0007009462 0.216%
19 0.3398437500 0.3443276133 0.0044838633 1.302%
20 0.3671875000 0.3639702343 0.0032172657 0.884%
21 0.3828125000 0.3838640350 0.0010515350 0.274%
22 0.4023437500 0.4040262258 0.0016824758 0.416%
23 0.4218750000 0.4244748162 0.0025998162 0.612%
24 0.4414062500 0.4452286853 0.0038224353 0.859%
25 0.4687500000 0.4663076582 0.0024423418 0.524%
26 0.4882812500 0.4877325886 0.0005486614 0.112%
27 0.5078125000 0.5095254495 0.0017129495 0.336%
28 0.5273437500 0.5317094317 0.0043656817 0.821%
29 0.5546875000 0.5543090515 0.0003784485 0.068%
30 0.5781250000 0.5773502692 0.0007747308 0.134%
31 0.5976562500 0.6008606190 0.0032043690 0.533%
32 0.6171875000 0.6248693519 0.0076818519 1.229%
33 0.6523437500 0.6494075932 0.0029361568 0.452%
34 0.6718750000 0.6745085168 0.0026335168 0.390%
35 0.6953125000 0.7002075382 0.0048950382 0.699%
36 0.7226562500 0.7265425280 0.0038862780 0.535%
37 0.7539062500 0.7535540501 0.0003521999 0.047%
38 0.7812500000 0.7812856265 0.0000356265 0.005%
39 0.8085937500 0.8097840332 0.0011902832 0.147%
40 0.8320312500 0.8390996312 0.0070683812 0.842%
41 0.8750000000 0.8692867378 0.0057132622 0.657%
42 0.8984375000 0.9004040443 0.0019665443 0.218%
43 0.9296875000 0.9325150861 0.0028275861 0.303%
44 0.9570312500 0.9656887748 0.0086575248 0.897%
45 0.9921875000 1.0000000000 0.0078125000 0.781%
46 1.0390625000 1.0355303138 0.0035321862 0.341%
47 1.0742187500 1.0723687100 0.0018500400 0.173%
48 1.1093750000 1.1106125148 0.0012375148 0.111%
49 1.1445312500 1.1503684072 0.0058371572 0.507%
50 1.2031250000 1.1917535926 0.0113714074 0.954%
51 1.2421875000 1.2348971565 0.0072903435 0.590%
52 1.2812500000 1.2799416322 0.0013083678 0.102%
53 1.3203125000 1.3270448216 0.0067323216 0.507%
54 1.3867187500 1.3763819205 0.0103368295 0.751%
55 1.4335937500 1.4281480067 0.0054457433 0.381%
56 1.4843750000 1.4825609685 0.0018140315 0.122%
57 1.5351562500 1.5398649638 0.0047087138 0.306%
58 1.6132812500 1.6003345290 0.0129467210 0.809%
59 1.6718750000 1.6642794824 0.0075955176 0.456%
60 1.7265625000 1.7320508076 0.0054883076 0.317%
61 1.7929687500 1.8040477553 0.0110790053 0.614%
62 1.9023437500 1.8807264653 0.0216172847 1.149%
63 1.9765625000 1.9626105055 0.0139519945 0.711%
64 2.0585937500 2.0503038416 0.0082899084 0.404%
65 2.1406250000 2.1445069205 0.0038819205 0.181%
66 2.2656250000 2.2460367739 0.0195882261 0.872%
67 2.3789062500 2.3558523658 0.0230538842 0.979%
68 2.4882812500 2.4750868534 0.0131943966 0.533%
69 2.6015625000 2.6050890647 0.0035265647 0.135%
70 2.7265625000 2.7474774195 0.0209149195 0.761%
71 2.9296875000 2.9042108777 0.0254766223 0.877%
72 3.0820312500 3.0776835372 0.0043477128 0.141%
73 3.2460937500 3.2708526185 0.0247588685 0.757%
74 3.4570312500 3.4874144438 0.0303831938 0.871%
75 3.7695312500 3.7320508076 0.0374804424 1.004%
76 4.0468750000 4.0107809335 0.0360940665 0.900%
77 4.3203125000 4.3314758743 0.0111633743 0.258%
78 4.6289062500 4.7046301095 0.0757238595 1.610%
79 5.2382812500 5.1445540160 0.0937272340 1.822%
80 5.7656250000 5.6712818196 0.0943431804 1.664%
81 6.3046875000 6.3137515147 0.0090640147 0.144%
82 7.0703125000 7.1153697224 0.0450572224 0.633%
83 8.3710937500 8.1443464280 0.2267473220 2.784%
84 9.7500000000 9.5143644542 0.2356355458 2.477%
85 11.3476562500 11.4300523028 0.0823960528 0.721%
86 14.0117187500 14.3006662567 0.2889475067 2.021%
87 20.0468750000 19.0811366877 0.9657383123 5.061%
88 30.0937500000 28.6362532829 1.4574967171 5.090%
89 60.2109375000 57.2899616308 2.9209758692 5.099%

5. Výpočet tangenty pomocí algoritmu CORDIC – optimalizovaná ver­ze

Pokud se nad výše uvedenou funkcí fx_tan_cordic() zamyslíme, pravděpodobně zjistíme, že se v ní zbytečně provádí některé aritmetické operace. Zejména se jedná o průběžné násobení koeficienty uloženými v tabulce pows[]. Tato tabulka měla svůj význam při práci s FP formátem, u formátu FX však ztrácí prakticky všechny své výhody, protože násobení zápornou mocninou čísla 2 je vlastně totožné s bitovým posunem doprava. Optimalizovaná funkce pro výpočet tangenty úhlu algoritmem CORDIC vypadá následovně:

/* výpočet funkce tan() pro zadaný úhel delta */
/* (optimalizovaná verze) */
fx fx_tan_cordic_optim(fx delta)
{
    int i;
    /* nastavení počátečních podmínek */
    fx x0=int2fx(1);
    fx y0=0;
    fx xn;
    if (delta==0) return 0;                 /* ošetření nulového úhlu */
    for (i=0; i<MAX_ITER; i++) {            /* iterační smyčka */
        if (delta<0) {                      /* úhel je záporný => rotace doleva */
            xn=fx_add(x0, y0>>i);           /* místo násobení bitový posuv */
            y0=fx_sub(y0, x0>>i);
            delta=fx_add(delta, atans[i]);
        }
        else {                              /* úhel je kladný => rotace doprava */
            xn=fx_sub(x0, y0>>i);
            y0=fx_add(y0, x0>>i);
            delta=fx_sub(delta, atans[i]);
        }
        x0=xn;
    }
    if (x0==0)                              /* ošetření tangenty pravého úhlu */
        if (y0<0) return 0;
        else return 0;
    else
        return fx_div(y0,x0);               /* vrátit výsledek operace */
} 

(Pokud by někdo pocítil potřebu provádět další optimalizace, nabízí se zde použití MMX instrukcí s paralelním výpočtem obou alternativních větví). Za povšimnutí stojí fakt, že se v iterační smyčce používají pouze aritmetické operace součtu, rozdílu a bitového posuvu. Kromě toho se i při inicializaci hodnot x0 a y0 používá pouze celočíselná aritmetika, tj. ve výpočtech se vůbec nevyskytuje volání instrukcí matematického koprocesoru. Tato absence složitých a implementačně náročných operací představuje další z důvodů velké oblíbenosti algoritmu CORDIC v komunitě vývojářů pro malé mikrořadiče a mikroprocesory.

Pokud výše uvedenou funkci fx_tan_cordic_op­tim() podrobíme stejnému testu, jako funkci fx_tan_cordic(), dostaneme následující (do značné míry podobné) výsledky:

Úhel tan podle CORDIC tan podle FPU absolutní chyba relativní chyba
00 0.0000000000 0.0000000000 0.0000000000 0.000%
01 0.0156250000 0.0174550649 0.0018300649 10.484%
02 0.0312500000 0.0349207695 0.0036707695 10.512%
03 0.0507812500 0.0524077793 0.0016265293 3.104%
04 0.0664062500 0.0699268119 0.0035205619 5.035%
05 0.0859375000 0.0874886635 0.0015511635 1.773%
06 0.1015625000 0.1051042353 0.0035417353 3.370%
07 0.1210937500 0.1227845609 0.0016908109 1.377%
08 0.1406250000 0.1405408347 0.0000841653 0.060%
09 0.1562500000 0.1583844403 0.0021344403 1.348%
10 0.1757812500 0.1763269807 0.0005457307 0.309%
11 0.1914062500 0.1943803091 0.0029740591 1.530%
12 0.2109375000 0.2125565617 0.0016190617 0.762%
13 0.2304687500 0.2308681911 0.0003994411 0.173%
14 0.2460937500 0.2493280028 0.0032342528 1.297%
15 0.2656250000 0.2679491924 0.0023241924 0.867%
16 0.2851562500 0.2867453858 0.0015891358 0.554%
17 0.3046875000 0.3057306815 0.0010431815 0.341%
18 0.3242187500 0.3249196962 0.0007009462 0.216%
19 0.3437500000 0.3443276133 0.0005776133 0.168%
20 0.3632812500 0.3639702343 0.0006889843 0.189%
21 0.3828125000 0.3838640350 0.0010515350 0.274%
22 0.4023437500 0.4040262258 0.0016824758 0.416%
23 0.4218750000 0.4244748162 0.0025998162 0.612%
24 0.4453125000 0.4452286853 0.0000838147 0.019%
25 0.4648437500 0.4663076582 0.0014639082 0.314%
26 0.4882812500 0.4877325886 0.0005486614 0.112%
27 0.5078125000 0.5095254495 0.0017129495 0.336%
28 0.5312500000 0.5317094317 0.0004594317 0.086%
29 0.5546875000 0.5543090515 0.0003784485 0.068%
30 0.5742187500 0.5773502692 0.0031315192 0.542%
31 0.5976562500 0.6008606190 0.0032043690 0.533%
32 0.6250000000 0.6248693519 0.0001306481 0.021%
33 0.6484375000 0.6494075932 0.0009700932 0.149%
34 0.6718750000 0.6745085168 0.0026335168 0.390%
35 0.6992187500 0.7002075382 0.0009887882 0.141%
36 0.7265625000 0.7265425280 0.0000199720 0.003%
37 0.7539062500 0.7535540501 0.0003521999 0.047%
38 0.7812500000 0.7812856265 0.0000356265 0.005%
39 0.8085937500 0.8097840332 0.0011902832 0.147%
40 0.8398437500 0.8390996312 0.0007441188 0.089%
41 0.8671875000 0.8692867378 0.0020992378 0.241%
42 0.8984375000 0.9004040443 0.0019665443 0.218%
43 0.9296875000 0.9325150861 0.0028275861 0.303%
44 0.9648437500 0.9656887748 0.0008450248 0.088%
45 1.0000000000 1.0000000000 0.0000000000 0.000%
46 1.0351562500 1.0355303138 0.0003740638 0.036%
47 1.0742187500 1.0723687100 0.0018500400 0.173%
48 1.1093750000 1.1106125148 0.0012375148 0.111%
49 1.1523437500 1.1503684072 0.0019753428 0.172%
50 1.1953125000 1.1917535926 0.0035589074 0.299%
51 1.2343750000 1.2348971565 0.0005221565 0.042%
52 1.2812500000 1.2799416322 0.0013083678 0.102%
53 1.3281250000 1.3270448216 0.0010801784 0.081%
54 1.3789062500 1.3763819205 0.0025243295 0.183%
55 1.4296875000 1.4281480067 0.0015394933 0.108%
56 1.4843750000 1.4825609685 0.0018140315 0.122%
57 1.5429687500 1.5398649638 0.0031037862 0.202%
58 1.6015625000 1.6003345290 0.0012279710 0.077%
59 1.6640625000 1.6642794824 0.0002169824 0.013%
60 1.7382812500 1.7320508076 0.0062304424 0.360%
61 1.8046875000 1.8040477553 0.0006397447 0.035%
62 1.8867187500 1.8807264653 0.0059922847 0.319%
63 1.9648437500 1.9626105055 0.0022332445 0.114%
64 2.0585937500 2.0503038416 0.0082899084 0.404%
65 2.1445312500 2.1445069205 0.0000243295 0.001%
66 2.2500000000 2.2460367739 0.0039632261 0.176%
67 2.3632812500 2.3558523658 0.0074288842 0.315%
68 2.4882812500 2.4750868534 0.0131943966 0.533%
69 2.6054687500 2.6050890647 0.0003796853 0.015%
70 2.7500000000 2.7474774195 0.0025225805 0.092%
71 2.9062500000 2.9042108777 0.0020391223 0.070%
72 3.0820312500 3.0776835372 0.0043477128 0.141%
73 3.2773437500 3.2708526185 0.0064911315 0.198%
74 3.4921875000 3.4874144438 0.0047730562 0.137%
75 3.7343750000 3.7320508076 0.0023241924 0.062%
76 4.0468750000 4.0107809335 0.0360940665 0.900%
77 4.3671875000 4.3314758743 0.0357116257 0.824%
78 4.7382812500 4.7046301095 0.0336511405 0.715%
79 5.1718750000 5.1445540160 0.0273209840 0.531%
80 5.6835937500 5.6712818196 0.0123119304 0.217%
81 6.4023437500 6.3137515147 0.0885922353 1.403%
82 7.1953125000 7.1153697224 0.0799427776 1.124%
83 8.2031250000 8.1443464280 0.0587785720 0.722%
84 9.5273437500 9.5143644542 0.0129792958 0.136%
85 11.6640625000 11.4300523028 0.2340101972 2.047%
86 14.5000000000 14.3006662567 0.1993337433 1.394%
87 19.1328125000 19.0811366877 0.0516758123 0.271%
88 30.0898437500 28.6362532829 1.4535904671 5.076%
89 60.2109375000 57.2899616308 2.9209758692 5.099%

6. Výpis celého demonstračního příkladu

V této kapitole je uveden výpis celého demonstračního příkladu, který byl odladěn pro třicetidvoubitové překladače jazyka C. Otestování bylo provedeno pomocí GCC 3.4.2 (uznávám, je to poněkud starší, zato však stabilní verze) a Borland C++ Compileru 5.5 (pro Windows) běžícím v režimu překladače C:

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

/* počet míst před a za binární řádovou tečkou */
#define A 16
#define B 16

/* Ludolfovo číslo */
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* maximální počet iterací při běhu algoritmu */
#define MAX_ITER 16

/* "zesílení" při rotacích */
#define K 0.6073

/* převody mezi stupni a radiány */
#define rad2deg(rad) ((rad)*180.0/M_PI)
#define deg2rad(deg) ((deg)/180.0*M_PI)

/* datový typ, se kterým budeme pracovat */
typedef signed int fx;

/* hlavičky použitých funkcí */
void   fx_print(fx x);
fx     fp2fx(double x);
double fx2fp(fx x);

/* tabulka arkustangentu úhlů */
fx atans[MAX_ITER];

/* tabulka záporných celočíselných mocnin hodnoty 2 */
fx pows[MAX_ITER];

/*
 * Tisk numerické hodnoty uložené ve formátu pevné
 * řádové binární čárky (FX)
 */
void fx_print(fx x)
{
    int i;
    int val=x;                         /* pomocná proměnná pro převod do dvojkové soustavy */
    printf("bin: ");
    for (i=0; i<A+B; i++) {            /* převod na řetězec bitů (do dvojkové soustavy) */
        putchar(!!(val & (1<<(A+B-1)))+'0'); /* výpis hodnoty aktuálně nejvyššího bitu */
        if (i==B-1) putchar('.');      /* po řádové binární čárce vypsat značku */
        val=val<<1;                    /* posun na další (méně významný) bit */
    }

    printf("   hex: %08x   fp: %+11.5f\n", x, fx2fp(x));
}

/*
 * Převod z formátu plovoucí řádové binární čárky (FP)
 * do formátu pevné řádové binární čárky (FX)
 */
fx fp2fx(double x)
{
    return (fx)(x*(2<<(B-1)));
}

/*
 * Převod z celočíselného formátu (integer)
 * do formátu pevné řádové binární čárky (FX)
 */
fx int2fx(int x)
{
    return (fx)(x<<B);
}

/*
 * Převod z formátu pevné řádové binární čárky (FX)
 * do formátu plovoucí řádové binární čárky (FP)
 */
double fx2fp(fx x)
{
    return (double)x/(2<<(B-1));
}

/*
 * Součet dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_add(fx x, fx y)
{
    return x+y;
}

/*
 * Rozdíl dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_sub(fx x, fx y)
{
    return x-y;
}

/*
 * Součin dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_mul(fx x, fx y)
{
    fx result=(x>>(B/2))*(y>>(B/2));
    return result;
}

/*
 * Podíl dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_div(fx x, fx y)
{
    fx result=x/(y>>(B/2));
    return result<<(B/2);
}

/*
 * Vytvoření tabulky pro výpočet goniometrických
 * funkcí pomocí algoritmu CORDIC
 */
void fx_create_tables(void)
{
    int i;
    for (i=0; i<MAX_ITER; i++) {
        double p=pow(2.0, -i);
        atans[i]=fp2fx(atan(p));
        pows[i]=fp2fx(p);
    }
}

/* výpočet funkce tan() pro zadaný úhel delta */
/* (neoptimalizovaná verze) */
fx fx_tan_cordic(fx delta)
{
    int i;
    /* nastavení počátečních podmínek */
    fx x0=fp2fx(1.0);
    fx y0=fp2fx(0.0);
    fx xn;
    if (delta==0) return 0;                 /* ošetření nulového úhlu */
    for (i=0; i<MAX_ITER; i++) {            /* iterační smyčka */
        if (delta<0) {                      /* úhel je záporný => rotace doleva */
            xn=fx_add(x0, fx_mul(y0, pows[i]));
            y0=fx_sub(y0, fx_mul(x0, pows[i]));
            delta=fx_add(delta, atans[i]);
        }
        else {                              /* úhel je kladný => rotace doprava */
            xn=fx_sub(x0, fx_mul(y0, pows[i]));
            y0=fx_add(y0, fx_mul(x0, pows[i]));
            delta=fx_sub(delta, atans[i]);
        }
        x0=xn;
/*        printf("%i\t%+f\t%+f\t%+f\n", i, fx2fp(x0), fx2fp(y0), fx2fp(delta)); */
    }
    if (x0==0)                              /* ošetření tangenty pravého úhlu */
        if (y0<0) return 0;
        else return 0;
    else
        return fx_div(y0,x0);               /* vrátit výsledek operace */
}

/* výpočet funkce tan() pro zadaný úhel delta */
/* (optimalizovaná verze) */
fx fx_tan_cordic_optim(fx delta)
{
    int i;
    /* nastavení počátečních podmínek */
    fx x0=int2fx(1);
    fx y0=0;
    fx xn;
    if (delta==0) return 0;                 /* ošetření nulového úhlu */
    for (i=0; i<MAX_ITER; i++) {            /* iterační smyčka */
        if (delta<0) {                      /* úhel je záporný => rotace doleva */
            xn=fx_add(x0, y0>>i);           /* místo násobení bitový posuv */
            y0=fx_sub(y0, x0>>i);
            delta=fx_add(delta, atans[i]);
        }
        else {                              /* úhel je kladný => rotace doprava */
            xn=fx_sub(x0, y0>>i);
            y0=fx_add(y0, x0>>i);
            delta=fx_sub(delta, atans[i]);
        }
        x0=xn;
    }
    if (x0==0)                              /* ošetření tangenty pravého úhlu */
        if (y0<0) return 0;
        else return 0;
    else
        return fx_div(y0,x0);               /* vrátit výsledek operace */
}


int main(void)
{
    int i;
    fx     tanfx;
    double delta;                       /* úhel, ze kterého se počítá tan */
    double tanval;                      /* vypočtené hodnoty */
    double tanerr;                      /* absolutní chyby */

    fx_create_tables();

    /* kontrolní výpis tabulky atans[] */
    puts("\nHodnoty uložené v tabulce atans[]:");
    for (i=0; i<MAX_ITER; i++)
        printf("%d\t%f\n", i, fx2fp(rad2deg(atans[i])));
        /*printf("%f\n", fx2fp(pows[i])); */

    /* výpočet tan() neooptimalizovanou metodou CORDIC */
    puts("\nVýpočet tan() neoptimalizovanou metodou CORDIC");
    for (i=0; i<90; i++) {                  /* výpočetní smyčka */
        delta=deg2rad(i);                   /* převod úhlu na radiány */
        tanfx=fx_tan_cordic(fp2fx(delta));  /* aplikace algoritmu CORDIC */
        tanval=fx2fp(tanfx);                /* výpočet funkce tan */
        tanerr=fabs(tanval-tan(delta));     /* výpočet absolutních chyb */
                                            /* tisk výsledků */
        printf("%02d\t%14.10f\t%14.10f\t%12.10f\t%8.3f%%\n",
                i,
                tanval,
                tan(delta),
                tanerr,
                tanerr==0 ? 0:100.0*tanerr/tan(delta));
    }

    /* výpočet tan() optimalizovanou metodou CORDIC */
    puts("\nVýpočet tan() optimalizovanou metodou CORDIC");
    for (i=0; i<90; i++) {                  /* výpočetní smyčka */
        delta=deg2rad(i);                   /* převod úhlu na radiány */
        tanfx=fx_tan_cordic_optim(fp2fx(delta));  /* aplikace algoritmu CORDIC */
        tanval=fx2fp(tanfx);                /* výpočet funkce tan */
        tanerr=fabs(tanval-tan(delta));     /* výpočet absolutních chyb */
                                            /* tisk výsledků */
        printf("%02d\t%14.10f\t%14.10f\t%12.10f\t%8.3f%%\n",
                i,
                tanval,
                tan(delta),
                tanerr,
                tanerr==0 ? 0:100.0*tanerr/tan(delta));
    }
    return 0;
} 

7. Obsah poslední části tohoto seriálu

Seriál o numerických formátech s pevnou a plovoucí binární řádovou čárkou se blíží ke svému závěru. V následujícím a současně i posledním dílu si ukážeme implementaci dalších goniometrických funkcí algoritmem CORDIC pracujícím nad FX formátem.

Našli jste v článku chybu?
120na80.cz: Pánové, pečujte o svoje přirození a prostatu

Pánové, pečujte o svoje přirození a prostatu

Lupa.cz: Proč firmy málo chrání data? Chovají se logicky

Proč firmy málo chrání data? Chovají se logicky

Podnikatel.cz: EET zvládneme, budou horší zákony

EET zvládneme, budou horší zákony

Podnikatel.cz: Přehledná titulka, průvodci, responzivita

Přehledná titulka, průvodci, responzivita

DigiZone.cz: ČRo rozšiřuje DAB do Berouna

ČRo rozšiřuje DAB do Berouna

Vitalia.cz: Tesco: Chudá rodina si koupí levné polské kuře

Tesco: Chudá rodina si koupí levné polské kuře

Podnikatel.cz: Víme první výsledky doby odezvy #EET

Víme první výsledky doby odezvy #EET

Vitalia.cz: Jsou čajové sáčky toxické?

Jsou čajové sáčky toxické?

Měšec.cz: Finančním poradcům hrozí vracení provizí

Finančním poradcům hrozí vracení provizí

Vitalia.cz: Baletky propagují zdravotní superpostel

Baletky propagují zdravotní superpostel

120na80.cz: Rakovina oka. Jak ji poznáte?

Rakovina oka. Jak ji poznáte?

Podnikatel.cz: Na poslední chvíli šokuje vyjímkami v EET

Na poslední chvíli šokuje vyjímkami v EET

Lupa.cz: Google měl výpadek, nejel Gmail ani YouTube

Google měl výpadek, nejel Gmail ani YouTube

Lupa.cz: UX přestává pro firmy být magie

UX přestává pro firmy být magie

Měšec.cz: U levneELEKTRO.cz už reklamaci nevyřídíte

U levneELEKTRO.cz už reklamaci nevyřídíte

Vitalia.cz: To není kašel! Správná diagnóza zachrání život

To není kašel! Správná diagnóza zachrání život

DigiZone.cz: ČT má dalšího zástupce v EBU

ČT má dalšího zástupce v EBU

Podnikatel.cz: Prodává přes internet. Kdy platí zdravotko?

Prodává přes internet. Kdy platí zdravotko?

Měšec.cz: Kdy vám stát dá na stěhování 50 000 Kč?

Kdy vám stát dá na stěhování 50 000 Kč?

DigiZone.cz: Recenze Westworld: zavraždit a...

Recenze Westworld: zavraždit a...