This is a sub-page of a longer article, The classical Chinese version of Horner’s method: Technical considerations.

Appendix 3
Computer program to illustrate classical
Chinese root extraction:
Wang Xiaotong cubics with fractional coefficients

100 ' Initialization.

110 '

120 epsilon = 0.001

130 root = 0

140 '

150 ' Subroutine to reduce the roots of the equation by the value of a

151 ' proposed digit.

160 '

170 sub chu(dv)

180     a = a+dv

190     b = b+a*dv

200     c = c-b*dv

210     a = a+dv

220     b = b+a*dv

230     a = a+dv

240     end sub

250 '

260 ' Ask for the coefficients of the equation.

270 '

280 print "coefficients: integer, numerator, denominator"

290 print "a:"

300 input a,anum,aden

310 print "b:"

320 input b,bnum,bden

330 print "c:"

340 input c,cnum,cden

350 '

360 ' Form improper fractions using a common multiple of the denominators.

370 '

380 if aden <> bden or aden <> cden then

390     anum = anum*bden*cden

400     bnum = bnum*aden*cden

410     cnum = cnum*aden*bden

420     divisor = aden*bden*cden

430 else

440     divisor = aden

450 endif

460 bnum = bnum*divisor

470 cnum = cnum*divisor^2

480 a = a*divisor+anum

490 b = b*divisor^2+bnum

500 c = c*divisor^3+cnum

510 '

520 ' Propose a digit.

530 '

540 cq = c^(1/3)

550 order = 10^floor(log10(cq)+1)

560 d = 1

570 '

580 ' Reduce the roots of the equation by the value of the proposed digit.

590 '

600 chu(d*order)

610 '

620 ' If the proposed digit was too large, undo the previous action, propose a

621 ' smaller digit, and go back to try again.

630 '

640 if c < 0

650     chu(-d*order)

660     d = d-1

670     if d = 0

680         order = order/10

690         d = 9

700         endif

710     goto 600

720     endif

730 '

740 ' Add the value of the digit to the accumulated result.

750 '

760 root = root+d*order

770 '

780 ' If the desired precision has not been reached, continue proposing digits.

790 '

800 if c > epsilon then goto 540

810 '

820 ' Print the final result and exit.

830 '

840 print "root=", root/divisor

850 end