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
