Solving a algebraic equation with Durand-Kerner-Aberth method

Last Update: Feb. 9, 2006


MPFR

Sample code

algebraic_mpf.c

#include <stdio.h>
#include "bnc.h"

#define MAX_LENGTH 1024

main()
{
    long int itimes;
    CMPFArray cmpfans, cmpfinit;
    MPFPoly mpff;
    mpf_t mpfabs_eps, mpfrel_eps;

    set_bnc_default_prec_decimal(100);

    /* init */
    mpf_init_set_d(mpfabs_eps, 1.0e-300);
    mpf_init_set_d(mpfrel_eps, 1.0e-25);
    mpff = init_mpfpoly(MAX_LENGTH);
    cmpfans = init_cmpfarray(20);
    cmpfinit = init_cmpfarray(20);

    // ff = (x-1)(x-2)...(x-20) 
    set_mpfpoly_i_str(mpff, 0, "2432902008176640000", 10);
    set_mpfpoly_i_str(mpff, 1, "-8752948036761600000", 10);
    set_mpfpoly_i_str(mpff, 2, "13803759753640704000", 10);
    set_mpfpoly_i_str(mpff, 3, "-12870931245150988800", 10);
    set_mpfpoly_i_str(mpff, 4, "8037811822645051776", 10);
    set_mpfpoly_i_str(mpff, 5, "-3599979517947607200", 10);
    set_mpfpoly_i_str(mpff, 6, "1206647803780373360", 10);
    set_mpfpoly_i_str(mpff, 7, "-311333643161390640", 10);
    set_mpfpoly_i_str(mpff, 8, "63030812099294896", 10);
    set_mpfpoly_i_str(mpff, 9, "-10142299865511450", 10);
    set_mpfpoly_i_str(mpff,10, "1307535010540395", 10);
    set_mpfpoly_i_str(mpff,11, "-135585182899530", 10);
    set_mpfpoly_i_str(mpff,12, "11310276995381", 10);
    set_mpfpoly_i_str(mpff,13, "-756111184500", 10);
    set_mpfpoly_i_str(mpff,14, "40171771630", 10);
    set_mpfpoly_i_str(mpff,15, "-1672280820", 10);
    set_mpfpoly_i_str(mpff,16, "53327946", 10);
    set_mpfpoly_i_str(mpff,17, "-1256850", 10);
    set_mpfpoly_i_str(mpff,18, "20615", 10);
    set_mpfpoly_i_str(mpff,19, "-210", 10);
    set_mpfpoly_i_str(mpff,20, "1", 10);

    print_mpfpoly(mpff);fflush(stdout);

    /* set Aberth's initial value */
    mpf_dka_init(cmpfinit, mpff);

    /* DKA method */
    itimes = mpf_dka(
        cmpfans,
        cmpfinit,
        mpff,
        1000,
        mpfabs_eps,
        mpfrel_eps
    );

    /* print answer */
    printf("Iterative times: %d\n", itimes);
    print_cmpfarray(cmpfans);

    /* clear */
    free_mpfpoly(mpff);
    free_cmpfarray(cmpfans);
    free_cmpfarray(cmpfinit);

    return 0;
}

Example execution

$ ./algebraic_mpf
-------------------------------------------------------------------------------
BNC Default Precision : 333 bits(100.2 decimal digits)
BNC Default Rounding Mode: Round to Nearest
-------------------------------------------------------------------------------
0 2.43290200817664000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e18
1 -8.75294803676160000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e18
2 1.38037597536407040000000000000000000000000000000000000000000000000000000000000000000000000000000000000e19
3 -1.28709312451509888000000000000000000000000000000000000000000000000000000000000000000000000000000000000e19
4 8.03781182264505177600000000000000000000000000000000000000000000000000000000000000000000000000000000000e18
5 -3.59997951794760720000000000000000000000000000000000000000000000000000000000000000000000000000000000000e18
6 1.20664780378037336000000000000000000000000000000000000000000000000000000000000000000000000000000000000e18
7 -3.11333643161390640000000000000000000000000000000000000000000000000000000000000000000000000000000000000e17
8 6.30308120992948960000000000000000000000000000000000000000000000000000000000000000000000000000000000000e16
9 -1.01422998655114500000000000000000000000000000000000000000000000000000000000000000000000000000000000000e16
10 1.30753501054039500000000000000000000000000000000000000000000000000000000000000000000000000000000000000e15
11 -1.35585182899530000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e14
12 1.13102769953810000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e13
13 -7.56111184500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e11
14 4.01717716300000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e10
15 -1.67228082000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e9
16 5.33279460000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e7
17 -1.25685000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e6
18 2.06150000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e4
19 -2.10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e2
20 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
ans->prec: 333
Iterative times: 143
0 (2.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000049249877108e1, -6.10914566986067981957943038530087397424424613934072483147443971467119376891552469614358781025000049365e-101)
1 (1.80000000000000000000000000000000000000000000000000000000000000000000000000000000000000003667709745089e1, -2.52695083787722710614922159880128403923778033391085292035781741138531507679327415606682087304874985586e-91)
2 (1.59999999999999999999999999999999999999999999999999999999999999999999999999999999999993920708592678750e1, -7.75965770080597591039635677913104515914294013093872965680154714166633432185322153895069962842292984721e-85)
3 (1.39999999999999999999999999999999999999999999999999999999999999999999999999999999649918157092808887228e1, -8.04143779896671793318300191021710720858148829850784813262850304315300619187553419537107798119152900286e-80)
4 (1.20000000000000000000000000000000000000000000000000000000000000000000000000025487333112413448084805549e1, -8.68871423159623493922182389746556695116007811703155930938717400675915257718224453540370967573402392270e-76)
5 (9.99999999999999999999999999999999999999999999999999999999999999999999999970050098118627535486795306303, -1.22087539938663158596561109757058325784625577471161609224220072203496571590862469941982436652525168120e-73)
6 (8.00000000000000000000000000000000000000000000000000000000000000000000000000003933034672441776909859553, -6.39953518295172508364563581627255385143419288736389241696643076738236142495694641515759357786212594767e-78)
7 (6.00000000000000000000000000000000000000000000000000000000000000000000000000000000032005680737796820423, 2.89183143467922120167266408537409366469463894819256874415439827244521767420488577196095418731280792088e-83)
8 (4.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000094928798168028, -1.94644249235878830991564458108015170768730699592778271652509023541731872468634907488941776545148774027e-88)
9 (2.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001438038, -9.66865623864486959919121052020096624611646850185255282228216623274744005063587468042422663139348734037e-96)
10 (1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000149, 6.10914566986067983035688094515758976520032898960291903122734907363773141140529548203977506926937991447e-101)
11 (3.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000017581024532, 2.52695083787722711080552830435161644525987525676551299410895236289684781074513955470813337349639011508e-91)
12 (5.00000000000000000000000000000000000000000000000000000000000000000000000000000000000061641179977117393, 7.75965770080597592043432022428989144465645848275508573135869366647628310214007913230462270476567928003e-85)
13 (7.00000000000000000000000000000000000000000000000000000000000000000000000000000003500819450215080175420, 8.04143779896671793752232124325929697625761744764772034239337575961235282431853343236405103633248527334e-80)
14 (8.99999999999999999999999999999999999999999999999999999999999999999999999999745126668875495797517288010, 8.68871423159623490306842055473971689243633310236201092228194750357621454598443476896702797567361995862e-76)
15 (1.10000000000000000000000000000000000000000000000000000000000000000000000002994990188137249795291497476e1, 1.22087539938663158600151407304377741642889409917639613242189558179648898501709559736290798452821875284e-73)
16 (1.29999999999999999999999999999999999999999999999999999999999999999999999999999606696532736955806537136e1, 6.39953518295172510864489031242169735701871947699999572125186616782965639730829171946049537344934602551e-78)
17 (1.49999999999999999999999999999999999999999999999999999999999999999999999999999999996799447544189410376e1, -2.89183143467922118784671690158375239961169789337402923717460860262878593676074967530456361464273660147e-83)
18 (1.69999999999999999999999999999999999999999999999999999999999999999999999999999999999999981903247187975e1, 1.94644249235878831412498607936740538493630567820526816133160961066112377238230469170693009974284627204e-88)
19 (1.89999999999999999999999999999999999999999999999999999999999999999999999999999999999999999002405459330e1, 9.66865623864486960467531007205987562440299278517241083179523641942955350056521883214077317793707010617e-96)

IEEE754 double

Sample code

algebraic_d.c

#include <stdio.h>
#include "bnc.h"

#define MAX_LENGTH 1024

main()
{
    long int itimes;
    CDArray cdans, cdinit;
    DPoly df;
    double dabs_eps, drel_eps;

    /* init */
    dabs_eps = 1.0e-100;
    drel_eps = 1.0e-14;
    df = init_dpoly(MAX_LENGTH);
    cdans = init_cdarray(10);
    cdinit = init_cdarray(10);

    /* ff = (x-1)(x-2)...(x-10) */
    set_dpoly_i(df, 0, (double)3628800);
    set_dpoly_i(df, 1, (double)-10628640);
    set_dpoly_i(df, 2, (double)12753576);
    set_dpoly_i(df, 3, (double)-8409500);
    set_dpoly_i(df, 4, (double)3416930);
    set_dpoly_i(df, 5, (double)-902055);
    set_dpoly_i(df, 6, (double)157773);
    set_dpoly_i(df, 7, (double)-18150);
    set_dpoly_i(df, 8, (double)1320);
    set_dpoly_i(df, 9, (double)-55);
    set_dpoly_i(df,10, (double)1);

    print_dpoly(df);

    /* set Aberth's initial value */
    ddka_init(cdinit, df);
    print_cdarray(cdinit);

    /* DKA method */
    itimes = ddka(
        cdans,
        cdinit,
        df,
        1000,
        dabs_eps,
        drel_eps
    );

    /* print answer */
    printf("Iterative times: %d\n", itimes);
    print_cdarray(cdans);

    /* clear */
    free_dpoly(df);
    free_cdarray(cdans);
    free_cdarray(cdinit);

    return 0;
}

Example execution

$ ./algebraic_d
0 3.62880000000000000e+06
1 -1.06286400000000000e+07
2 1.27535760000000000e+07
3 -8.40950000000000000e+06
4 3.41693000000000000e+06
5 -9.02055000000000000e+05
6 1.57773000000000000e+05
7 -1.81500000000000000e+04
8 1.32000000000000000e+03
9 -5.50000000000000000e+01
10 1.00000000000000000e+00
0 ( 6.03706502151305585e+02, 9.04100701465275307e+01)
1 ( 4.36317520495142901e+02, 4.24760243001174899e+02)
2 ( 1.04370888958790118e+02, 5.96866440099038073e+02)
3 ( -2.65341061661903836e+02, 5.40989943823221893e+02)
4 ( -5.31600932276856611e+02, 2.78473676578831260e+02)
5 ( -5.92706502151305585e+02, -9.04100701465274028e+01)
6 ( -4.25317520495143015e+02, -4.24760243001174786e+02)
7 ( -9.33708889587899193e+01, -5.96866440099038073e+02)
8 ( 2.76341061661904007e+02, -5.40989943823221893e+02)
9 ( 5.42600932276856724e+02, -2.78473676578831089e+02)
Iterative times: 1001
0 ( 9.99999999999513456e+00, 0.00000000000000000e+00)
1 ( 7.99999999995502709e+00, 0.00000000000000000e+00)
2 ( 5.99999999986361932e+00, 0.00000000000000000e+00)
3 ( 4.00000000001738076e+00, 0.00000000000000000e+00)
4 ( 2.00000000000000622e+00, 0.00000000000000000e+00)
5 ( 1.00000000000000089e+00, 0.00000000000000000e+00)
6 ( 3.00000000000059242e+00, 0.00000000000000000e+00)
7 ( 4.99999999998342304e+00, 0.00000000000000000e+00)
8 ( 7.00000000006232526e+00, 0.00000000000000000e+00)
9 ( 8.99999999990950705e+00, 0.00000000000000000e+00)

<- Go back