The following script gives a closed form for the solution of a cubic equation

with real coefficients. This is the method published by Vieta in 1615. The

basis of the solution is due originally to Scipio del Ferro, who died in

1526, but who never published his solution. It was rediscovered by Tartaglia

in 1535, who also kept the method to himself. He divulged it under an oath of

secrecy to Cardano, who published it in 1545, crediting Tartaglia. After this

Cardano and Tartaglia debated bitterly and exchanged insults. Cardano's

method didn't handle the case where three real roots were obtained as the

sums of complex numbers. Bombelli in 1572 developed a method for working

freely with complex numbers, and so was able to solve the cubic equation in

greater generality. Vieta converted the solution to the trigonometric one

given below.

NB. Cubic equation: see Numerical Recipes in C, Press et al., Section 5.5, p.

157,

NB. or Handbook of Mathematical Functions by Abramowitz and Stegun, Section

3.8.2, p. 17.

NB. For quartic equation, see Abramowitz and Stegun, op. cit., Section 3.8.3,

pp 17-18.

NB. Uses arccos and cos from 'trig' utility script

NB. 'cubic' applies to a 4-item list of real number coefficients.

NB. The leading item is the constant term, then the linear term, etc.

NB. If last item is not 1, apply 'N' to argument before applying 'cubic'.

NB. For example, cubic _30 31 _10 1 is 2 5 3

NB. and cubic N _60 62 _20 2 is also 2 5 3

NB. and Cubic _60 62 _20 2 is also 2 5 3.

NB. Note that 2+5+3 is 10 (-a2, the quadratic term),

NB. (2*5)+(2*3)+(5*3) is 31 (a1, the linear term),

NB. and 2*5*3 is 30 (-a0, the constant term).

NB. This gives a simple way of checking the result.

N =: % {: NB. divide by tail to normalize coefficients with 1 as

high-order term

a0 =: {. NB. constant term

a1 =: 1&{ NB. linear term

a2 =: 2&{ NB. quadratic term

cubic =: (_2: * ([: %: Q) * [: cos (0p1 2p1 4p1"_ + theta) % 3:) - a2 % 3:

A simpler, more general, and much more efficient closed form is also

available in J, namely the polynomial primitive, which solves polynomials of

arbitrary order and with complex coefficients.

p.

Eugene McDonnell

Palo Alto