Sorry if I misattributed. The one time I met Roger face to face was in a London pub when he explained the idea to a few APLers, at the time he was collaborating with John Scholes, I think.
I haven’t explored the remainder issue, since divisibility is assured in the cyclotomic routines. Yes, I’d seen the complex stuff, but stuck with the pedestrian approach for starters. Given that “ctlist” may have several polynomials, there might be a better way of achieving fft “1 . Cheers, Mike Sent from my iPad > On 1 Mar 2022, at 16:31, Henry Rich <henryhr...@gmail.com> wrote: > > I supplied the idea for FFT multiplication (which was not original) and Roger > wrote the essay. > > Your idea for division is inspired. I have never seen division done this > way. I wonder if it works only when there is no remainder? If it works, it > will be fast for large polynomials. > > If speed matters, look at the fftw addon. If you have multiple FFTs to do > simultaneously, look into combining two real FFTs into one complex FFT; or > folding one real FFT into a half-size complex FFT. > > Henry Rich > > > >> On 3/1/2022 11:18 AM, 'Michael Day' via Programming 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 AVG. > https://www.avg.com > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
