Using my caching algorithm, and resetting the cache inside the
timespacex sentence (using [cache=:2{.cache on the right hand side), I
see that your fft polynomial division approach is slightly over twice
as fast as my slightly optimized version of your long division
approach for taskorder 10 (or taskorderu 10). There's a bit of
variance here (2.2 times faster or 2.5 times faster), and it's slow
enough that I haven't run a lot of tests.I am also wondering if the fft division mechanism here can be optimized to take advantage of the rather small coefficients in these polynomials. Anyways, lots of tangents to explore. Thanks, -- Raul On Tue, Mar 1, 2022 at 11:19 AM 'Michael Day' via Programming <programm...@jsoftware.com> 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
