-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathpoly_mp.cpp
80 lines (79 loc) · 5.9 KB
/
poly_mp.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include "./rpoly.hpp"
#include<stdlib.h>
#include<stdio.h>
#include "./rpoly.hpp"
#include <complex>
#define WP 200
// N.B. you can use either CPP, GMP or MPC backend by
// defining CPP_MP, GMP_MP or MPC_MP
#define NDEG 20
#define MPC_MP
//#define GMP_MP
#ifdef CPP_MP
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_complex.hpp>
using namespace boost;
using namespace boost::multiprecision;
using namespace boost::multiprecision::backends;
using mpreal = number<cpp_bin_float<WP>>;
using mpcmplx = cpp_complex<WP>;
#elif defined(GMP_MP)
#include <boost/multiprecision/gmp.hpp>
#include <boost/multiprecision/complex_adaptor.hpp>
using namespace boost;
using namespace boost::multiprecision;
using namespace boost::multiprecision::backends;
using mpreal=number<gmp_float<WP>>;
using mpcmplx=number<complex_adaptor<gmp_float<WP>>>;
#elif defined(MPC_MP)
#include <boost/multiprecision/mpc.hpp>
#include <boost/multiprecision/mpfr.hpp>
using namespace boost;
using namespace boost::multiprecision;
using namespace boost::multiprecision::backends;
//we set 100 digits working precision!
using mpreal=number<mpfr_float_backend<WP>>;
using mpcmplx=number<mpc_complex_backend<WP>>;
#endif
int main(void)
{
rpoly<mpreal,NDEG,mpcmplx> P;
pvector<mpreal,NDEG+1> c;
pvector<mpcmplx,NDEG> r;
for (int i=0; i <= NDEG; i++)
c[i] = mpreal("1.0");
//c << mpreal("1E-06"), mpreal("-1000000000.000000003"), mpreal("3000000.000000000003"), mpreal("-3000.000000000000001"), mpreal("1.0");
#if 0
c << mpreal("3552713678800500903314963761379186629513463.280024451608251784379834516701253674307814857321148620648696049533158239403637926943122293988022731143994247458624314256025439881708959214876568165618129543"),
mpreal("-1000000000000000060171766412330792359273479264521267797403.871279328051461590840020033783734502413524223268837512539114455969113858472891070026756896810006098814698950731670032325433882912589604105128"),
mpreal("19000000000000000475577756347068607974614726050871621895.897999642271722782834130677208690431525889402826775405255708196627181945508929459369162872936886501258843643017715073901754016325587313251990268"),
mpreal("-171000000000000002321121151571240479911807900169495443.17038587502153766695111267745821610033109017784548996320381597329139171686108588316762685108510772852986859595390259291296733916290933590475643866"),
mpreal("969000000000000007789206735369702057959683200537013.65957142133955256393174072105413701280690677451714239062925283529236644229023570441116442883323607294256595343955536275006796314361714670816003566924"),
mpreal("-3876000000000000018881786593363326587085865770190.5572702623597528291962625882816849835515318093843077549713752112902939219644534644588909191308182489178693306475878957501133789778528223808182614971242"),
mpreal("11628000000000000033590556954774276329978334384.802486448049730157703252000761004768207061577300926579277495673352342000667457423256952394793929731661775950210897657518717903870694275722326688333297713"),
mpreal("-27132000000000000042943850594110992102667489.59554635755440191461070699208915859314133797833600589512908740609462425934464143373299767778165435118653922586670041114052572351316939245783123176643307621"),
mpreal("50388000000000000035257857408318920806100.036184633594502758921052057418107480539890151078871754226530185330557565454595079821427126979007736018764880153814906746915408652860325582925406733930668832265"),
mpreal("-75582000000000000007455480908909431990.717078402121917147771083775112964228820726244164039555823265035356904092999785356395945174490107904980608884915506235629953140012776256447894752460407424023144006"),
mpreal("92377999999999999971818024985429071.871341500121210893777630301850749565252475325383696836617238224747084098613553747969301812427039045667703260210103317674141344629097881312819427894298393705555461446"),
mpreal("-92377999999999999947784101959776.218694905434949670173498983025271340482542945891702795085256494216919815151507075581165150120702176910783704799621913182262290689921633609359077782132280683634820737441"),
mpreal("75581999999999999945750788984.619100868006821556698111985007739171511504320775588725426985343657356764920508842195727967758536696634557199072415316805887052116437194994130183836294736887578713616958976"),
mpreal("-50387999999999999960321941.279926542937695676940132270546371448067299700505531439040283618043415744466914358208900804841897952109539863063015758641219006128764429754224674645027768707655971729250153195"),
mpreal("27131999999999999978509.284078651310228072972892515132173571232693029603499782174348663454315038527362148844897185884527773068008622311426087159629837801884342985454260521002848830154355118670111051582"),
mpreal("-11627999999999999991.341150155301843541324800917039618245728562597519580084054303521536393390121476315675133827476425792437402639672510349162484855856958946906690461566221719010113125726377175611042131"),
mpreal("3875999999999999.9974614680710966459181787744982270831603612460457822615650269345631807164104951930040477763098215608341954829249017252563137240076790073008706573241397553810076733991175617631812093514"),
mpreal("-968999999999.99999948605563588844075603802421434314061374046209096282057086914180815041406701662512141239940402720710588149512418949328046963083941377462521749868004512657045736992939160887496799033111"),
mpreal("170999999.99999999993556527716939247026355467408651069635060047699798973873052157440860323522696521311758736742082743419055175883953711399423348158914520498198390679470214384330602508824246504019303109"),
mpreal("-18999.99999999999999622292019392177607415503404502562359995833823385258878281556219946231999877439424863678071810253744633456562951323576271533966064453125"),
mpreal("1");
#endif
P.set_coeff(c);
P.find_roots(r);
int cc=0;
for (auto& r0: r)
{
cout << setprecision(WP) << "root #" << cc << "=" << r0 << "\n";
cout << setprecision(WP) << "p(#" << cc << ")=" << P.evalpoly(r0) << "\n\n";
cc++;
}
return 0;
}