I've updated the rosettacode page to include this implementation. And,
I plan on updating the j wiki polynomial division essay, also, once I
think through a bit how to describe this approach. (Descriptive
phrases and sentences would be welcome. Documentation needs good
perspectives, above all else.)
For the purpose of the rosettacode page, I used a different phrasing of rou:
rou =: {{(,j.) r,(%:0j1),j.+|.}.r=.^o.j.(i.@%&8 % -:) y}} NB.
roots of unity
This seems (to me) to be more concise, and straightforward, than the
purely tacit expression.
(That said, I'm not sure if this change in notation would be a good
fit for the FFT essay.)
Thanks,
--
Raul
On Tue, Mar 1, 2022 at 11:19 AM 'Michael Day' via Programming
<[email protected]> wrote:
>
> Raul Miller's and RE Boss's (revised) entries here refer:
> http://rosettacode.org/wiki/Cyclotomic_polynomial#J
>
> I think I'm partly responsible for pDiv in this (otherwise!) nicely
> concise and elegant
> presentation of J routines to derive these polynomials.
> Typically,
> NB. (1+5x+10x^2+10x^3+5x^4+x^5 )%(1+3x+3x^2+x^3) = 1+2x+x^2
> 1 5 10 10 5 1 pDiv 1 3 3 1
> 1 2 1
>
> pDiv is essentially elementary long division realised for polynomials
> rather than numbers;
> it's rather inefficient for high-order polynomials, and I was wondering
> how to speed it up.
> I found a paper by Madhu Sudan which addresses fast polynomial division:
> http://people.csail.mit.edu/madhu/ST12/scribe/lect06.pdf
> but couldn't understand his hints at an algorithm. I _think_ it
> successively improves
> approximations to the solution, working up from a constant term, then
> modulo x^2,
> then x^4 , but it's very muddled to my mind, and my mind couldn't see
> how to set about it!
>
> Then I realised a compromise is to exploit polynomial multiplication
> using Fast Fourier Transforms,
> as described first for this purpose, I believe, by Roger Hui; Henry's
> essay is to be found at
> https://code.jsoftware.com/wiki/Essays/FFT
>
> That article desribes polynomial multiplication, where we multiply the
> ffts of two
> polynomials extended to an appropriate length, a suitable power of 2.
> The inverse fft
> is then applied to yield the required polynomial product.
> If instead we divide one fft by another, and take the inverse fft, we
> get the polynomial
> quotient!
>
> So here's a crib of Raul's functions; I've merged cyclotomic and
> cyclotomic000, dispensing
> with "cache" and using the memoisation afforded by M. :
> (OK in fixed width font)
>
> ctfft=: {{ assert.0<y
> if. y = 1 do. _1 1 return. end.
> 'q p'=. __ q: y
> if. 1=#q do.
> ,(y%*/q) {."0 q#1
> elseif.2={.q do.
> ,(y%*/q) {."0 (* 1 _1 $~ #) ct */}.q
> elseif. 1 e. 1 < p do.
> ,(y%*/q) {."0 ct */q
> else.
> lgl =. {:$ ctlist =. ct "0 }:*/@>,{1,each q NB. ctlist is
> 2-d table of polynomial divisors
> lgd =. # dividend =. _1,(-y){.1 NB. (x^n) - 1,
> and its size
> lg =. >.&.(2&^.) lgl >. lgd NB. required
> lengths of all polynomials for fft transforms
> NB. really, "divisor" is the fft of the divisor!
> divisor =. */ fft"1 lg{."1 ctlist NB. FFT article
> doesn't deal with lists of multiplicands
> unpad roundreal ifft"1 divisor %~ fft lg{.dividend NB. similar to
> article's multiplication
> end.
> }} M.
>
> NB. fft, roundreal, ifft from FFT article.
>
> NB. Performance - poor for lower orders...
> NB. taskorderu is adverbial form of taskorder with required verb as left
> argument.
>
> timer'cyclotomic taskorderu 5' NB. timer returns elpased time,
> before cache/memoisation established!
> ┌────────┬───────────────────┐
> │1.163002│1 105 385 1365 1785│
> └────────┴───────────────────┘
> timer'ctfft taskorderu 5' NB. not too good here cf cyclotomic!
> ┌─────────┬───────────────────┐
> │3.5719986│1 105 385 1365 1785│
> └─────────┴───────────────────┘
> NB. continuing with cache/memoisation as above
> timer'cyclotomic taskorderu 10'
> ┌───────────┬─────────────────────────────────────────────┐
> │5 7.8180008│1 105 385 1365 1785 2805 3135 6545 6545 10465│
> └───────────┴─────────────────────────────────────────────┘
> timer'ctfft taskorderu 10' NB. some improvement!
> ┌────────┬─────────────────────────────────────────────┐
> │3 52.945│1 105 385 1365 1785 2805 3135 6545 6545 10465│
> └────────┴─────────────────────────────────────────────┘
>
> So - not an order of magnitude effect, but of some interest, I hope!
>
> Mike
>
> On 20/02/2022 18:32, Raul Miller wrote:
> > Oh, bleah... of course.
> >
> > y does not change. So that should be a test on {.x -- and that reveals
> > that I was using the wrong values throughout that section of code.
> >
> > Here's a fixed version:
> >
> > pDiv=: {{
> > q=. $j=. 2 + x -&# y
> > 'x y'=. x,:y
> > while. j=. j-1 do.
> > if. 0={.x do. j=. j-<:i=. 0 i.~ 0=x
> > q=. q,i#0
> > x=. i |.!.0 x
> > else.
> > q=. q, r=. x %&{. y
> > x=. 1 |.!.0 x - y*r
> > end.
> > end.q
> > }}
> >
> > Thanks,
> >
>
>
> --
> This email has been checked for viruses by Avast antivirus software.
> https://www.avast.com/antivirus
>
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
