I've found in 'addons/math/misc/brent.ijs'
implementation of Brent's method as an adverb. I would like to build a Newton's method as an adverb too but it's much harder than building tacit verbs.
Here is a explicit version of Newton's iteration:
newton_i =: 1 : '] - u % u d.1'
With such usage:
2&o. newton_i^:_ (1) NB. (-: 1p1) must be found
1.5708
2 o. 1.5708 NB. after substitution we get almost 0
_3.67321e_6
And of course, for convenience:
newton =: 1 : 'u newton_i^:_'
What's a tacit equivalent?
Per the comments, a short answer; the tacit equivalent to the original, explicit newton_i
and newton
are, respectively:
n_i =: d.0 1 (%/@:) (]`-`) (`:6)
newton =: n_i (^:_)
Some techniques for how such translations are obtained, in general, can be found on the J Forums.
The key insights here are that (a) that a function is identical to its own "zeroeth derivative", and that (b) we can calculate the "zeroeth" and first derivative of a function in J simultaneously, thanks to the language's array-oriented nature. The rest is mere stamp-collecting.
In an ideal world, given a function f
, we'd like to produce a verb train like (] - f % f d. 1)
. The problem is that tacit adverbial programming in J constrains us to produce a verb which mentions the input function (f
) once and only once.
So, instead, we use a sneaky trick: we calculate two derivatives of f
at the same time: the "zeroth" derivative (which is an identity function) and the first derivative.
load 'trig'
sin NB. Sine function (special case of the "circle functions", o.)
1&o.
sin d. 1 f. NB. First derivative of sine, sin'.
2&o.
sin d. 0 f. NB. "Zeroeth" derivative of sine, i.e. sine.
1&o."0
sin d. 0 1 f. NB. Both, resulting in two outputs.
(1&o. , 2&o.)"0
znfd =: d. 0 1 NB. Packaged up as a re-usable name.
sin znfd f.
(1&o. , 2&o.)"0
Then we simply insert a division between them:
dh =: znfd (%/@) NB. Quotient of first-derivative over 0th-derivattive
sin dh
%/@(sin d.0 1)
sin dh f.
%/@((1&o. , 2&o.)"0)
sin dh 1p1 NB. 0
_1.22465e_16
sin 1p1 NB. sin(pi) = 0
1.22465e_16
sin d. 1 ] 1p1 NB. sin'(pi) = -1
_1
sin dh 1p1 NB. sin(pi)/sin'(pi) = 0/-1 = 0
_1.22465e_16
The (%/@)
comes to the right of the znfd
because tacit adverbial programming in J is LIFO (i.e. left-to-right, where as "normal" J is right-to-left).
As I said, the remaining code is mere stamp collecting, using the standard tools to construct a verb-train which subtracts this quotient from the original input:
ssub =: (]`-`) (`:6) NB. x - f(x)
+: ssub NB. x - double(x)
] - +:
-: ssub NB. x - halve(x)
] - -:
-: ssub 16 NB. 16 - halve(16)
8
+: ssub 16 NB. 16 - double(16)
_16
*: ssub 16 NB. 16 - square(16)
_240
%: ssub 16 NB. 16 - sqrt(16)
12
Thus:
n_i =: znfd ssub NB. x - f'(x)/f(x)
And, finally, using "apply until fixed point" feature of ^:_
, we have:
newton =: n_i (^:_)
Voila.