用BigDecimals实现Java中的小错误

Question

我用Java实现了CORDIC算法，在第一次迭代时，我在http://en.wikipedia.org/wiki/CORDIC上以这个例子为例，用Java重写了它。

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

所以这些测试失败了：

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];
    }
}

Answer 1

你的atan值是错的。精度较高的岩浆告诉我，第一个条目应该是

0.78539816339744830961566084581987572104929234984377645524373614807695410157155224

而在你的桌子上

0.7853981633974482789994908671360463

麦格玛的全套桌子是

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] ];