Phrases/Polynomials

Polynomials in J are represented as rank 1 arrays of numbers. These numbers represent coefficients from zero power up. For example, _1 _4 0 3 represents the polynomial $-1-4x+3x^{3}$. Primitives p. and p.. and t. provide some polynomial-related functions.

Operations on polynomials

eval=: ([: +`*/ [: }:@, ,"0)"1 0	Evaluation of a polynomial, P&eval is equivalent to P&p.
normalize=: #~ 1: 0} +./\.@(0&~:)	Normalization (trimming nonessential trailing zeroes)
ppr=: +//.@(*/)	Polynomial product (from JPhrases)
ppr=: ([: p. 1 ; ])@,&(1 {:: p.)	Polynomial product (alternative form)
divmod=: >:@-&# ({. ,&< }.) ([ (] -/@,:&}. (* {:)) ] , %&{.~)^:(>:@-~&#)&.|.~	Division with remainder
atop=: [ +/ .* [:ppr/\ 1 0, 1}. ]"0 1	Composition: f atop g is equivalent to f(g(x)) (from JPhrases)
atop=: 4 : 'x&p. @ (y&p.) t. i. x+&#y'	Composition, alternative
roots=: 1 {:: p.	roots

Lagrange interpolation

Lagrange interpolation polynomial is a polynomial that takes specified values y in specified points x. The name is slightly misleading -- Lagrange polynomial is not used to interpolate experimental data.

lagrange=: ] +/ .* [ (] % p.~) 1 p.@;~&1\. [

For example, polynomial ${141 \over 130}+{47 \over 520}x+{3 \over 520}x^{2}$ takes values y=1 2 3 in points x=_1 7 12

   ]p=._1 7 12 lagrange 1 2 3
141r130 47r520 3r520
   p p. _1 7 12
1 2 3

Alternatively,

   lagrange1=: %.@(^/ i.@#)@[ +/ .* ]
   _1 7 12 lagrange1&x: 1 2 3
141r130 47r520 3r520

The lagrange polynomial approach may also be used directly, to evaluate a polynomial directly from pre-determined X and Y values. This is not very efficient, but can be illustrative.

   lagrangep=: ({:@[ +/ .* ] (- */@(% #~ 0~:])"1 -/~@]) {.@[)"2 0
   (_1 7 12,:1 2 3) lagrangep 7 12
2 3

Roots

Improvements to p. for J6.01 have improved the robustness of J's root finding facilities:

   p.1 5e4 1            NB. J version < J6.01
|limit error
|       p.1 50000 1
   p.1 5e4 1            NB. J version >: J6.01
┌─┬────────────┐
│1│_50000 _2e_5│
└─┴────────────┘

Nevertheless it may be necessary to use specialized polynomial root finding libraries or write your own routines to solve some specialized cases. For example, to solve quadratic equations you may use

quadratic=: 3 : 0
  'c b a'=.y.
  assert. 0~:a
  t=.b%2*a
  (-t) +`-`:0 %:(*:t)-c%a
)
   ]r=.quadratic 1 5e4 1
_2e_5 _50000
   0j22":r
_0.0000199999994947575030 _49999.9999800000010000000000
   1 5e4 1 p. r
2.56621e_8 2.56621e_8

You may also use the functions in GSL (GNU Scientific Library) as shown in the example at the Jprogramming forum thread Roots robustness. For an even more robust algorithm (especially if there are roots close to each other), you generate the companion matrix of the polynomial and use the lapack library to find its eigenvalues.