2

我用 Java 实现了 CORDIC 算法，在第一次迭代中，我只是在http://en.wikipedia.org/wiki/CORDIC上举了一个例子，然后用 Java 重写了它。

0.32719469679615224417334408526762060但似乎我的实现确实返回了 sin 和 cos 并带有一个小错误：例如，当计算我的实现的 1/3 的正弦时， 0.3271946967961523696204423482988852它只正确到小数点后 14 位，而不是定义的 34MathContext.DECIMAL128

所以这些测试失败：

BigDecimal oneThird = BigDecimal.ONE.divide(new BigDecimal(3, MathContext.DECIMAL128),MathContext.DECIMAL128);
BigDecimal tmp = Cordic.sin(oneThird, MathContext.DECIMAL128);
// tmp = 0.3271946967961523696204423482988852
assertEquals(0, new BigDecimal("0.32719469679615224417334408526762060").compareTo(tmp));

tmp = Cordic.cos(oneThird, MathContext.DECIMAL128);
// tmp = 0.9449569463147379497146525865916989
assertEquals(0, new BigDecimal("0.944956946314737664388284007675880609").compareTo(tmp));

我的实现如下所示：

package net.objecthunter.math;

import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;

import javax.print.attribute.standard.PresentationDirection;

public class Cordic {
    // in GNU octae this table can be generated using the command
    // "printf("%.34f\n", cumprod(1./sqrt(1.+2.^-(2.*i))))"
    public static BigDecimal K[] = new BigDecimal[]{
            new BigDecimal("0.7071067811865474617150084668537602",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6324555320336757713306496953009628",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6135719910778962837838435007142834",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6088339125177524291387953780940734",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6076482562561682509993943313020281",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6073517701412960434481647098436952",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072776440935261366149688910809346",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072591122988928447057332959957421",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072544793325624912228022367344238",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072533210898752864537186724191997",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072530315291344571448917122324929",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529591389449477034645497042220",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529410413972650317759871541057",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529365170102888527026152587496",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529353859135170523586566559970",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529351031393796134238982631359",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529350324458174981145930360071",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529350147723992137116511003114",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529350103540446426109156163875",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529350092494837554113473743200",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529350089733712891870709427167",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529350089043154170553862059023",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529350088871069601736962795258",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529350088827770903776581690181",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529350088816668673530330124777",
                    MathContext.DECIMAL128),
            new BigDecimal("0.6072529350088814448227481079811696",
                    MathContext.DECIMAL128)};

    // In GNU octave this table is generated using this command:
    // printf("%.34f\n",atan(2.^-(0:40)))
    public static BigDecimal A[] = new BigDecimal[]{
            new BigDecimal("0.7853981633974482789994908671360463",
                    MathContext.DECIMAL128),
            new BigDecimal("0.4636476090008060935154787784995278",
                    MathContext.DECIMAL128),
            new BigDecimal("0.2449786631268641434733268624768243",
                    MathContext.DECIMAL128),
            new BigDecimal("0.1243549945467614381566789916178095",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0624188099959573500230547438150097",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0312398334302682774421544564802389",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0156237286204768312941615349132007",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0078123410601011111439873069173245",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0039062301319669717573901390750279",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0019531225164788187584341550007139",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0009765621895593194594364927496599",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0004882812111948982899262139412144",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0002441406201493617712447448120372",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0001220703118936702078530659454358",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000610351561742087725935014541623",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000305175781155260957271547345160",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000152587890613157615423778681873",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000076293945311019699810389967098",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000038146972656064961417507561819",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000019073486328101869647792853193",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000009536743164059608441276310632",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000004768371582030888422810640821",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000002384185791015579736676881098",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000001192092895507806808997385635",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000596046447753905522081060953",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000298023223876953025738326494",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000149011611938476545956387749",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000074505805969238281250000000",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000037252902984619140625000000",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000018626451492309570312500000",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000009313225746154785156250000",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000004656612873077392578125000",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000002328306436538696289062500",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000001164153218269348144531250",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000000582076609134674072265625",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000000291038304567337036132812",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000000145519152283668518066406",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000000072759576141834259033203",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000000036379788070917129516602",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000000018189894035458564758301",
                    MathContext.DECIMAL128),
            new BigDecimal("0.0000000000009094947017729282379150",
                    MathContext.DECIMAL128)};

    public static BigDecimal[] apply(BigDecimal arg, MathContext mc) {
        // TODO: reduce the arg to 0 <= arg <= 2*pi

        if (arg.compareTo(BigDecimal.ZERO) == 0) {
            return new BigDecimal[]{BigDecimal.ONE, BigDecimal.ZERO};
        }

        if (arg.compareTo(new BigDecimal(Math.PI)) == 0) {
            return new BigDecimal[]{BigDecimal.ONE.negate(), BigDecimal.ZERO};
        }

        if (arg.compareTo(new BigDecimal(2 * Math.PI)) == 0) {
            return new BigDecimal[]{BigDecimal.ONE, BigDecimal.ZERO};
        }

        if (arg.compareTo(new BigDecimal(Math.PI / 2d)) == 0) {
            return new BigDecimal[]{BigDecimal.ZERO, BigDecimal.ONE};
        }

        if (arg.compareTo(new BigDecimal(3d * Math.PI / 2d)) == 0) {
            return new BigDecimal[]{BigDecimal.ZERO, BigDecimal.ONE.negate()};
        }

        int i = 0;
        BigDecimal beta = arg;
        BigDecimal sigma;
        BigDecimal factor;
        BigDecimal powerOfTwo = BigDecimal.ONE;
        BigDecimal x = BigDecimal.ONE;
        BigDecimal y = BigDecimal.ZERO;
        BigDecimal angle = A[0];
        BigDecimal two = BigDecimal.ONE.add(BigDecimal.ONE);
        BigDecimal precision = new BigDecimal(BigInteger.ONE,
                mc.getPrecision() * 2);
        BigDecimal tmpx;
        BigDecimal tmpy;
        for (; i < Integer.MAX_VALUE; i++) {
            if (beta.compareTo(BigDecimal.ZERO) < 0) {
                sigma = BigDecimal.ONE.negate();
            } else {
                sigma = BigDecimal.ONE;
            }
            factor = sigma.multiply(powerOfTwo, mc);
            tmpx = x.subtract(y.multiply(factor, mc), mc);
            tmpy = y.add(x.multiply(factor, mc), mc);
            x = tmpx;
            y = tmpy;
            beta = beta.subtract(sigma.multiply(angle));
            powerOfTwo = powerOfTwo.divide(two, mc);
            if (i + 2 > A.length) {
                angle = angle.divide(two);
            } else {
                angle = A[i + 1];
            }
            if (beta.abs().compareTo(precision) < 0) {
                i++;
                break;
            }
        }
        return new BigDecimal[]{x.multiply(K[Math.min(K.length - 1, i)], mc),
                y.multiply(K[Math.min(K.length - 1, i)], mc)};
    }

    public static BigDecimal cos(BigDecimal arg, MathContext mc) {
        return apply(arg, mc)[0];
    }

    public static BigDecimal sin(BigDecimal arg, MathContext mc) {
        return apply(arg, mc)[1];
    }
}
4

1 回答 1

3

您的 atan 值是错误的。精度更高的岩浆告诉我第一个条目应该是

0.78539816339744830961566084581987572104929234984377645524373614807695410157155224

在你的桌子上是

0.7853981633974482789994908671360463

Magma 的完整表格是

0.78539816339744830961566084581987572104929234984377645524373614807695410157155224, 
0.4636476090008061162142562314612144020285370542861202638109330887201978641657417, 
0.24497866312686415417208248121127581091414409838118406712737591466735511958764209, 
0.12435499454676143503135484916387102557317019176980408991511411911572226742756675, 
0.06241880999595734847397911298550511360627388779749919460752781689869026721680345, 
0.03123983343026827625371174489249097703249566372540004025531558625579642101943244, 
0.01562372862047683080280152125657031891111413980090541788141050739666477417640177, 
0.00781234106010111129646339184219928162122281172501472355745390224838987204533523, 
0.00390623013196697182762866531142438714035749011520285621521309514901134416395438, 
0.00195312251647881868512148262507671393161074677723351033905753396043108530313709, 
0.00097656218955931943040343019971729085163419701581008759004900725226763752035508, 
0.00048828121119489827546923962564484866619236113313500303710940335348751213674327, 
0.00024414062014936176401672294325965998621241779097061761180790046091017847021746, 
0.00012207031189367020423905864611795630093082940901578749845193983784664259022045, 
6.1035156174208775021662569173829153785143536833346179337671134316586565776889807e-5, 
3.0517578115526096861825953438536019750949675119437837531021156883611630486111094e-5, 
1.525878906131576210723193581269788513742923814457587484624118640744586426707683e-5, 
7.629394531101970263388482340105090586350743918468077157763830696533368540109726e-6, 
3.814697265606496282923075616372993722805257303968866310187439250393888463610412e-6, 
1.907348632810187035365369305917244168714342165450153366670057723467064463709843e-6, 
9.53674316405960879420670689923112390019634124498790160133611802076003329888112e-7, 
4.7683715820308885992758382144924707587049404378664196740053215887142363814443e-7, 
2.3841857910155798249094797721893269783096898769063155913766911372217648282103e-7, 
1.19209289550780685311368497137922112645967587664586735576738225215437199588955e-7, 
5.9604644775390554413921062141788874250030195782366297314294565710005108461658e-8, 
2.9802322387695303676740132767709503349043907067445107249258477840843557260847e-8, 
1.4901161193847655147092516595963247108248930025964720012170057805491014206727e-8, 
7.4505805969238279871365645744953921132066925545665870075947601416173711836981947e-9, 
3.7252902984619140452670705718119235836719483287370405242319982692391073974358196e-9, 
1.8626451492309570290958838214764904345065282835738863513491050124951302594430928e-9, 
9.313225746154785153557354776845613038929264961492906739437685424219745532957262e-10, 
4.656612873077392577788419347105701629734786389156161742132349255441464969391566e-10, 
2.328306436538696289020427418388212703712742932049818605254866622806071463876315e-10, 
1.164153218269348144525990927298526587963964573800142900265849791708846857314265e-10, 
5.82076609134674072264967615912315823495491562577952724239762061671471623655995e-11, 
2.91038304567337036132730326989039477936936320036398304958299345250291480963431e-11, 
1.45519152283668518066395978373629934742117036089367107320672702133070941642532e-11, 
7.27595761418342590332018410467037418427646293888214296401117528908389857511e-12, 
3.6379788070917129516601402005837967730345578669779258118296083646485740987638e-12, 
1.8189894035458564758300761188229745966293197333602925371452076535033553005523e-12, 
9.0949470177292823791503881172787182457866496666966318622647928818549076982886781e-13

或者自己尝试一下（http://magma.maths.usyd.edu.au/calc/）：

RR:=RealField(80);
[ Arctan(RR!2^(-k)) : k in [0..40] ];
于 2014-02-17T17:24:16.353 回答