Metoda CORDIC a výpočet funkcí tan, atan a length

Obsah

1. Matematické funkce, jejichž hodnoty lze spočítat metodou CORDIC
2. Výpočet tangenty pomocí algoritmu CORDIC
3. První demonstrační příklad – výpočet tangenty
4. Výpočet arkustangenty pomocí algoritmu CORDIC
5. Druhý demonstrační příklad – výpočet arkustangenty
6. Výpočet délky vektoru pomocí algoritmu CORDIC
7. Třetí demonstrační příklad – výpočet délky vektoru
8. Převod z kartézských do polárních souřadnic pomocí algoritmu CORDIC
9. Obsah dalšího pokračování tohoto seriálu

1. Matematické funkce, jejichž hodnoty lze spočítat metodou CORDIC

V předchozím pokračování tohoto seriálu jsme si odvodili vztahy používané iterační metodou CORDIC pro výpočet goniometrických i jiných matematických funkcí. Jednalo se o upravené rovnice určené pro rotaci vektoru okolo počátku souřadné soustavy o úhly, které jsou vybrány s ohledem na hodnotu jejich tangenty (ta musí být rovna dvojce umocněné na záporný celočíselný exponent). Po několika krocích úprav jsme došli ke vztahům:

x_r=K_i (x₀ – y₀ d_i 2^-i)
y_r=K_i (y₀ + x₀ d_i 2^-i)

kde:

K_i=cos (arctan 2^-i)=1/(1+2^-2i)^1/2

Veškerá práce algoritmu CORDIC spočívá v tom, že se nastaví počáteční souřadnice vektoru r a iterativně se provádí rotace o předem známé úhly δ₁…δ_n tak, aby se dosáhlo požadované hodnoty rotace δ.

Také jsme si uvedli postup uplatňovaný při výpočtu goniometrických funkcí sinus a kosinus. Na tento výpočet v dalších kapitolách navážeme a ukážeme si, jakým způsobem je možné vypočítat další goniometrické a jiné matematické funkce.

2. Výpočet tangenty pomocí algoritmu CORDIC

Výpočet goniometrické funkce tangent je odvozen z minule popsaného algoritmu pro výpočet funkcí sin() a cos(), jelikož je možné použít známý vztah:

tan α=sin α / cos α

Funkce sin() a cos() se pomocí algoritmu CORDIC počítají současně, výpočet dokonce není možné žádným způsobem rozdělit, neboť obě hodnoty jsou na sobě závislé. To znamená, že tangentu je možné vypočítat s podobnou složitostí jako tyto dvě funkce. Jediný rozdíl spočívá v tom, že není zapotřebí obě vypočtené hodnoty násobit konstantou K_i, protože se hodnota této konstanty vzájemným vydělením vypočtených hodnot vyruší. Na druhou stranu se musí aplikovat dělení, které se v některých případech (jednodušší FPU bez děličky) opět řeší iteračními metodami. Vzhledem k tomu, že se vzájemně dělí dvě obecně nepřesné hodnoty, relativní chyba výsledku roste, přesněji řečeno, je součtem relativních chyb obou mezivýsledků. Proto se při výpočtu tangenty používá větší počet iterací, než je nutné pro samostatný výpočet funkcí sin() a cos(). Praktický příklad na výpočet tangenty pomocí algoritmu CORDIC je uveden ve třetí kapitole.

3. První demonstrační příklad – výpočet tangenty

V této kapitole je uveden výpis demonstračního příkladu, který slouží pro výpočet goniometrické funkce tangent. Všimněte si, že před vlastním iteračním výpočtem je proveden test na vstup nulové hodnoty. Je to z toho důvodu, že pro nulový úhel (a také úhel blízký úhlu pravému) je relativní chyba velká a test na nulovou vstupní hodnotu je implementačně velmi jednoduchý: postačuje hradlo AND s více vstupy. Na konci iteračního výpočtu je ošetřen stav, kdy je vypočtená tangenta nekonečná, tj. jedná se buď o úhel +90° nebo –90° (π/2 resp. -π/2). V tomto případě je vrácena hodnota +∞ nebo -∞ (jedná se o makra podle normy C99).

// --------------------------------------------------------
// Výpočet hodnot funkce tan() pomocí iteračního algoritmu
// CORDIC.
// --------------------------------------------------------

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

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

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

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

// tabulka arkustangentu úhlů
double atans[MAX_ITER];

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

// naplnění tabulek atans[] a pows[]
void createTables(void)
{
    int i;
    for (i=0; i<MAX_ITER; i++) {
        double p=pow(2.0, -i);
        atans[i]=atan(p);
        pows[i]=p;
    }
}

// výpočet funkce tan() pro zadaný úhel delta
double tan_cordic(double delta)
{
    int i;
    double x0=1.0;                          // nastavení počátečních podmínek
    double y0=0.0;
    double xn;
    if (delta==0) return 0.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=x0+y0*pows[i];
            y0-=x0*pows[i];
            delta+=atans[i];
        }
        else {                              // úhel je kladný => rotace doprava
            xn=x0-y0*pows[i];
            y0+=x0*pows[i];
            delta-=atans[i];
        }
        x0=xn;
    }
    if (x0==0)                              // ošetření tangenty pravého úhlu
        if (y0>0) return INFINITY;
        else return -INFINITY;
    else
        return y0/x0;
}

int main(void)
{
    int i;
    createTables();
    for (i=0; i<90; i++) {                  // výpočetní smyčka
        double delta;                       // úhel, ze kterého se počítá tan
        double tanval;                      // vypočtené hodnoty
        double tanerr;                      // absolutní chyby
        delta=i*M_PI/180.0;                 // převod úhlu na radiány
        tanval=tan_cordic(delta);           // 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;
}

// finito 

Výsledky běhu prvního demonstračního příkladu jsou vypsány v následující tabulce:

4. Výpočet arkustangenty pomocí algoritmu CORDIC

Výpočet arkustangenty se provádí poněkud odlišným způsobem než výše uvedený výpočet sinu, kosinu či tangenty, stále se však jedná o aplikaci velmi univerzálního principu, který CORDIC představuje. Využívá se zde takzvaný vektorový režim, ve kterém se – narozdíl od výše popsaného rotačního režimu – iterativně snažíme dosáhnout toho, aby se vynulovala hodnota uložená v registru y, tj. aby rotovaný vektor ležel na kladné či záporné horizontální poloose.

Výsledek běhu algoritmu CORDIC není v tomto případě uložen v registrech x nebo y (ty mají jiný význam), ale naopak v registru delta, který měl v rotačním režimu úlohu arbitru o směru rotace (ve vektorovém režimu tuto úlohu převzala hodnota v registru y,  resp. znaménko hodnoty v tomto registru uložené).

Vstupem algoritmu jsou počáteční hodnoty x a y (ty byly u předchozích metod nulové), přičemž se počítá arctan(y/x). To je velmi praktické, protože takto je možné vypočítat i arkustangentu kladného a záporného nekonečna (1/0 resp. –1/0), což jsou legální operace (stejným způsobem tento problém řeší standardní céčkovská knihovní funkce atan2()). Demonstrační příklad na výpočet arkustangenty je, spolu s tabulkou výsledků, uveden v páté kapitole.

5. Druhý demonstrační příklad – výpočet arkustangenty

V této kapitole je uveden výpis demonstračního příkladu, který provádí výpočet arkustangenty na základě čitatele a jmenovatele zlomku, tj. provádí se stejná operace, jaká je představována céčkovskou funkcí atan2(). Po překladu a spuštění příkladu se vypočte tabulka arkustangent pro hodnoty z intervalu 0,0 – 1,0, které odpovídají úhlům 0° – 45° (tj. 0 – π/4).

// --------------------------------------------------------
// Výpočet hodnot funkce atan() pomocí iteračního algoritmu
// CORDIC.
// --------------------------------------------------------

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

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

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

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

// tabulka arkustangentu úhlů
double atans[MAX_ITER];

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

// naplnění tabulek atans[] a pows[]
void createTables(void)
{
    int i;
    for (i=0; i<MAX_ITER; i++) {
        double p=pow(2.0, -i);
        atans[i]=atan(p);
        pows[i]=p;
    }
}

// výpočet funkce atan() pro zadané souřadnice x, y
double atan_cordic(double y, double x)
{
    int i;
    double x0=x;                            // nastavení počátečních podmínek
    double y0=y;
    double xn;
    double delta=0.0;
    for (i=0; i<MAX_ITER; i++) {            // iterační smyčka
        if (y0>0) {                         // kladná polorovina => rotace doleva
            xn=x0+y0*pows[i];
            y0-=x0*pows[i];
            delta+=atans[i];
        }
        else {                              // záporná polorovina => rotace doprava
            xn=x0-y0*pows[i];
            y0+=x0*pows[i];
            delta-=atans[i];
        }
        x0=xn;
    }
    return delta;                           // výsledek je uložen v akumulátoru úhlu
}

int main(void)
{
    double i;
    createTables();
    for (i=0.0; i<1.05; i+=0.05) {          // výpočetní smyčka
        double atanval;                     // vypočtená hodnoty
        double atanerr;                     // absolutní chyba
        double atanfloat;

        atanval=atan_cordic((double)i, 1.0)*180.0/M_PI; // výpočet funkce atan
        atanfloat=atan(i)*180.0/M_PI;
        atanerr=fabs(atanval-atanfloat);    // výpočet absolutní chyby
                                            // tisk výsledků
        printf("%3.2f\t%14.10f\t%14.10f\t%12.10f\t%8.3f%%\n",
                i,
                atanval,
                atanfloat,
                atanerr,
                atanfloat==0 ? 0:100.0*atanerr/atanfloat);
    }
    // důkaz, že atan se spočte i pro nekonečno, tj. pravý úhel:
    printf("\natan nekonecna: %f\n", atan_cordic(1.0, 0.0)*180.0/M_PI);
    return 0;
}

// finito 

V následující tabulce je ukázán výstup z předchozího programu spolu s absolutními a relativními chybami oproti hodnotě vypočtené pomocí standardní céčkovské funkce atan().

6. Výpočet délky vektoru pomocí algoritmu CORDIC

Délka vektoru se kupodivu počítá naprosto stejným způsobem jako funkce arkustangenty. Je to z toho důvodu, že se arkustangenta iterativně vypočítá pootočením vstupního vektoru do polohy, ve které má y-ovou souřadnici nulovou. To však logicky znamená, že jeho x-ová souřadnice odpovídá délce vektoru, protože platí:

d=sqrt(x²+y²)­=sqrt(x²+0)=sqrt(x²)=x

Proto se výše uvedená funkce pro výpočet arkustangenty může upravit na výpočet délky vektoru a to tak, že se výsledná hodnota x_i (vzniklá po i iteracích) vynásobí konstantou K_i stejným způsobem, jaký jsme prováděli u výpočtů funkcí sin() a cos() v předchozím pokračování tohoto seriálu.

Vlastnosti CORDICu, díky které je možné jednoduše spočítat délku vektoru, není široce známa, proto některé FPU pro tuto činnost ani neobsahují instrukci, což je škoda, zejména při práci s grafikou, kde se délky velmi často počítají a je pro ně nutné použít dvojici násobení, jedno sečítání a k tomu ještě druhou odmocninu, tj. poměrně zdlouhavé a složité operace.

7. Třetí demonstrační příklad – výpočet délky vektoru

Po překladu a spuštění třetího demonstračního příkladu se zobrazí tabulka s délkami vektorů, jejichž souřadnice postupně nabývají hodnot od nuly do desíti. Všimněte si, že přesnost výpočtů je poměrně velká, například pro vektor (6,8) je přesně vypočtena délka 10.

// --------------------------------------------------------
// Výpočet velikosti vektoru pomocí iteračního algoritmu
// CORDIC.
// --------------------------------------------------------

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

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

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

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

// tabulka arkustangentu úhlů
double atans[MAX_ITER];

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

// naplnění tabulek atans[] a pows[]
void createTables(void)
{
    int i;
    for (i=0; i<MAX_ITER; i++) {
        double p=pow(2.0, -i);
        atans[i]=atan(p);
        pows[i]=p;
    }
}

// výpočet velikosti vektoru pomocí algoritmu CORDIC
double mag_cordic(double y, double x)
{
    int i;
    double x0=x;                            // nastavení počátečních podmínek
    double y0=y;
    double xn;
    double delta=0.0;
    for (i=0; i<MAX_ITER; i++) {            // iterační smyčka
        if (y0>0) {                         // kladná polorovina => rotace doleva
            xn=x0+y0*pows[i];
            y0-=x0*pows[i];
            delta+=atans[i];
        }
        else {                              // záporná polorovina => rotace doprava
            xn=x0-y0*pows[i];
            y0+=x0*pows[i];
            delta-=atans[i];
        }
        x0=xn;
    }
    return x0*K;
}

int main(void)
{
    double x, y;
    createTables();
    printf("%f\n", mag_cordic(3, 4));       // výpočet Pythagorova trojúhelníka
    for (y=0.0; y<=10.0; y+=1.0) {          // tabulka velikostí různých vektorů
        for (x=0.0; x<=10.0; x+=1.0) {
            printf("%5.2f ", mag_cordic(x,y));
        }
        putchar('\n');
    }
    return 0;
}

// finito 

8. Převod z kartézských do polárních souřadnic pomocí algoritmu CORDIC

Vzhledem k tomu, že převod z kartézských souřadnic do souřadnic polárních vyžaduje spočtení úhlu vektoru a jeho délky (popř. úhlu a délky spojnice bodu s počátkem souřadné soustavy), je na tomto místě možné říci, že převod již máme vyřešen, protože jediným iteračním postupem CORDICu vyjádříme jak délku vektoru (ta je uložena v registru x), tak i jeho úhel (ten je uložen v registru delta). Demonstrační příklad na převod z kartézských souřadnic do souřadnic polárních zde proto nebudu uvádět, čtenář může „spojit“ obě výše uvedené céčkovské funkce atan_cordic() a mag_cordic() do funkce jediné a výsledkem bude zmíněná konverzní (mapovací) funkce.

9. Obsah dalšího pokračování tohoto seriálu

V následujícím pokračování tohoto seriálu konečně obrátíme svoji pozornost na formát pevné řádové binární čárky, protože si podrobně popíšeme čtyři základní způsoby uložení čísel v tomto formátu, včetně vysvětlení kódování záporných čísel. Také si popíšeme, jakým způsobem je možné implementovat základní matematické operace s hodnotami uloženými ve formátu pevné řádové binární čárky.