Polynomials.xtl

1⍝# Polynomials: polynomials -- evaluation, sums, products, derivatives,
2⍝# integrals, real roots, text.
3⍝# Import it with an alias of your choice: "py:" u_se< "Polynomials".
4⍝# Put libs/Polynomials/src on XETAL_PATH ("just path"); the reference is libs/Polynomials/docs.
5⍝# Names with l: are exported; those under h: are private to this file.
6⍝#
7⍝# A polynomial is its coefficients, highest power first, as APL's
8⍝# decode reads digits: 3 -2 1 is 3x^2 - 2x + 1. Results are Floats.
9
10ᵗ⁼u̲se< "Strings"
11ᶠ⁼u̲se< "Format"
12
13⍝## Evaluation
14
15⍝# p a_t x: p's value at each x: the coefficients decoded in radix x
16⍝# (APL's decode is Horner's rule, on any numbers).
17ˡa̲t ← { p x → c ← f̲loat p◆ '{ v → v d̲ecode c } e̲ach f̲loat x }
18
19⍝## Arithmetic
20
21⍝# p without its leading zeros (at least one coefficient kept).
22ʰt̲rim ← { p → n ← t̲ally w̲here '∧ s̲\ p = 0◆ (n m̲in (t̲ally p) − 1) d̲rop p }
23
24⍝# a p_lus b: the sum.
25ˡp̲lus ← { a b →
26  n ← (t̲ally a) m̲ax t̲ally b
27  ʰt̲rim ((n̲eg n) t̲ake f̲loat a) + (n̲eg n) t̲ake f̲loat b
28}
29
30⍝# a t_imes b: the product: every product of a coefficient of a with
31⍝# one of b, summed by the power it belongs to.
32ˡt̲imes ← { a b →
33  m ← (f̲loat a) '× t̲able f̲loat b
34  k ← r̲avel (o̲ffsets t̲ally a) '+ t̲able o̲ffsets t̲ally b
35  v ← r̲avel m
36  ʰt̲rim '{ i → '+ r̲/ (k = i) r̲eplicate v } e̲ach o̲ffsets (t̲ally a) + (t̲ally b) − 1
37}
38
39⍝# d_erivative p: the derivative (0 for a constant).
40ˡd̲erivative ← { p →
41  1 = t̲ally p ? 1 r̲eshape 0.0
42  -1 d̲rop (f̲loat p) × f̲loat r̲ev o̲ffsets t̲ally p
43}
44
45⍝# i_ntegral p: the integral that is 0 at 0.
46ˡi̲ntegral ← { p → ((f̲loat p) ÷ f̲loat 1 + r̲ev o̲ffsets t̲ally p) c̲at 0.0 }
47
48⍝## Roots
49
50⍝# One Newton step for every x at once (a flat spot is nudged, not
51⍝# divided by).
52ʰs̲tep ← { p x →
53  d ← (ˡd̲erivative p) ˡa̲t x
54  x − (p ˡa̲t x) ÷ d + 0.000000000001 × f̲loat d = 0.0
55}
56
57⍝# r_oots p: the real roots, found by Newton's method from 64 starting
58⍝# points spread over the interval holding them all (Cauchy's bound),
59⍝# kept where the value is 0 to within 1e-9, rounded to 9 decimals.
60ˡr̲oots ← { p →
61  q ← ʰt̲rim f̲loat p
62  1 = t̲ally q ? 0 r̲eshape 0.0
63  b ← 1.0 + 'm̲ax r̲/ a̲bs (1 d̲rop q) ÷ f̲irst q
64  x ← 60 '{ y → q ʰs̲tep y } p̲ower b × ((f̲loat o̲ffsets 64) ÷ 31.5) − 1.0
65  x ← (0.000000001 > a̲bs q ˡa̲t x) r̲eplicate x
66  s̲ort u̲nique (f̲loat f̲loor 0.5 + x × 1000000000.0) ÷ 1000000000.0
67}
68
69⍝## Text
70
71⍝# A coefficient's size as text: whole numbers without a point.
72ʰs̲ize ← { c → a ← a̲bs c◆ a = f̲loat f̲loor a ? f̲ormat f̲loor a◆ f̲ormat a }
73
74⍝# x to the power k, as text.
75ʰp̲ow ← { k → k = 0 ? ""◆ k = 1 ? "x"◆ "x^" c̲at f̲ormat k }
76
77⍝# One term, without its sign: the size (left out when it is 1 and a
78⍝# power of x follows) and the power.
79ʰt̲erm ← { c k → (k > 0) ∧ 1.0 = a̲bs c ? ʰp̲ow k◆ (ʰs̲ize c) c̲at ʰp̲ow k }
80
81⍝# The sign between terms.
82ʰs̲ep ← { c → c < 0.0 ? " - "◆ " + " }
83
84⍝# t_ext p: the polynomial as text, 3x^2 - 2x + 1.
85ˡt̲ext ← { p →
86  c0 ← f̲loat ʰt̲rim p
87  keep ← (c0 ≠ 0.0) ∨ 1 = t̲ally c0
88  c ← keep r̲eplicate c0
89  k ← keep r̲eplicate r̲ev o̲ffsets t̲ally c0
90  first ← (((f̲irst c) < 0.0) r̲eplicate "-") c̲at (f̲irst c) ʰt̲erm f̲irst k
91  rest ← '{ i → (ʰs̲ep f̲irst i s̲elect c) c̲at (f̲irst i s̲elect c) ʰt̲erm f̲irst i s̲elect k } m̲ap 1 d̲rop r̲ange t̲ally c
92  first c̲at "" ᵗj̲oin rest
93}