1implement FFTs; 2include "sys.m"; 3 sys: Sys; 4 print: import sys; 5include "math.m"; 6 math: Math; 7 cos, sin, Degree, Pi: import math; 8include "ffts.m"; 9 10# by r. c. singleton, stanford research institute, sept. 1968 11# translated to limbo by eric grosse, jan 1997 12# arrays at(maxf), ck(maxf), bt(maxf), sk(maxf), and np(maxp) 13# are used for temporary storage. if the available storage 14# is insufficient, the program exits. 15# maxf must be >= the maximum prime factor of n. 16# maxp must be > the number of prime factors of n. 17# in addition, if the square-free portion k of n has two or 18# more prime factors, then maxp must be >= k-1. 19# array storage in nfac for a maximum of 15 prime factors of n. 20# if n has more than one square-free factor, the product of the 21# square-free factors must be <= 210 22 23ffts(a,b:array of real, ntot,n,nspan,isn:int){ 24 maxp: con 209; 25 i,ii,inc,j,jc,jf,jj,k,k1,k2,k3,k4,kk:int; 26 ks,kspan,kspnn,kt,m,maxf,nn,nt:int; 27 aa,aj,ajm,ajp,ak,akm,akp,bb,bj,bjm,bjp,bk,bkm,bkp:real; 28 c1,c2,c3,c72,cd,rad,radf,s1,s2,s3,s72,s120,sd:real; 29 maxf = 23; 30 if(math == nil){ 31 sys = load Sys Sys->PATH; 32 math = load Math Math->PATH; 33 } 34 nfac := array[12] of int; 35 np := array[maxp] of int; 36 at := array[23] of real; 37 ck := array[23] of real; 38 bt := array[23] of real; 39 sk := array[23] of real; 40 41 if(n<2) return; 42 inc = isn; 43 c72 = cos(72.*Degree); 44 s72 = sin(72.*Degree); 45 s120 = sin(120.*Degree); 46 rad = 2.*Pi; 47 if(isn<0){ 48 s72 = -s72; 49 s120 = -s120; 50 rad = -rad; 51 inc = -inc; 52 } 53 nt = inc*ntot; 54 ks = inc*nspan; 55 kspan = ks; 56 nn = nt-inc; 57 jc = ks/n; 58 radf = rad*real(jc)*0.5; 59 i = 0; 60 jf = 0; 61 62 # determine the factors of n 63 m = 0; 64 k = n; 65 while(k==k/16*16){ 66 m = m+1; 67 nfac[m] = 4; 68 k = k/16; 69 } 70 j = 3; 71 jj = 9; 72 for(;;) 73 if(k%jj==0){ 74 m = m+1; 75 nfac[m] = j; 76 k = k/jj; 77 }else{ 78 j = j+2; 79 jj = j*j; 80 if(jj>k) 81 break; 82 } 83 if(k<=4){ 84 kt = m; 85 nfac[m+1] = k; 86 if(k!=1) 87 m = m+1; 88 }else{ 89 if(k==k/4*4){ 90 m = m+1; 91 nfac[m] = 2; 92 k = k/4; 93 } 94 kt = m; 95 j = 2; 96 do{ 97 if(k%j==0){ 98 m = m+1; 99 nfac[m] = j; 100 k = k/j; 101 } 102 j = ((j+1)/2)*2+1; 103 }while(j<=k); 104 } 105 if(kt!=0){ 106 j = kt; 107 do{ 108 m = m+1; 109 nfac[m] = nfac[j]; 110 j = j-1; 111 }while(j!=0); 112 } 113 114 for(;;){ # compute fourier transform 115 sd = radf/real(kspan); 116 cd = sin(sd); 117 cd = 2.0*cd*cd; 118 sd = sin(sd+sd); 119 kk = 1; 120 i = i+1; 121 if(nfac[i]==2){ # transform for factor of 2 (including rotation factor) 122 kspan = kspan/2; 123 k1 = kspan+2; 124 for(;;){ 125 k2 = kk+kspan; 126 ak = a[k2-1]; 127 bk = b[k2-1]; 128 a[k2-1] = a[kk-1]-ak; 129 b[k2-1] = b[kk-1]-bk; 130 a[kk-1] = a[kk-1]+ak; 131 b[kk-1] = b[kk-1]+bk; 132 kk = k2+kspan; 133 if(kk>nn){ 134 kk = kk-nn; 135 if(kk>jc) 136 break; 137 } 138 } 139 if(kk>kspan) 140 break; 141 do{ 142 c1 = 1.0-cd; 143 s1 = sd; 144 for(;;){ 145 k2 = kk+kspan; 146 ak = a[kk-1]-a[k2-1]; 147 bk = b[kk-1]-b[k2-1]; 148 a[kk-1] = a[kk-1]+a[k2-1]; 149 b[kk-1] = b[kk-1]+b[k2-1]; 150 a[k2-1] = c1*ak-s1*bk; 151 b[k2-1] = s1*ak+c1*bk; 152 kk = k2+kspan; 153 if(kk>=nt){ 154 k2 = kk-nt; 155 c1 = -c1; 156 kk = k1-k2; 157 if(kk<=k2){ 158 ak = c1-(cd*c1+sd*s1); 159 s1 = (sd*c1-cd*s1)+s1; 160 c1 = 2.0-(ak*ak+s1*s1); 161 s1 = c1*s1; 162 c1 = c1*ak; 163 kk = kk+jc; 164 if(kk>=k2) 165 break; 166 } 167 } 168 } 169 k1 = k1+inc+inc; 170 kk = (k1-kspan)/2+jc; 171 }while(kk<=jc+jc); 172 }else{ # transform for factor of 4 173 if(nfac[i]!=4){ 174 # transform for odd factors 175 k = nfac[i]; 176 kspnn = kspan; 177 kspan = kspan/k; 178 if(k==3) 179 for(;;){ 180 # transform for factor of 3 (optional code) 181 k1 = kk+kspan; 182 k2 = k1+kspan; 183 ak = a[kk-1]; 184 bk = b[kk-1]; 185 aj = a[k1-1]+a[k2-1]; 186 bj = b[k1-1]+b[k2-1]; 187 a[kk-1] = ak+aj; 188 b[kk-1] = bk+bj; 189 ak = -0.5*aj+ak; 190 bk = -0.5*bj+bk; 191 aj = (a[k1-1]-a[k2-1])*s120; 192 bj = (b[k1-1]-b[k2-1])*s120; 193 a[k1-1] = ak-bj; 194 b[k1-1] = bk+aj; 195 a[k2-1] = ak+bj; 196 b[k2-1] = bk-aj; 197 kk = k2+kspan; 198 if(kk>=nn){ 199 kk = kk-nn; 200 if(kk>kspan) 201 break; 202 } 203 } 204 else if(k==5){ 205 # transform for factor of 5 (optional code) 206 c2 = c72*c72-s72*s72; 207 s2 = 2.0*c72*s72; 208 for(;;){ 209 k1 = kk+kspan; 210 k2 = k1+kspan; 211 k3 = k2+kspan; 212 k4 = k3+kspan; 213 akp = a[k1-1]+a[k4-1]; 214 akm = a[k1-1]-a[k4-1]; 215 bkp = b[k1-1]+b[k4-1]; 216 bkm = b[k1-1]-b[k4-1]; 217 ajp = a[k2-1]+a[k3-1]; 218 ajm = a[k2-1]-a[k3-1]; 219 bjp = b[k2-1]+b[k3-1]; 220 bjm = b[k2-1]-b[k3-1]; 221 aa = a[kk-1]; 222 bb = b[kk-1]; 223 a[kk-1] = aa+akp+ajp; 224 b[kk-1] = bb+bkp+bjp; 225 ak = akp*c72+ajp*c2+aa; 226 bk = bkp*c72+bjp*c2+bb; 227 aj = akm*s72+ajm*s2; 228 bj = bkm*s72+bjm*s2; 229 a[k1-1] = ak-bj; 230 a[k4-1] = ak+bj; 231 b[k1-1] = bk+aj; 232 b[k4-1] = bk-aj; 233 ak = akp*c2+ajp*c72+aa; 234 bk = bkp*c2+bjp*c72+bb; 235 aj = akm*s2-ajm*s72; 236 bj = bkm*s2-bjm*s72; 237 a[k2-1] = ak-bj; 238 a[k3-1] = ak+bj; 239 b[k2-1] = bk+aj; 240 b[k3-1] = bk-aj; 241 kk = k4+kspan; 242 if(kk>=nn){ 243 kk = kk-nn; 244 if(kk>kspan) 245 break; 246 } 247 } 248 }else{ 249 if(k!=jf){ 250 jf = k; 251 s1 = rad/real(k); 252 c1 = cos(s1); 253 s1 = sin(s1); 254 if(jf>maxf){ 255 sys->fprint(sys->fildes(2),"too many primes for fft"); 256 exit; 257 } 258 ck[jf-1] = 1.0; 259 sk[jf-1] = 0.0; 260 j = 1; 261 do{ 262 ck[j-1] = ck[k-1]*c1+sk[k-1]*s1; 263 sk[j-1] = ck[k-1]*s1-sk[k-1]*c1; 264 k = k-1; 265 ck[k-1] = ck[j-1]; 266 sk[k-1] = -sk[j-1]; 267 j = j+1; 268 }while(j<k); 269 } 270 for(;;){ 271 k1 = kk; 272 k2 = kk+kspnn; 273 aa = a[kk-1]; 274 bb = b[kk-1]; 275 ak = aa; 276 bk = bb; 277 j = 1; 278 k1 = k1+kspan; 279 do{ 280 k2 = k2-kspan; 281 j = j+1; 282 at[j-1] = a[k1-1]+a[k2-1]; 283 ak = at[j-1]+ak; 284 bt[j-1] = b[k1-1]+b[k2-1]; 285 bk = bt[j-1]+bk; 286 j = j+1; 287 at[j-1] = a[k1-1]-a[k2-1]; 288 bt[j-1] = b[k1-1]-b[k2-1]; 289 k1 = k1+kspan; 290 }while(k1<k2); 291 a[kk-1] = ak; 292 b[kk-1] = bk; 293 k1 = kk; 294 k2 = kk+kspnn; 295 j = 1; 296 do{ 297 k1 = k1+kspan; 298 k2 = k2-kspan; 299 jj = j; 300 ak = aa; 301 bk = bb; 302 aj = 0.0; 303 bj = 0.0; 304 k = 1; 305 do{ 306 k = k+1; 307 ak = at[k-1]*ck[jj-1]+ak; 308 bk = bt[k-1]*ck[jj-1]+bk; 309 k = k+1; 310 aj = at[k-1]*sk[jj-1]+aj; 311 bj = bt[k-1]*sk[jj-1]+bj; 312 jj = jj+j; 313 if(jj>jf) 314 jj = jj-jf; 315 }while(k<jf); 316 k = jf-j; 317 a[k1-1] = ak-bj; 318 b[k1-1] = bk+aj; 319 a[k2-1] = ak+bj; 320 b[k2-1] = bk-aj; 321 j = j+1; 322 }while(j<k); 323 kk = kk+kspnn; 324 if(kk>nn){ 325 kk = kk-nn; 326 if(kk>kspan) 327 break; 328 } 329 } 330 } 331 # multiply by rotation factor (except for factors of 2 and 4) 332 if(i==m) 333 break; 334 kk = jc+1; 335 do{ 336 c2 = 1.0-cd; 337 s1 = sd; 338 do{ 339 c1 = c2; 340 s2 = s1; 341 kk = kk+kspan; 342 for(;;){ 343 ak = a[kk-1]; 344 a[kk-1] = c2*ak-s2*b[kk-1]; 345 b[kk-1] = s2*ak+c2*b[kk-1]; 346 kk = kk+kspnn; 347 if(kk>nt){ 348 ak = s1*s2; 349 s2 = s1*c2+c1*s2; 350 c2 = c1*c2-ak; 351 kk = kk-nt+kspan; 352 if(kk>kspnn) 353 break; 354 } 355 } 356 c2 = c1-(cd*c1+sd*s1); 357 s1 = s1+(sd*c1-cd*s1); 358 c1 = 2.0-(c2*c2+s1*s1); 359 s1 = c1*s1; 360 c2 = c1*c2; 361 kk = kk-kspnn+jc; 362 }while(kk<=kspan); 363 kk = kk-kspan+jc+inc; 364 }while(kk<=jc+jc); 365 }else{ 366 kspnn = kspan; 367 kspan = kspan/4; 368 do{ 369 c1 = 1.; 370 s1 = 0.; 371 for(;;){ 372 k1 = kk+kspan; 373 k2 = k1+kspan; 374 k3 = k2+kspan; 375 akp = a[kk-1]+a[k2-1]; 376 akm = a[kk-1]-a[k2-1]; 377 ajp = a[k1-1]+a[k3-1]; 378 ajm = a[k1-1]-a[k3-1]; 379 a[kk-1] = akp+ajp; 380 ajp = akp-ajp; 381 bkp = b[kk-1]+b[k2-1]; 382 bkm = b[kk-1]-b[k2-1]; 383 bjp = b[k1-1]+b[k3-1]; 384 bjm = b[k1-1]-b[k3-1]; 385 b[kk-1] = bkp+bjp; 386 bjp = bkp-bjp; 387 do10 := 0; 388 if(isn<0){ 389 akp = akm+bjm; 390 akm = akm-bjm; 391 bkp = bkm-ajm; 392 bkm = bkm+ajm; 393 if(s1!=0.) do10 = 1; 394 }else{ 395 akp = akm-bjm; 396 akm = akm+bjm; 397 bkp = bkm+ajm; 398 bkm = bkm-ajm; 399 if(s1!=0.) do10 = 1; 400 } 401 if(do10){ 402 a[k1-1] = akp*c1-bkp*s1; 403 b[k1-1] = akp*s1+bkp*c1; 404 a[k2-1] = ajp*c2-bjp*s2; 405 b[k2-1] = ajp*s2+bjp*c2; 406 a[k3-1] = akm*c3-bkm*s3; 407 b[k3-1] = akm*s3+bkm*c3; 408 kk = k3+kspan; 409 if(kk<=nt) 410 continue; 411 }else{ 412 a[k1-1] = akp; 413 b[k1-1] = bkp; 414 a[k2-1] = ajp; 415 b[k2-1] = bjp; 416 a[k3-1] = akm; 417 b[k3-1] = bkm; 418 kk = k3+kspan; 419 if(kk<=nt) 420 continue; 421 } 422 c2 = c1-(cd*c1+sd*s1); 423 s1 = (sd*c1-cd*s1)+s1; 424 c1 = 2.0-(c2*c2+s1*s1); 425 s1 = c1*s1; 426 c1 = c1*c2; 427 c2 = c1*c1-s1*s1; 428 s2 = 2.0*c1*s1; 429 c3 = c2*c1-s2*s1; 430 s3 = c2*s1+s2*c1; 431 kk = kk-nt+jc; 432 if(kk>kspan) 433 break; 434 } 435 kk = kk-kspan+inc; 436 }while(kk<=jc); 437 if(kspan==jc) 438 break; 439 } 440 } 441 } # end "compute fourier transform" 442 443 # permute the results to normal order---done in two stages 444 # permutation for square factors of n 445 np[0] = ks; 446 if(kt!=0){ 447 k = kt+kt+1; 448 if(m<k) 449 k = k-1; 450 j = 1; 451 np[k] = jc; 452 do{ 453 np[j] = np[j-1]/nfac[j]; 454 np[k-1] = np[k]*nfac[j]; 455 j = j+1; 456 k = k-1; 457 }while(j<k); 458 k3 = np[k]; 459 kspan = np[1]; 460 kk = jc+1; 461 k2 = kspan+1; 462 j = 1; 463 if(n!=ntot){ 464 for(;;){ 465 # permutation for multivariate transform 466 k = kk+jc; 467 do{ 468 ak = a[kk-1]; 469 a[kk-1] = a[k2-1]; 470 a[k2-1] = ak; 471 bk = b[kk-1]; 472 b[kk-1] = b[k2-1]; 473 b[k2-1] = bk; 474 kk = kk+inc; 475 k2 = k2+inc; 476 }while(kk<k); 477 kk = kk+ks-jc; 478 k2 = k2+ks-jc; 479 if(kk>=nt){ 480 k2 = k2-nt+kspan; 481 kk = kk-nt+jc; 482 if(k2>=ks) 483 permm: for(;;){ 484 k2 = k2-np[j-1]; 485 j = j+1; 486 k2 = np[j]+k2; 487 if(k2<=np[j-1]){ 488 j = 1; 489 do{ 490 if(kk<k2) 491 break permm; 492 kk = kk+jc; 493 k2 = kspan+k2; 494 }while(k2<ks); 495 if(kk>=ks) 496 break permm; 497 } 498 } 499 } 500 } 501 jc = k3; 502 }else{ 503 for(;;){ 504 # permutation for single-variate transform (optional code) 505 ak = a[kk-1]; 506 a[kk-1] = a[k2-1]; 507 a[k2-1] = ak; 508 bk = b[kk-1]; 509 b[kk-1] = b[k2-1]; 510 b[k2-1] = bk; 511 kk = kk+inc; 512 k2 = kspan+k2; 513 if(k2>=ks) 514 perms: for(;;){ 515 k2 = k2-np[j-1]; 516 j = j+1; 517 k2 = np[j]+k2; 518 if(k2<=np[j-1]){ 519 j = 1; 520 do{ 521 if(kk<k2) 522 break perms; 523 kk = kk+inc; 524 k2 = kspan+k2; 525 }while(k2<ks); 526 if(kk>=ks) 527 break perms; 528 } 529 } 530 } 531 jc = k3; 532 } 533 } 534 if(2*kt+1>=m) 535 return; 536 kspnn = np[kt]; 537 # permutation for square-free factors of n 538 j = m-kt; 539 nfac[j+1] = 1; 540 do{ 541 nfac[j] = nfac[j]*nfac[j+1]; 542 j = j-1; 543 }while(j!=kt); 544 kt = kt+1; 545 nn = nfac[kt]-1; 546 if(nn<=maxp){ 547 jj = 0; 548 j = 0; 549 for(;;){ 550 k2 = nfac[kt]; 551 k = kt+1; 552 kk = nfac[k]; 553 j = j+1; 554 if(j>nn) 555 break; 556 for(;;){ 557 jj = kk+jj; 558 if(jj<k2) 559 break; 560 jj = jj-k2; 561 k2 = kk; 562 k = k+1; 563 kk = nfac[k]; 564 } 565 np[j-1] = jj; 566 } 567 # determine the permutation cycles of length greater than 1 568 j = 0; 569 for(;;){ 570 j = j+1; 571 kk = np[j-1]; 572 if(kk>=0) 573 if(kk==j){ 574 np[j-1] = -j; 575 if(j==nn) 576 break; 577 }else{ 578 do{ 579 k = kk; 580 kk = np[k-1]; 581 np[k-1] = -kk; 582 }while(kk!=j); 583 k3 = kk; 584 } 585 } 586 maxf = inc*maxf; 587 for(;;){ 588 j = k3+1; 589 nt = nt-kspnn; 590 ii = nt-inc+1; 591 if(nt<0) 592 break; 593 for(;;){ 594 j = j-1; 595 if(np[j-1]>=0){ 596 jj = jc; 597 do{ 598 kspan = jj; 599 if(jj>maxf) 600 kspan = maxf; 601 jj = jj-kspan; 602 k = np[j-1]; 603 kk = jc*k+ii+jj; 604 k1 = kk+kspan; 605 k2 = 0; 606 do{ 607 k2 = k2+1; 608 at[k2-1] = a[k1-1]; 609 bt[k2-1] = b[k1-1]; 610 k1 = k1-inc; 611 }while(k1!=kk); 612 do{ 613 k1 = kk+kspan; 614 k2 = k1-jc*(k+np[k-1]); 615 k = -np[k-1]; 616 do{ 617 a[k1-1] = a[k2-1]; 618 b[k1-1] = b[k2-1]; 619 k1 = k1-inc; 620 k2 = k2-inc; 621 }while(k1!=kk); 622 kk = k2; 623 }while(k!=j); 624 k1 = kk+kspan; 625 k2 = 0; 626 do{ 627 k2 = k2+1; 628 a[k1-1] = at[k2-1]; 629 b[k1-1] = bt[k2-1]; 630 k1 = k1-inc; 631 }while(k1!=kk); 632 }while(jj!=0); 633 if(j==1) 634 break; 635 } 636 } 637 } 638 } 639} 640