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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890

// rpoly_ak1.cpp - Program for calculating the roots of a polynomial of real coefficients.
// Written in Microsoft Visual Studio Express 2013 for Windows Desktop
// 27 May 2014
//
// The sub-routines listed below are translations of the FORTRAN routines included in RPOLY.FOR,
// posted off the NETLIB site as TOMS/493:
//
// http://www.netlib.org/toms/493
//
// TOMS/493 is based on the Jenkins-Traub algorithm.
//
// To distinguish the routines posted below from others, an _ak1 suffix has been appended to them.
//
// Following is a list of the major changes made in the course of translating the TOMS/493 routines
// to the C++ versions posted below:
// 1) All global variables have been eliminated.
// 2) The "FAIL" parameter passed into RPOLY.FOR has been eliminated.
// 3) RPOLY.FOR solves polynomials of degree up to 100, but does not explicitly state this limit.
//     rpoly_ak1 explicitly states this limit; uses the macro name MAXDEGREE to specify this limit;
//     and does a check to ensure that the user input variable Degree is not greater than MAXDEGREE
//     (if it is, an error message is output and rpoly_ak1 terminates). If a user wishes to compute
//     roots of polynomials of degree greater than MAXDEGREE, using a macro name like MAXDEGREE provides
//     the simplest way of offering this capability.
// 4) All "GO TO" statements have been eliminated.
//
// A small main program is included also, to provide an example of how to use rpoly_ak1. In this
// example, data is input from a file to eliminate the need for a user to type data in via
// the console.

#include <iostream>
#include <fstream>
#include <cctype>
#include <cmath>
#include <cfloat>

using namespace std;

#define MAXDEGREE 100
#define MDP1 MAXDEGREE+1

void rpoly_ak1(double op[MDP1], int* Degree, double zeror[MAXDEGREE], double zeroi[MAXDEGREE]);
void Fxshfr_ak1(int L2, int* NZ, double sr, double bnd, double K[MDP1], int N, double p[MDP1], int NN, double qp[MDP1], double* lzi, double* lzr, double* szi, double* szr);
void QuadSD_ak1(int NN, double u, double v, double p[MDP1], double q[MDP1], double* a, double* b);
int calcSC_ak1(int N, double a, double b, double* a1, double* a3, double* a7, double* c, double* d, double* e, double* f, double* g, double* h, double K[MDP1], double u, double v, double qk[MDP1]);
void nextK_ak1(int N, int tFlag, double a, double b, double a1, double* a3, double* a7, double K[MDP1], double qk[MDP1], double qp[MDP1]);
void newest_ak1(int tFlag, double* uu, double* vv, double a, double a1, double a3, double a7, double b, double c, double d, double f, double g, double h, double u, double v, double K[MDP1], int N, double p[MDP1]);
void QuadIT_ak1(int N, int* NZ, double uu, double vv, double* szr, double* szi, double* lzr, double* lzi, double qp[MDP1], int NN, double* a, double* b, double p[MDP1], double qk[MDP1], double* a1, double* a3, double* a7, double* d, double* e, double* f, double* g, double* h, double K[MDP1]);
void RealIT_ak1(int* iFlag, int* NZ, double* sss, int N, double p[MDP1], int NN, double qp[MDP1], double* szr, double* szi, double K[MDP1], double qk[MDP1]);
void Quad_ak1(double a, double b1, double c, double* sr, double* si, double* lr, double* li);

void rpoly_ak1(double op[MDP1], int* Degree, double zeror[MAXDEGREE], double zeroi[MAXDEGREE]){

int i, j, jj, l, N, NM1, NN, NZ, zerok;

double K[MDP1], p[MDP1], pt[MDP1], qp[MDP1], temp[MDP1];
double bnd, df, dx, factor, ff, moduli_max, moduli_min, sc, x, xm;
double aa, bb, cc, lzi, lzr, sr, szi, szr, t, xx, xxx, yy;

const double RADFAC = 3.14159265358979323846/180; // Degrees-to-radians conversion factor = pi/180
const double lb2 = log(2.0); // Dummy variable to avoid re-calculating this value in loop below
const double lo = FLT_MIN/DBL_EPSILON;
const double cosr = cos(94.0*RADFAC); // = -0.069756474
const double sinr = sin(94.0*RADFAC); // = 0.99756405

if ((*Degree) > MAXDEGREE){
    cout << "\nThe entered Degree is greater than MAXDEGREE. Exiting rpoly. No further action taken.\n";
    *Degree = -1;
    return;
} // End ((*Degree) > MAXDEGREE)

//Do a quick check to see if leading coefficient is 0
if (op[0] != 0){

N = *Degree;
xx = sqrt(0.5); // = 0.70710678
yy = -xx;

// Remove zeros at the origin, if any
j = 0;
while (op[N] == 0){
    zeror[j] = zeroi[j] = 0.0;
    N--;
    j++;
} // End while (op[N] == 0)

NN = N + 1;

// Make a copy of the coefficients
for (i = 0; i < NN; i++)   p[i] = op[i];

while (N >= 1){ // Main loop
    // Start the algorithm for one zero
    if (N <= 2){
    // Calculate the final zero or pair of zeros
        if (N < 2){
            zeror[(*Degree) - 1] = -(p[1]/p[0]);
            zeroi[(*Degree) - 1] = 0.0;
        } // End if (N < 2)
        else { // else N == 2
            Quad_ak1(p[0], p[1], p[2], &zeror[(*Degree) - 2], &zeroi[(*Degree) - 2], &zeror[(*Degree) - 1], &zeroi[(*Degree) - 1]);
        } // End else N == 2
        break;
    } // End if (N <= 2)

    // Find the largest and smallest moduli of the coefficients

    moduli_max = 0.0;
    moduli_min = FLT_MAX;

    for (i = 0; i < NN; i++){
        x = fabs(p[i]);
        if (x > moduli_max)   moduli_max = x;
        if ((x != 0) && (x < moduli_min))   moduli_min = x;
    } // End for i

    // Scale if there are large or very small coefficients
    // Computes a scale factor to multiply the coefficients of the polynomial. The scaling
    // is done to avoid overflow and to avoid undetected underflow interfering with the
    // convergence criterion.
    // The factor is a power of the base.

    sc = lo/moduli_min;

    if (((sc <= 1.0) && (moduli_max >= 10)) || ((sc > 1.0) && (FLT_MAX/sc >= moduli_max))){
        sc = ((sc == 0) ? FLT_MIN : sc);
        l = (int)(log(sc)/lb2 + 0.5);
        factor = pow(2.0, l);
        if (factor != 1.0){
            for (i = 0; i < NN; i++)   p[i] *= factor;
        } // End if (factor != 1.0)
    } // End if (((sc <= 1.0) && (moduli_max >= 10)) || ((sc > 1.0) && (FLT_MAX/sc >= moduli_max)))

    // Compute lower bound on moduli of zeros

    for (i = 0; i < NN; i++)   pt[i] = fabs(p[i]);
    pt[N] = -(pt[N]);

    NM1 = N - 1;

    // Compute upper estimate of bound

    x = exp((log(-pt[N]) - log(pt[0]))/(double)N);

    if (pt[NM1] != 0) {
        // If Newton step at the origin is better, use it
        xm = -pt[N]/pt[NM1];
        x = ((xm < x) ? xm : x);
    } // End if (pt[NM1] != 0)

    // Chop the interval (0, x) until ff <= 0

    xm = x;
    do {
        x = xm;
        xm = 0.1*x;
        ff = pt[0];
        for (i = 1; i < NN; i++)   ff = ff *xm + pt[i];
    } while (ff > 0); // End do-while loop

    dx = x;

    // Do Newton iteration until x converges to two decimal places

    while (fabs(dx/x) > 0.005) {
        df = ff = pt[0];
        for (i = 1; i < N; i++){
            ff = x*ff + pt[i];
            df = x*df + ff;
        } // End for i
        ff = x*ff + pt[N];
        dx = ff/df;
        x -= dx;
    } // End while loop

    bnd = x;

    // Compute the derivative as the initial K polynomial and do 5 steps with no shift

    for (i = 1; i < N; i++)   K[i] = (double)(N - i)*p[i]/((double)N);
    K[0] = p[0];

    aa = p[N];
    bb = p[NM1];
    zerok = ((K[NM1] == 0) ? 1 : 0);

    for (jj = 0; jj < 5; jj++) {
        cc = K[NM1];
        if (zerok){
            // Use unscaled form of recurrence
            for (i = 0; i < NM1; i++){
                j = NM1 - i;
                K[j] = K[j - 1];
            } // End for i
            K[0] = 0;
           zerok = ((K[NM1] == 0) ? 1 : 0);
        } // End if (zerok)

        else { // else !zerok
            // Used scaled form of recurrence if value of K at 0 is nonzero
            t = -aa/cc;
            for (i = 0; i < NM1; i++){
                j = NM1 - i;
                K[j] = t*K[j - 1] + p[j];
            } // End for i
            K[0] = p[0];
            zerok = ((fabs(K[NM1]) <= fabs(bb)*DBL_EPSILON*10.0) ? 1 : 0);
        } // End else !zerok

    } // End for jj

    // Save K for restarts with new shifts
    for (i = 0; i < N; i++)   temp[i] = K[i];

    // Loop to select the quadratic corresponding to each new shift

    for (jj = 1; jj <= 20; jj++){

        // Quadratic corresponds to a double shift to a non-real point and its
        // complex conjugate. The point has modulus BND and amplitude rotated
        // by 94 degrees from the previous shift.

        xxx = -(sinr*yy) + cosr*xx;
        yy = sinr*xx + cosr*yy;
        xx = xxx;
        sr = bnd*xx;

        // Second stage calculation, fixed quadratic

        Fxshfr_ak1(20*jj, &NZ, sr, bnd, K, N, p, NN, qp, &lzi, &lzr, &szi, &szr);

        if (NZ != 0){

            // The second stage jumps directly to one of the third stage iterations and
            // returns here if successful. Deflate the polynomial, store the zero or
            // zeros, and return to the main algorithm.

            j = (*Degree) - N;
            zeror[j] = szr;
            zeroi[j] = szi;
            NN = NN - NZ;
            N = NN - 1;
            for (i = 0; i < NN; i++)   p[i] = qp[i];
            if (NZ != 1){
                zeror[j + 1] = lzr;
                zeroi[j + 1] = lzi;
            } // End if (NZ != 1)
            break;
        } // End if (NZ != 0)
        else { // Else (NZ == 0)

            // If the iteration is unsuccessful, another quadratic is chosen after restoring K
            for (i = 0; i < N; i++)   K[i] = temp[i];
        } // End else (NZ == 0)

    } // End for jj

    // Return with failure if no convergence with 20 shifts

    if (jj > 20) {
        cout << "\nFailure. No convergence after 20 shifts. Program terminated.\n";
        *Degree -= N;
        break;
    } // End if (jj > 20)

} // End while (N >= 1)

} // End if op[0] != 0
else { // else op[0] == 0
    cout << "\nThe leading coefficient is zero. No further action taken. Program terminated.\n";
    *Degree = 0;
} // End else op[0] == 0

return;
} // End rpoly_ak1

void Fxshfr_ak1(int L2, int* NZ, double sr, double bnd, double K[MDP1], int N, double p[MDP1], int NN, double qp[MDP1], double* lzi, double* lzr, double* szi, double* szr){

// Computes up to L2 fixed shift K-polynomials, testing for convergence in the linear or
// quadratic case. Initiates one of the variable shift iterations and returns with the
// number of zeros found.

// L2 limit of fixed shift steps
// NZ number of zeros found

int fflag, i, iFlag, j, spass, stry, tFlag, vpass, vtry;
double a, a1, a3, a7, b, betas, betav, c, d, e, f, g, h, oss, ots, otv, ovv, s, ss, ts, tss, tv, tvv, u, ui, v, vi, vv;
double qk[MDP1], svk[MDP1];

*NZ = 0;
betav = betas = 0.25;
u = -(2.0*sr);
oss = sr;
ovv = v = bnd;

//Evaluate polynomial by synthetic division
QuadSD_ak1(NN, u, v, p, qp, &a, &b);

tFlag = calcSC_ak1(N, a, b, &a1, &a3, &a7, &c, &d, &e, &f, &g, &h, K, u, v, qk);

for (j = 0; j < L2; j++){

    //Calculate next K polynomial and estimate v
    nextK_ak1(N, tFlag, a, b, a1, &a3, &a7, K, qk, qp);
    tFlag = calcSC_ak1(N, a, b, &a1, &a3, &a7, &c, &d, &e, &f, &g, &h, K, u, v, qk);
    newest_ak1(tFlag, &ui, &vi, a, a1, a3, a7, b, c, d, f, g, h, u, v, K, N, p);

    vv = vi;

    // Estimate s

    ss = ((K[N - 1] != 0.0) ? -(p[N]/K[N - 1]) : 0.0);

    ts = tv = 1.0;

    if ((j != 0) && (tFlag != 3)){

       // Compute relative measures of convergence of s and v sequences

        tv = ((vv != 0.0) ? fabs((vv - ovv)/vv) : tv);
        ts = ((ss != 0.0) ? fabs((ss - oss)/ss) : ts);

        // If decreasing, multiply the two most recent convergence measures

        tvv = ((tv < otv) ? tv*otv : 1.0);
        tss = ((ts < ots) ? ts*ots : 1.0);

        // Compare with convergence criteria

        vpass = ((tvv < betav) ? 1 : 0);
        spass = ((tss < betas) ? 1 : 0);

        if ((spass) || (vpass)){

            // At least one sequence has passed the convergence test.
            // Store variables before iterating

            for (i = 0; i < N; i++)   svk[i] = K[i];

            s = ss;

            // Choose iteration according to the fastest converging sequence

            stry = vtry = 0;
            fflag = 1;

            do {

                iFlag = 1; // Begin each loop by assuming RealIT will be called UNLESS iFlag changed below

                if ((fflag && ((fflag = 0) == 0)) && ((spass) && (!vpass || (tss < tvv)))){
                    ; // Do nothing. Provides a quick "short circuit".
                } // End if (fflag)

                else { // else !fflag
                    QuadIT_ak1(N, NZ, ui, vi, szr, szi, lzr, lzi, qp, NN, &a, &b, p, qk, &a1, &a3, &a7, &d, &e, &f, &g, &h, K);

                    if ((*NZ) > 0)   return;

                    // Quadratic iteration has failed. Flag that it has been tried and decrease the
                    // convergence criterion

                    vtry = 1;
                    betav *= 0.25;

                    // Try linear iteration if it has not been tried and the s sequence is converging
                    if (stry || (!spass)){
                        iFlag = 0;
                    } // End if (stry || (!spass))
                    else {
                        for (i = 0; i < N; i++)   K[i] = svk[i];
                    } // End if (stry || !spass)

                } // End else !fflag

                if (iFlag != 0){
                    RealIT_ak1(&iFlag, NZ, &s, N, p, NN, qp, szr, szi, K, qk);

                    if ((*NZ) > 0)   return;

                    // Linear iteration has failed. Flag that it has been tried and decrease the
                    // convergence criterion

                    stry = 1;
                    betas *= 0.25;

                    if (iFlag != 0){

                        // If linear iteration signals an almost double real zero, attempt quadratic iteration

                        ui = -(s + s);
                        vi = s*s;
                        continue;

                    } // End if (iFlag != 0)
                } // End if (iFlag != 0)

                // Restore variables
                for (i = 0; i < N; i++)   K[i] = svk[i];

                // Try quadratic iteration if it has not been tried and the v sequence is converging

            } while (vpass && !vtry); // End do-while loop

            // Re-compute qp and scalar values to continue the second stage

            QuadSD_ak1(NN, u, v, p, qp, &a, &b);
            tFlag = calcSC_ak1(N, a, b, &a1, &a3, &a7, &c, &d, &e, &f, &g, &h, K, u, v, qk);

        } // End if ((spass) || (vpass))

    } // End if ((j != 0) && (tFlag != 3))

    ovv = vv;
    oss = ss;
    otv = tv;
    ots = ts;
} // End for j

return;
} // End Fxshfr_ak1

void QuadSD_ak1(int NN, double u, double v, double p[MDP1], double q[MDP1], double* a, double* b){

// Divides p by the quadratic 1, u, v placing the quotient in q and the remainder in a, b

int i;

q[0] = *b = p[0];
q[1] = *a = -((*b)*u) + p[1];

for (i = 2; i < NN; i++){
    q[i] = -((*a)*u + (*b)*v) + p[i];
    *b = (*a);
    *a = q[i];
} // End for i

return;
} // End QuadSD_ak1

int calcSC_ak1(int N, double a, double b, double* a1, double* a3, double* a7, double* c, double* d, double* e, double* f, double* g, double* h, double K[MDP1], double u, double v, double qk[MDP1]){

// This routine calculates scalar quantities used to compute the next K polynomial and
// new estimates of the quadratic coefficients.

// calcSC - integer variable set here indicating how the calculations are normalized
// to avoid overflow.

int dumFlag = 3; // TYPE = 3 indicates the quadratic is almost a factor of K

// Synthetic division of K by the quadratic 1, u, v
QuadSD_ak1(N, u, v, K, qk, c, d);

if (fabs((*c)) <= (100.0*DBL_EPSILON*fabs(K[N - 1]))) {
    if (fabs((*d)) <= (100.0*DBL_EPSILON*fabs(K[N - 2])))   return dumFlag;
} // End if (fabs(c) <= (100.0*DBL_EPSILON*fabs(K[N - 1])))

*h = v*b;
if (fabs((*d)) >= fabs((*c))){
    dumFlag = 2; // TYPE = 2 indicates that all formulas are divided by d
    *e = a/(*d);
    *f = (*c)/(*d);
    *g = u*b;
    *a3 = (*e)*((*g) + a) + (*h)*(b/(*d));
    *a1 = -a + (*f)*b;
    *a7 = (*h) + ((*f) + u)*a;
} // End if(fabs(d) >= fabs(c))
else {
    dumFlag = 1; // TYPE = 1 indicates that all formulas are divided by c;
    *e = a/(*c);
    *f = (*d)/(*c);
    *g = (*e)*u;
    *a3 = (*e)*a + ((*g) + (*h)/(*c))*b;
    *a1 = -(a*((*d)/(*c))) + b;
    *a7 = (*g)*(*d) + (*h)*(*f) + a;
} // End else

return dumFlag;
} // End calcSC_ak1

void nextK_ak1(int N, int tFlag, double a, double b, double a1, double* a3, double* a7, double K[MDP1], double qk[MDP1], double qp[MDP1]){

// Computes the next K polynomials using the scalars computed in calcSC_ak1

int i;
double temp;

if (tFlag == 3){ // Use unscaled form of the recurrence
    K[1] = K[0] = 0.0;

    for (i = 2; i < N; i++)   K[i] = qk[i - 2];

    return;
} // End if (tFlag == 3)

temp = ((tFlag == 1) ? b : a);

if (fabs(a1) > (10.0*DBL_EPSILON*fabs(temp))){
    // Use scaled form of the recurrence

    (*a7) /= a1;
    (*a3) /= a1;
    K[0] = qp[0];
    K[1] = -((*a7)*qp[0]) + qp[1];

    for (i = 2; i < N; i++)   K[i] = -((*a7)*qp[i - 1]) + (*a3)*qk[i - 2] + qp[i];

} // End if (fabs(a1) > (10.0*DBL_EPSILON*fabs(temp)))
else {
    // If a1 is nearly zero, then use a special form of the recurrence

    K[0] = 0.0;
    K[1] = -(*a7)*qp[0];

    for (i = 2; i < N; i++)   K[i] = -((*a7)*qp[i - 1]) + (*a3)*qk[i - 2];
} // End else

return;

} // End nextK_ak1

void newest_ak1(int tFlag, double* uu, double* vv, double a, double a1, double a3, double a7, double b, double c, double d, double f, double g, double h, double u, double v, double K[MDP1], int N, double p[MDP1]){
// Compute new estimates of the quadratic coefficients using the scalars computed in calcSC_ak1

double a4, a5, b1, b2, c1, c2, c3, c4, temp;

(*vv) = (*uu) = 0.0; // The quadratic is zeroed

if (tFlag != 3){

    if (tFlag != 2){
        a4 = a + u*b + h*f;
        a5 = c + (u + v*f)*d;
    } // End if (tFlag != 2)
    else { // else tFlag == 2
        a4 = (a + g)*f + h;
        a5 = (f + u)*c + v*d;
    } // End else tFlag == 2

    // Evaluate new quadratic coefficients

    b1 = -K[N - 1]/p[N];
    b2 = -(K[N - 2] + b1*p[N - 1])/p[N];
    c1 = v*b2*a1;
    c2 = b1*a7;
    c3 = b1*b1*a3;
    c4 = -(c2 + c3) + c1;
    temp = -c4 + a5 + b1*a4;
    if (temp != 0.0) {
        *uu= -((u*(c3 + c2) + v*(b1*a1 + b2*a7))/temp) + u;
        *vv = v*(1.0 + c4/temp);
    } // End if (temp != 0)

} // End if (tFlag != 3)

return;
} // End newest_ak1

void QuadIT_ak1(int N, int* NZ, double uu, double vv, double* szr, double* szi, double* lzr, double* lzi, double qp[MDP1], int NN, double* a, double* b, double p[MDP1], double qk[MDP1], double* a1, double* a3, double* a7, double* d, double* e, double* f, double* g, double* h, double K[MDP1]){

// Variable-shift K-polynomial iteration for a quadratic factor converges only if the
// zeros are equimodular or nearly so.

int i, j = 0, tFlag, triedFlag = 0;
double c, ee, mp, omp, relstp, t, u, ui, v, vi, zm;

*NZ = 0; // Number of zeros found
u = uu; // uu and vv are coefficients of the starting quadratic
v = vv;

do {
    Quad_ak1(1.0, u, v, szr, szi, lzr, lzi);

    // Return if roots of the quadratic are real and not close to multiple or nearly
    // equal and of opposite sign.

    if (fabs(fabs(*szr) - fabs(*lzr)) > 0.01*fabs(*lzr))   break;

    // Evaluate polynomial by quadratic synthetic division

    QuadSD_ak1(NN, u, v, p, qp, a, b);

    mp = fabs(-((*szr)*(*b)) + (*a)) + fabs((*szi)*(*b));

    // Compute a rigorous bound on the rounding error in evaluating p

    zm = sqrt(fabs(v));
    ee = 2.0*fabs(qp[0]);
    t = -((*szr)*(*b));

    for (i = 1; i < N; i++)   ee = ee*zm + fabs(qp[i]);

    ee = ee*zm + fabs((*a) + t);
    ee = (9.0*ee + 2.0*fabs(t) - 7.0*(fabs((*a) + t) + zm*fabs((*b))))*DBL_EPSILON;

    // Iteration has converged sufficiently if the polynomial value is less than 20 times this bound

    if (mp <= 20.0*ee){
        *NZ = 2;
        break;
    } // End if (mp <= 20.0*ee)

    j++;

    // Stop iteration after 20 steps
    if (j > 20)   break;

    if (j >= 2){
        if ((relstp <= 0.01) && (mp >= omp) && (!triedFlag)){
        // A cluster appears to be stalling the convergence. Five fixed shift
        // steps are taken with a u, v close to the cluster.

        relstp = ((relstp < DBL_EPSILON) ? sqrt(DBL_EPSILON) : sqrt(relstp));

        u -= u*relstp;
        v += v*relstp;

        QuadSD_ak1(NN, u, v, p, qp, a, b);

        for (i = 0; i < 5; i++){
            tFlag = calcSC_ak1(N, *a, *b, a1, a3, a7, &c, d, e, f, g, h, K, u, v, qk);
            nextK_ak1(N, tFlag, *a, *b, *a1, a3, a7, K, qk, qp);
        } // End for i

        triedFlag = 1;
        j = 0;

        } // End if ((relstp <= 0.01) && (mp >= omp) && (!triedFlag))

    } // End if (j >= 2)

    omp = mp;

    // Calculate next K polynomial and new u and v

    tFlag = calcSC_ak1(N, *a, *b, a1, a3, a7, &c, d, e, f, g, h, K, u, v, qk);
    nextK_ak1(N, tFlag, *a, *b, *a1, a3, a7, K, qk, qp);
    tFlag = calcSC_ak1(N, *a, *b, a1, a3, a7, &c, d, e, f, g, h, K, u, v, qk);
    newest_ak1(tFlag, &ui, &vi, *a, *a1, *a3, *a7, *b, c, *d, *f, *g, *h, u, v, K, N, p);

    // If vi is zero, the iteration is not converging
    if (vi != 0){
        relstp = fabs((-v + vi)/vi);
        u = ui;
        v = vi;
    } // End if (vi != 0)
} while (vi != 0); // End do-while loop

return;

} //End QuadIT_ak1

void RealIT_ak1(int* iFlag, int* NZ, double* sss, int N, double p[MDP1], int NN, double qp[MDP1], double* szr, double* szi, double K[MDP1], double qk[MDP1]){

// Variable-shift H-polynomial iteration for a real zero

// sss - starting iterate
// NZ - number of zeros found
// iFlag - flag to indicate a pair of zeros near real axis

int i, j = 0, nm1 = N - 1;
double ee, kv, mp, ms, omp, pv, s, t;

*iFlag = *NZ = 0;
s = *sss;

for ( ; ; ) {
    qp[0] = pv = p[0];

    // Evaluate p at s
    for (i = 1; i < NN; i++)   qp[i] = pv = pv*s + p[i];

    mp = fabs(pv);

    // Compute a rigorous bound on the error in evaluating p

    ms = fabs(s);
    ee = 0.5*fabs(qp[0]);
    for (i = 1; i < NN; i++)   ee = ee*ms + fabs(qp[i]);

    // Iteration has converged sufficiently if the polynomial value is less than
    // 20 times this bound

    if (mp <= 20.0*DBL_EPSILON*(2.0*ee - mp)){
        *NZ = 1;
        *szr = s;
        *szi = 0.0;
        break;
    } // End if (mp <= 20.0*DBL_EPSILON*(2.0*ee - mp))

    j++;

    // Stop iteration after 10 steps

    if (j > 10)   break;

    if (j >= 2){
        if ((fabs(t) <= 0.001*fabs(-t + s)) && (mp > omp)){
            // A cluster of zeros near the real axis has been encountered;
            // Return with iFlag set to initiate a quadratic iteration

            *iFlag = 1;
            *sss = s;
            break;
        } // End if ((fabs(t) <= 0.001*fabs(s - t)) && (mp > omp))

    } //End if (j >= 2)

    // Return if the polynomial value has increased significantly

    omp = mp;

    // Compute t, the next polynomial and the new iterate
    qk[0] = kv = K[0];
    for (i = 1; i < N; i++)   qk[i] = kv = kv*s + K[i];

    if (fabs(kv) > fabs(K[nm1])*10.0*DBL_EPSILON){
        // Use the scaled form of the recurrence if the value of K at s is non-zero
        t = -(pv/kv);
        K[0] = qp[0];
        for (i = 1; i < N; i++)   K[i] = t*qk[i - 1] + qp[i];
    } // End if (fabs(kv) > fabs(K[nm1])*10.0*DBL_EPSILON)
    else { // else (fabs(kv) <= fabs(K[nm1])*10.0*DBL_EPSILON)
        // Use unscaled form
        K[0] = 0.0;
        for (i = 1; i < N; i++)   K[i] = qk[i - 1];
    } // End else (fabs(kv) <= fabs(K[nm1])*10.0*DBL_EPSILON)

    kv = K[0];
    for (i = 1; i < N; i++)   kv = kv*s + K[i];

    t = ((fabs(kv) > (fabs(K[nm1])*10.0*DBL_EPSILON)) ? -(pv/kv) : 0.0);

    s += t;

} // End infinite for loop

return;

} // End RealIT_ak1

void Quad_ak1(double a, double b1, double c, double* sr, double* si, double* lr, double* li) {
// Calculates the zeros of the quadratic a*Z^2 + b1*Z + c
// The quadratic formula, modified to avoid overflow, is used to find the larger zero if the
// zeros are real and both zeros are complex. The smaller real zero is found directly from
// the product of the zeros c/a.

double b, d, e;

*sr = *si = *lr = *li = 0.0;

if (a == 0) {
    *sr = ((b1 != 0) ? -(c/b1) : *sr);
    return;
} // End if (a == 0))

if (c == 0){
    *lr = -(b1/a);
    return;
} // End if (c == 0)

// Compute discriminant avoiding overflow

b = b1/2.0;
if (fabs(b) < fabs(c)){
    e = ((c >= 0) ? a : -a);
    e = -e + b*(b/fabs(c));
    d = sqrt(fabs(e))*sqrt(fabs(c));
} // End if (fabs(b) < fabs(c))
else { // Else (fabs(b) >= fabs(c))
    e = -((a/b)*(c/b)) + 1.0;
    d = sqrt(fabs(e))*(fabs(b));
} // End else (fabs(b) >= fabs(c))

if (e >= 0) {
    // Real zeros

    d = ((b >= 0) ? -d : d);
    *lr = (-b + d)/a;
    *sr = ((*lr != 0) ? (c/(*lr))/a : *sr);
} // End if (e >= 0)
else { // Else (e < 0)
    // Complex conjugate zeros

    *lr = *sr = -(b/a);
    *si = fabs(d/a);
    *li = -(*si);
} // End else (e < 0)

return;
} // End Quad_ak1

int main()
{char rflag = 0; //Readiness flag

cout << "                                           rpoly_ak1 (27 May 2014)\n";
cout << "=========================================================================== \n";
cout << "This program calculates the roots of a polynomial of real coefficients:\n";
cout << "\nop[0]*x^N + op[1]*x^(N-1) + op[2]*x^(N-2) + . . . + op[N]*x^0 = 0 \n";
cout << "\n--------------------------------------------------------------------------- \n";
cout << "\nThis program can accept polynomials of degree 100 or less, specified by the\n";
cout << "constant MAXDEGREE. If a higher order polynomial is to be input, redefine\n";
cout << "MAXDEGREE and re-compile the program.\n";
cout << "\n--------------------------------------------------------------------------- \n";
cout << "\nAll relevant data for the polynomial whose roots will be sought should have\n";
cout << "been saved beforehand in a file named rpoly_ak1dat.txt.\n";
cout << "rpoly_ak1dat.txt should be in the same folder as the rpoly_ak1 executable. \n";
cout << "--------------------------------------------------------------------------- \n";
cout << "\nThe first entry of this file must be the degree, N, of the polynomial for\n";
cout << "which the roots are to be calculated.\n";
cout << "Entries for the coefficients of the polynomial should follow, starting with\n";
cout << "the coefficient for the highest power of x and working down to the coefficient\n";
cout << "for the x^0 term.\n";
cout << "\nThe data is assumed to be of type double. Variables used within this program\n";
cout << "are type double.\n";
cout << "\n--------------------------------------------------------------------------- \n";
cout << "\nThe output is written to the file rpoly_ak1out.txt.\n";
cout << "\nNote the returned value of the variable Degree.\n";
cout << "If Degree > 0, it specifies the number of zeros found.\n";
cout << "If Degree = 0, the leading coefficient of the input polynomial was 0.\n";
cout << "If Degree = -1, the input value of Degree was greater than MAXDEGREE.\n";
cout << "\n--------------------------------------------------------------------------- \n";
cout << "\nAdditional information is posted at the following URL:\n";
cout << "http://www.akiti.ca/rpoly_ak1_Intro.html\n";
cout << "--------------------------------------------------------------------------- \n";
cout << "\nIs everything ready (are you ready to continue?)? If yes, Enter y. \n";
cout << "Otherwise Enter any other key. \n";
cin >> rflag;

if (toupper(rflag) == 'Y') {

    int Degree; // The degree of the polynomial to be solved
    cout << "Appear to be ready. \n";

    ifstream in("rpoly_ak1dat.txt", ios::in);

    if (!in) {
        cout << "Cannot open the input file.\n";
        return 0;
    }

    in >> Degree; //Input the polynomial degree from the file
    if (Degree < 0) {
        cout << "Invalid polynomial degree entered. Program terminated. \n";
        in.close(); //Close the input file before terminating
        return 0;
    }

    ofstream out("rpoly_ak1out.txt", ios::out);
    if (!out) {
        cout << "Cannot open the output file. Program terminated.\n";
        in.close(); //Close the input file before terminating
        return 0;
    }

    double op[MDP1], zeroi[MAXDEGREE], zeror[MAXDEGREE]; // Coefficient vectors
    int i; // vector index

    //Input the polynomial coefficients from the file and put them in the op vector
    for (i = 0; i < (Degree+1); i++){
        in >> op[i];
    }//End for i

    in.close(); //Close the input file

    rpoly_ak1(op, &Degree, zeror, zeroi);

    out << "Degree = " << Degree << ".\n";
    out << "\n";

    if (Degree <= 0){
        cout << "\nReturned from rpoly_ak1 and Degree had a value <= 0.\n";
    } // End if (Degree <= 0)
    else { // else Degree > 0
        out.precision(DBL_DIG);
        out << "The roots follow:\n";
        out << "\n";
        for (i = 0; i < Degree; i++){
            out << zeror[i] << " + " << zeroi[i] << "i" << " \n";
        }//End for i
    } // End else Degree > 0

    out.close(); // Close the output file
} //End if rflag = 'Y'
else cout << "\nNot ready. Try again when ready with information. \n";
cout << "\nEnter any key to continue. \n";
cin >> rflag;
return 0;
} // End main program.