[Jprogramming] more fun with parallelism

[Henry Rich]

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