f/\ is quadratic but f/\. shouldn't be.

Henry Rich

On 1/9/2023 11:20 AM, Marshall Lochbaum wrote:
It's not just Fibonacci numbers. This is equivalent to a general method
for computing continued fractions convergents in linear time, using
matrix multiplication. There's a nice explanation of why it works in
here:

https://perl.plover.com/classes/cftalk/INFO/gosper.txt

And a J version, with phi's continued fraction 1 1 1... :

    %&:{./"_1 +/ .*/\ (1 0,:~,&1)"0 ]10$1
1 2 1.5 1.66667 1.6 1.625 1.61538 1.61905 1.61765 1.61818

Here's the continued fraction sequence for e:

    2,1}.,1 1&,"0]2*1+i.3
2 1 2 1 1 4 1 1 6

    %&:{./"_1 +/ .*/\ (1 0,:~,&1)"0 ]2,1}.,1 1&,"0]2*1+i.3
2 3 2.66667 2.75 2.71429 2.71875 2.71795 2.71831 2.71828

My timings show J taking quadratic time for the scan, so this
formulation is pretty slow but shows the principle.

Marshall

On Sun, Jan 08, 2023 at 07:05:46PM -0800, Elijah Stone wrote:
Can we generate phi, the golden ratio, in parallel?

Of course we can!  Follows is an exposition of a classic method for doing
it, which may be of interest; I have not seen it satisfactorily described
elsewhere.

The classic method for generating phi in j uses a continued fraction:

    (+%)/10#1
1.61818

(It can be equivalently spelt (1+%)^:n]1.)

Using rational numbers clarifies slightly:

    (+%)/10#1x
89r55

Unsurprisingly, it's generating ratios of successive fibonacci numbers.  Can
it be parallelised?  The operation used for reduction is not associative:

    ((1 +% 1) +% 1) -: (1 +% (1 +% 1))
0

So it's not obvious how we would parallelise it.  Taking a step back, it
performs repeated division, where we don't know the divisor a priori,
instead generating it recursively, so attacking it thus seems hopeless.

Here's another method, based more directly on fibonacci numbers, which is
more promising:

    (u=. {:,+/)^:9]1 1
55 89

Here we generate fibonacci numbers via a recurrence relation, using two
successive terms to generate the next.

At first glance, this seems just as hopelessly sequential as the first
solution, but the use of ^: is a tell.  u^:n is the same as u@:u@:u@: ... n
times.  And _@:_ is associative!  So if we can somehow express the operation
of u as 'data', and ditto the composition of any number of us, then we can
create a big array of n copies of u, and then reduce @: over it (in
parallel), and finally apply the reduced u^:n to our y.

Of course, this all depends on finding an efficient way of representing
compositions of u.  If we can't do that, then we'll waste a bunch of
parallel work constructing compositions, only to _still_ need linear
sequential work at the end, once we apply our composed u.

Let's look at the operation of u _symbolically_:

u a,b is b,a+b
u u a,b is (a+b),(a+2*b)
u u u a,b is (a+2*b),((2*a)+(3*b))

In other words, u gives two results, each of which is a linear combination
of its two inputs.  And we have a convenient way of representing such
functions: matrices!  Hence obtains an alternate implementation of u, as a
matrix product:

    u2=. (0 1,:1 1)&(+/ . *)
    u2^:9]1 1 NB.same result as u
55 89

Hence, we can generate an efficient representation for u^:n in parallel,
with only logarithmic span, for matrix multiplication is isomorphic to
composition (and it is associative).  The resultant function always takes
the form of a 2x2 matrix, so we have only a constant amount of additional
work to do at the end.

    (+/ .*/9#,:0 1,:1 1) +/ .* 1 1
55 89

The astute may notice that this is wasted parallelism.  We use a parallel
reduction to find +/ .*/n#,:0 1,:1 1 in logarithmic time, but this is in
fact just the nth power of matrix 0 1,:1 1, which can be calculated
_sequentially_ in logarithmic time, by repeated squaring.  But this method
has a key advantage over that: if we want the entire fibonacci sequence, we
can generate it, still with logarithmic span and linear work, by simply
replacing our reduction with a scan:

    (+/ .*/\9#,:0 1,:1 1) +/ .* 1 1
  1  2
  2  3
  3  5
  5  8
  8 13
13 21
21 34
34 55
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm

----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm

Reply via email to