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
