Something went wrong. Try again.
Monorepo for Aesthetic.Computer aesthetic.computer
Something went wrong. Try again.
C++
1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653//-----------------------------------------------------------------------------// Copyright 2012 Masanori Morise// Author: mmorise [at] meiji.ac.jp (Masanori Morise)// Last update: 2021/02/15//// This file represents the functions about FFT (Fast Fourier Transform)// implemented by Mr. Ooura, and wrapper functions implemented by M. Morise.// We can use these wrapper functions as well as the FFTW functions.// Please see the FFTW web-page to show the usage of the wrapper functions.// Ooura FFT:// (Japanese) http://www.kurims.kyoto-u.ac.jp/~ooura/index-j.html// (English) http://www.kurims.kyoto-u.ac.jp/~ooura/index.html// FFTW:// (English) http://www.fftw.org/// 2012/08/24 by M. Morise//-----------------------------------------------------------------------------#include "world/fft.h"
#include <math.h>#include <stdlib.h>
void cdft(int n, int isgn, double *a, int *ip, double *w);void rdft(int n, int isgn, double *a, int *ip, double *w);
namespace {static void BackwardFFT(fft_plan p) { if (p.c_out == NULL) { // c2r p.input[0] = p.c_in[0][0]; p.input[1] = p.c_in[p.n / 2][0]; for (int i = 1; i < p.n / 2; ++i) { p.input[i * 2] = p.c_in[i][0]; p.input[i * 2 + 1] = -p.c_in[i][1]; } rdft(p.n, -1, p.input, p.ip, p.w); for (int i = 0; i < p.n; ++i) p.out[i] = p.input[i] * 2.0; } else { // c2c for (int i = 0; i < p.n; ++i) { p.input[i * 2] = p.c_in[i][0]; p.input[i * 2 + 1] = p.c_in[i][1]; } cdft(p.n * 2, -1, p.input, p.ip, p.w); for (int i = 0; i < p.n; ++i) { p.c_out[i][0] = p.input[i * 2]; p.c_out[i][1] = -p.input[i * 2 + 1]; } }}
static void ForwardFFT(fft_plan p) { if (p.c_in == NULL) { // r2c for (int i = 0; i < p.n; ++i) p.input[i] = p.in[i]; rdft(p.n, 1, p.input, p.ip, p.w); p.c_out[0][0] = p.input[0]; p.c_out[0][1] = 0.0; for (int i = 1; i < p.n / 2; ++i) { p.c_out[i][0] = p.input[i * 2]; p.c_out[i][1] = -p.input[i * 2 + 1]; } p.c_out[p.n / 2][0] = p.input[1]; p.c_out[p.n / 2][1] = 0.0; } else { // c2c for (int i = 0; i < p.n; ++i) { p.input[i * 2] = p.c_in[i][0]; p.input[i * 2 + 1] = p.c_in[i][1]; } cdft(p.n * 2, 1, p.input, p.ip, p.w); for (int i = 0; i < p.n; ++i) { p.c_out[i][0] = p.input[i * 2]; p.c_out[i][1] = -p.input[i * 2 + 1]; } }}
} // namespace
fft_plan fft_plan_dft_1d(int n, fft_complex *in, fft_complex *out, int sign, unsigned int flags) { void makewt(int nw, int *ip, double *w);
fft_plan output = {0}; output.n = n; output.in = NULL; output.c_in = in; output.out = NULL; output.c_out = out; output.sign = sign; output.flags = flags; output.input = new double[n * 2]; output.ip = new int[n]; output.w = new double[n * 5 / 4];
output.ip[0] = 0; makewt(output.n >> 1, output.ip, output.w); return output;}
fft_plan fft_plan_dft_c2r_1d(int n, fft_complex *in, double *out, unsigned int flags) { void makewt(int nw, int *ip, double *w); void makect(int nc, int *ip, double *c);
fft_plan output = {0}; output.n = n; output.in = NULL; output.c_in = in; output.out = out; output.c_out = NULL; output.sign = FFT_BACKWARD; output.flags = flags; output.input = new double[n]; output.ip = new int[n]; output.w = new double[n * 5 / 4];
output.ip[0] = 0; makewt(output.n >> 2, output.ip, output.w); makect(output.n >> 2, output.ip, output.w + (output.n >> 2)); return output;}
fft_plan fft_plan_dft_r2c_1d(int n, double *in, fft_complex *out, unsigned int flags) { void makewt(int nw, int *ip, double *w); void makect(int nc, int *ip, double *c);
fft_plan output = {0}; output.n = n; output.in = in; output.c_in = NULL; output.out = NULL; output.c_out = out; output.sign = FFT_FORWARD; output.flags = flags; output.input = new double[n]; output.ip = new int[n]; output.w = new double[n * 5 / 4];
output.ip[0] = 0; makewt(output.n >> 2, output.ip, output.w); makect(output.n >> 2, output.ip, output.w + (output.n >> 2)); return output;}
void fft_execute(fft_plan p) { if (p.sign == FFT_FORWARD) { ForwardFFT(p); } else { // ifft BackwardFFT(p); }}
void fft_destroy_plan(fft_plan p) { p.n = 0; p.in = NULL; p.c_in = NULL; p.out = NULL; p.c_out = NULL; p.sign = 0; p.flags = 0; delete[] p.input; delete[] p.ip; delete[] p.w;}
//-----------------------------------------------------------------------// The following functions are reffered by// http://www.kurims.kyoto-u.ac.jp/~ooura/index.html
void cdft(int n, int isgn, double *a, int *ip, double *w) { void cftfsub(int n, double *a, int *ip, int nw, double *w); void cftbsub(int n, double *a, int *ip, int nw, double *w); int nw;
nw = ip[0]; if (isgn >= 0) { cftfsub(n, a, ip, nw, w); } else { cftbsub(n, a, ip, nw, w); }}
void rdft(int n, int isgn, double *a, int *ip, double *w) { void cftfsub(int n, double *a, int *ip, int nw, double *w); void cftbsub(int n, double *a, int *ip, int nw, double *w); void rftfsub(int n, double *a, int nc, double *c); void rftbsub(int n, double *a, int nc, double *c); double xi;
int nw = ip[0]; int nc = ip[1];
if (isgn >= 0) { if (n > 4) { cftfsub(n, a, ip, nw, w); rftfsub(n, a, nc, w + nw); } else if (n == 4) { cftfsub(n, a, ip, nw, w); } xi = a[0] - a[1]; a[0] += a[1]; a[1] = xi; } else { a[1] = 0.5 * (a[0] - a[1]); a[0] -= a[1]; if (n > 4) { rftbsub(n, a, nc, w + nw); cftbsub(n, a, ip, nw, w); } else if (n == 4) { cftbsub(n, a, ip, nw, w); } }}
void makewt(int nw, int *ip, double *w) { void makeipt(int nw, int *ip); int j, nwh, nw0, nw1; double delta, wn4r, wk1r, wk1i, wk3r, wk3i;
ip[0] = nw; ip[1] = 1; if (nw > 2) { nwh = nw >> 1; delta = atan(1.0) / nwh; wn4r = cos(delta * nwh); w[0] = 1; w[1] = wn4r; if (nwh == 4) { w[2] = cos(delta * 2); w[3] = sin(delta * 2); } else if (nwh > 4) { makeipt(nw, ip); w[2] = 0.5 / cos(delta * 2); w[3] = 0.5 / cos(delta * 6); for (j = 4; j < nwh; j += 4) { w[j] = cos(delta * j); w[j + 1] = sin(delta * j); w[j + 2] = cos(3 * delta * j); w[j + 3] = -sin(3 * delta * j); } } nw0 = 0; while (nwh > 2) { nw1 = nw0 + nwh; nwh >>= 1; w[nw1] = 1; w[nw1 + 1] = wn4r; if (nwh == 4) { wk1r = w[nw0 + 4]; wk1i = w[nw0 + 5]; w[nw1 + 2] = wk1r; w[nw1 + 3] = wk1i; } else if (nwh > 4) { wk1r = w[nw0 + 4]; wk3r = w[nw0 + 6]; w[nw1 + 2] = 0.5 / wk1r; w[nw1 + 3] = 0.5 / wk3r; for (j = 4; j < nwh; j += 4) { wk1r = w[nw0 + 2 * j]; wk1i = w[nw0 + 2 * j + 1]; wk3r = w[nw0 + 2 * j + 2]; wk3i = w[nw0 + 2 * j + 3]; w[nw1 + j] = wk1r; w[nw1 + j + 1] = wk1i; w[nw1 + j + 2] = wk3r; w[nw1 + j + 3] = wk3i; } } nw0 = nw1; } }}
void makeipt(int nw, int *ip) { int j, l, m, m2, p, q;
ip[2] = 0; ip[3] = 16; m = 2; for (l = nw; l > 32; l >>= 2) { m2 = m << 1; q = m2 << 3; for (j = m; j < m2; j++) { p = ip[j] << 2; ip[m + j] = p; ip[m2 + j] = p + q; } m = m2; }}
void makect(int nc, int *ip, double *c) { int j, nch; double delta;
ip[1] = nc; if (nc > 1) { nch = nc >> 1; delta = atan(1.0) / nch; c[0] = cos(delta * nch); c[nch] = 0.5 * c[0]; for (j = 1; j < nch; j++) { c[j] = 0.5 * cos(delta * j); c[nc - j] = 0.5 * sin(delta * j); } }}
// -------- child routines --------
void cftfsub(int n, double *a, int *ip, int nw, double *w) { void bitrv2(int n, int *ip, double *a); void bitrv216(double *a); void bitrv208(double *a); void cftf1st(int n, double *a, double *w); void cftrec4(int n, double *a, int nw, double *w); void cftleaf(int n, int isplt, double *a, int nw, double *w); void cftfx41(int n, double *a, int nw, double *w); void cftf161(double *a, double *w); void cftf081(double *a, double *w); void cftf040(double *a); void cftx020(double *a);
if (n > 8) { if (n > 32) { cftf1st(n, a, &w[nw - (n >> 2)]); if (n > 512) { cftrec4(n, a, nw, w); } else if (n > 128) { cftleaf(n, 1, a, nw, w); } else { cftfx41(n, a, nw, w); } bitrv2(n, ip, a); } else if (n == 32) { cftf161(a, &w[nw - 8]); bitrv216(a); } else { cftf081(a, w); bitrv208(a); } } else if (n == 8) { cftf040(a); } else if (n == 4) { cftx020(a); }}
void cftbsub(int n, double *a, int *ip, int nw, double *w) { void bitrv2conj(int n, int *ip, double *a); void bitrv216neg(double *a); void bitrv208neg(double *a); void cftb1st(int n, double *a, double *w); void cftrec4(int n, double *a, int nw, double *w); void cftleaf(int n, int isplt, double *a, int nw, double *w); void cftfx41(int n, double *a, int nw, double *w); void cftf161(double *a, double *w); void cftf081(double *a, double *w); void cftb040(double *a); void cftx020(double *a);
if (n > 8) { if (n > 32) { cftb1st(n, a, &w[nw - (n >> 2)]); if (n > 512) { cftrec4(n, a, nw, w); } else if (n > 128) { cftleaf(n, 1, a, nw, w); } else { cftfx41(n, a, nw, w); } bitrv2conj(n, ip, a); } else if (n == 32) { cftf161(a, &w[nw - 8]); bitrv216neg(a); } else { cftf081(a, w); bitrv208neg(a); } } else if (n == 8) { cftb040(a); } else if (n == 4) { cftx020(a); }}
void bitrv2(int n, int *ip, double *a) { int j, j1, k, k1, l, m, nh, nm; double xr, xi, yr, yi;
m = 1; for (l = n >> 2; l > 8; l >>= 2) { m <<= 1; } nh = n >> 1; nm = 4 * m; if (l == 8) { for (k = 0; k < m; k++) { for (j = 0; j < k; j++) { j1 = 4 * j + 2 * ip[m + k]; k1 = 4 * k + 2 * ip[m + j]; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 -= nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nh; k1 += 2; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= 2 * nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 += nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= 2 * nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += 2; k1 += nh; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 -= nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nh; k1 -= 2; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= 2 * nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 += nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= 2 * nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; } k1 = 4 * k + 2 * ip[m + k]; j1 = k1 + 2; k1 += nh; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 -= nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= 2; k1 -= nh; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nh + 2; k1 += nh + 2; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nh - nm; k1 += 2 * nm - 2; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; } } else { for (k = 0; k < m; k++) { for (j = 0; j < k; j++) { j1 = 4 * j + ip[m + k]; k1 = 4 * k + ip[m + j]; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nh; k1 += 2; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += 2; k1 += nh; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nh; k1 -= 2; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; } k1 = 4 * k + ip[m + k]; j1 = k1 + 2; k1 += nh; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; } }}
void bitrv2conj(int n, int *ip, double *a) { int j, j1, k, k1, l, m, nh, nm; double xr, xi, yr, yi;
m = 1; for (l = n >> 2; l > 8; l >>= 2) { m <<= 1; } nh = n >> 1; nm = 4 * m; if (l == 8) { for (k = 0; k < m; k++) { for (j = 0; j < k; j++) { j1 = 4 * j + 2 * ip[m + k]; k1 = 4 * k + 2 * ip[m + j]; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 -= nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nh; k1 += 2; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= 2 * nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 += nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= 2 * nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += 2; k1 += nh; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 -= nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nh; k1 -= 2; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= 2 * nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 += nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= 2 * nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; } k1 = 4 * k + 2 * ip[m + k]; j1 = k1 + 2; k1 += nh; a[j1 - 1] = -a[j1 - 1]; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; a[k1 + 3] = -a[k1 + 3]; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 -= nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= 2; k1 -= nh; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nh + 2; k1 += nh + 2; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nh - nm; k1 += 2 * nm - 2; a[j1 - 1] = -a[j1 - 1]; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; a[k1 + 3] = -a[k1 + 3]; } } else { for (k = 0; k < m; k++) { for (j = 0; j < k; j++) { j1 = 4 * j + ip[m + k]; k1 = 4 * k + ip[m + j]; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nh; k1 += 2; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += 2; k1 += nh; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nh; k1 -= 2; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 -= nm; k1 -= nm; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; } k1 = 4 * k + ip[m + k]; j1 = k1 + 2; k1 += nh; a[j1 - 1] = -a[j1 - 1]; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; a[k1 + 3] = -a[k1 + 3]; j1 += nm; k1 += nm; a[j1 - 1] = -a[j1 - 1]; xr = a[j1]; xi = -a[j1 + 1]; yr = a[k1]; yi = -a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; a[k1 + 3] = -a[k1 + 3]; } }}
void bitrv216(double *a) { double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i, x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
x1r = a[2]; x1i = a[3]; x2r = a[4]; x2i = a[5]; x3r = a[6]; x3i = a[7]; x4r = a[8]; x4i = a[9]; x5r = a[10]; x5i = a[11]; x7r = a[14]; x7i = a[15]; x8r = a[16]; x8i = a[17]; x10r = a[20]; x10i = a[21]; x11r = a[22]; x11i = a[23]; x12r = a[24]; x12i = a[25]; x13r = a[26]; x13i = a[27]; x14r = a[28]; x14i = a[29]; a[2] = x8r; a[3] = x8i; a[4] = x4r; a[5] = x4i; a[6] = x12r; a[7] = x12i; a[8] = x2r; a[9] = x2i; a[10] = x10r; a[11] = x10i; a[14] = x14r; a[15] = x14i; a[16] = x1r; a[17] = x1i; a[20] = x5r; a[21] = x5i; a[22] = x13r; a[23] = x13i; a[24] = x3r; a[25] = x3i; a[26] = x11r; a[27] = x11i; a[28] = x7r; a[29] = x7i;}
void bitrv216neg(double *a) { double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i, x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i, x15r, x15i;
x1r = a[2]; x1i = a[3]; x2r = a[4]; x2i = a[5]; x3r = a[6]; x3i = a[7]; x4r = a[8]; x4i = a[9]; x5r = a[10]; x5i = a[11]; x6r = a[12]; x6i = a[13]; x7r = a[14]; x7i = a[15]; x8r = a[16]; x8i = a[17]; x9r = a[18]; x9i = a[19]; x10r = a[20]; x10i = a[21]; x11r = a[22]; x11i = a[23]; x12r = a[24]; x12i = a[25]; x13r = a[26]; x13i = a[27]; x14r = a[28]; x14i = a[29]; x15r = a[30]; x15i = a[31]; a[2] = x15r; a[3] = x15i; a[4] = x7r; a[5] = x7i; a[6] = x11r; a[7] = x11i; a[8] = x3r; a[9] = x3i; a[10] = x13r; a[11] = x13i; a[12] = x5r; a[13] = x5i; a[14] = x9r; a[15] = x9i; a[16] = x1r; a[17] = x1i; a[18] = x14r; a[19] = x14i; a[20] = x6r; a[21] = x6i; a[22] = x10r; a[23] = x10i; a[24] = x2r; a[25] = x2i; a[26] = x12r; a[27] = x12i; a[28] = x4r; a[29] = x4i; a[30] = x8r; a[31] = x8i;}
void bitrv208(double *a) { double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
x1r = a[2]; x1i = a[3]; x3r = a[6]; x3i = a[7]; x4r = a[8]; x4i = a[9]; x6r = a[12]; x6i = a[13]; a[2] = x4r; a[3] = x4i; a[6] = x6r; a[7] = x6i; a[8] = x1r; a[9] = x1i; a[12] = x3r; a[13] = x3i;}
void bitrv208neg(double *a) { double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, x5r, x5i, x6r, x6i, x7r, x7i;
x1r = a[2]; x1i = a[3]; x2r = a[4]; x2i = a[5]; x3r = a[6]; x3i = a[7]; x4r = a[8]; x4i = a[9]; x5r = a[10]; x5i = a[11]; x6r = a[12]; x6i = a[13]; x7r = a[14]; x7i = a[15]; a[2] = x7r; a[3] = x7i; a[4] = x3r; a[5] = x3i; a[6] = x5r; a[7] = x5i; a[8] = x1r; a[9] = x1i; a[10] = x6r; a[11] = x6i; a[12] = x2r; a[13] = x2i; a[14] = x4r; a[15] = x4i;}
void cftf1st(int n, double *a, double *w) { int j, j0, j1, j2, j3, k, m, mh; double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i; double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
mh = n >> 3; m = 2 * mh; j1 = m; j2 = j1 + m; j3 = j2 + m; x0r = a[0] + a[j2]; x0i = a[1] + a[j2 + 1]; x1r = a[0] - a[j2]; x1i = a[1] - a[j2 + 1]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; a[0] = x0r + x2r; a[1] = x0i + x2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i - x2i; a[j2] = x1r - x3i; a[j2 + 1] = x1i + x3r; a[j3] = x1r + x3i; a[j3 + 1] = x1i - x3r; wn4r = w[1]; csc1 = w[2]; csc3 = w[3]; wd1r = 1; wd1i = 0; wd3r = 1; wd3i = 0; k = 0; for (j = 2; j < mh - 2; j += 4) { k += 4; wk1r = csc1 * (wd1r + w[k]); wk1i = csc1 * (wd1i + w[k + 1]); wk3r = csc3 * (wd3r + w[k + 2]); wk3i = csc3 * (wd3i + w[k + 3]); wd1r = w[k]; wd1i = w[k + 1]; wd3r = w[k + 2]; wd3i = w[k + 3]; j1 = j + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j] + a[j2]; x0i = a[j + 1] + a[j2 + 1]; x1r = a[j] - a[j2]; x1i = a[j + 1] - a[j2 + 1]; y0r = a[j + 2] + a[j2 + 2]; y0i = a[j + 3] + a[j2 + 3]; y1r = a[j + 2] - a[j2 + 2]; y1i = a[j + 3] - a[j2 + 3]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; y2r = a[j1 + 2] + a[j3 + 2]; y2i = a[j1 + 3] + a[j3 + 3]; y3r = a[j1 + 2] - a[j3 + 2]; y3i = a[j1 + 3] - a[j3 + 3]; a[j] = x0r + x2r; a[j + 1] = x0i + x2i; a[j + 2] = y0r + y2r; a[j + 3] = y0i + y2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i - x2i; a[j1 + 2] = y0r - y2r; a[j1 + 3] = y0i - y2i; x0r = x1r - x3i; x0i = x1i + x3r; a[j2] = wk1r * x0r - wk1i * x0i; a[j2 + 1] = wk1r * x0i + wk1i * x0r; x0r = y1r - y3i; x0i = y1i + y3r; a[j2 + 2] = wd1r * x0r - wd1i * x0i; a[j2 + 3] = wd1r * x0i + wd1i * x0r; x0r = x1r + x3i; x0i = x1i - x3r; a[j3] = wk3r * x0r + wk3i * x0i; a[j3 + 1] = wk3r * x0i - wk3i * x0r; x0r = y1r + y3i; x0i = y1i - y3r; a[j3 + 2] = wd3r * x0r + wd3i * x0i; a[j3 + 3] = wd3r * x0i - wd3i * x0r; j0 = m - j; j1 = j0 + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j0] + a[j2]; x0i = a[j0 + 1] + a[j2 + 1]; x1r = a[j0] - a[j2]; x1i = a[j0 + 1] - a[j2 + 1]; y0r = a[j0 - 2] + a[j2 - 2]; y0i = a[j0 - 1] + a[j2 - 1]; y1r = a[j0 - 2] - a[j2 - 2]; y1i = a[j0 - 1] - a[j2 - 1]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; y2r = a[j1 - 2] + a[j3 - 2]; y2i = a[j1 - 1] + a[j3 - 1]; y3r = a[j1 - 2] - a[j3 - 2]; y3i = a[j1 - 1] - a[j3 - 1]; a[j0] = x0r + x2r; a[j0 + 1] = x0i + x2i; a[j0 - 2] = y0r + y2r; a[j0 - 1] = y0i + y2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i - x2i; a[j1 - 2] = y0r - y2r; a[j1 - 1] = y0i - y2i; x0r = x1r - x3i; x0i = x1i + x3r; a[j2] = wk1i * x0r - wk1r * x0i; a[j2 + 1] = wk1i * x0i + wk1r * x0r; x0r = y1r - y3i; x0i = y1i + y3r; a[j2 - 2] = wd1i * x0r - wd1r * x0i; a[j2 - 1] = wd1i * x0i + wd1r * x0r; x0r = x1r + x3i; x0i = x1i - x3r; a[j3] = wk3i * x0r + wk3r * x0i; a[j3 + 1] = wk3i * x0i - wk3r * x0r; x0r = y1r + y3i; x0i = y1i - y3r; a[j3 - 2] = wd3i * x0r + wd3r * x0i; a[j3 - 1] = wd3i * x0i - wd3r * x0r; } wk1r = csc1 * (wd1r + wn4r); wk1i = csc1 * (wd1i + wn4r); wk3r = csc3 * (wd3r - wn4r); wk3i = csc3 * (wd3i - wn4r); j0 = mh; j1 = j0 + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j0 - 2] + a[j2 - 2]; x0i = a[j0 - 1] + a[j2 - 1]; x1r = a[j0 - 2] - a[j2 - 2]; x1i = a[j0 - 1] - a[j2 - 1]; x2r = a[j1 - 2] + a[j3 - 2]; x2i = a[j1 - 1] + a[j3 - 1]; x3r = a[j1 - 2] - a[j3 - 2]; x3i = a[j1 - 1] - a[j3 - 1]; a[j0 - 2] = x0r + x2r; a[j0 - 1] = x0i + x2i; a[j1 - 2] = x0r - x2r; a[j1 - 1] = x0i - x2i; x0r = x1r - x3i; x0i = x1i + x3r; a[j2 - 2] = wk1r * x0r - wk1i * x0i; a[j2 - 1] = wk1r * x0i + wk1i * x0r; x0r = x1r + x3i; x0i = x1i - x3r; a[j3 - 2] = wk3r * x0r + wk3i * x0i; a[j3 - 1] = wk3r * x0i - wk3i * x0r; x0r = a[j0] + a[j2]; x0i = a[j0 + 1] + a[j2 + 1]; x1r = a[j0] - a[j2]; x1i = a[j0 + 1] - a[j2 + 1]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; a[j0] = x0r + x2r; a[j0 + 1] = x0i + x2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i - x2i; x0r = x1r - x3i; x0i = x1i + x3r; a[j2] = wn4r * (x0r - x0i); a[j2 + 1] = wn4r * (x0i + x0r); x0r = x1r + x3i; x0i = x1i - x3r; a[j3] = -wn4r * (x0r + x0i); a[j3 + 1] = -wn4r * (x0i - x0r); x0r = a[j0 + 2] + a[j2 + 2]; x0i = a[j0 + 3] + a[j2 + 3]; x1r = a[j0 + 2] - a[j2 + 2]; x1i = a[j0 + 3] - a[j2 + 3]; x2r = a[j1 + 2] + a[j3 + 2]; x2i = a[j1 + 3] + a[j3 + 3]; x3r = a[j1 + 2] - a[j3 + 2]; x3i = a[j1 + 3] - a[j3 + 3]; a[j0 + 2] = x0r + x2r; a[j0 + 3] = x0i + x2i; a[j1 + 2] = x0r - x2r; a[j1 + 3] = x0i - x2i; x0r = x1r - x3i; x0i = x1i + x3r; a[j2 + 2] = wk1i * x0r - wk1r * x0i; a[j2 + 3] = wk1i * x0i + wk1r * x0r; x0r = x1r + x3i; x0i = x1i - x3r; a[j3 + 2] = wk3i * x0r + wk3r * x0i; a[j3 + 3] = wk3i * x0i - wk3r * x0r;}
void cftb1st(int n, double *a, double *w) { int j, j0, j1, j2, j3, k, m, mh; double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i; double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
mh = n >> 3; m = 2 * mh; j1 = m; j2 = j1 + m; j3 = j2 + m; x0r = a[0] + a[j2]; x0i = -a[1] - a[j2 + 1]; x1r = a[0] - a[j2]; x1i = -a[1] + a[j2 + 1]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; a[0] = x0r + x2r; a[1] = x0i - x2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i + x2i; a[j2] = x1r + x3i; a[j2 + 1] = x1i + x3r; a[j3] = x1r - x3i; a[j3 + 1] = x1i - x3r; wn4r = w[1]; csc1 = w[2]; csc3 = w[3]; wd1r = 1; wd1i = 0; wd3r = 1; wd3i = 0; k = 0; for (j = 2; j < mh - 2; j += 4) { k += 4; wk1r = csc1 * (wd1r + w[k]); wk1i = csc1 * (wd1i + w[k + 1]); wk3r = csc3 * (wd3r + w[k + 2]); wk3i = csc3 * (wd3i + w[k + 3]); wd1r = w[k]; wd1i = w[k + 1]; wd3r = w[k + 2]; wd3i = w[k + 3]; j1 = j + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j] + a[j2]; x0i = -a[j + 1] - a[j2 + 1]; x1r = a[j] - a[j2]; x1i = -a[j + 1] + a[j2 + 1]; y0r = a[j + 2] + a[j2 + 2]; y0i = -a[j + 3] - a[j2 + 3]; y1r = a[j + 2] - a[j2 + 2]; y1i = -a[j + 3] + a[j2 + 3]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; y2r = a[j1 + 2] + a[j3 + 2]; y2i = a[j1 + 3] + a[j3 + 3]; y3r = a[j1 + 2] - a[j3 + 2]; y3i = a[j1 + 3] - a[j3 + 3]; a[j] = x0r + x2r; a[j + 1] = x0i - x2i; a[j + 2] = y0r + y2r; a[j + 3] = y0i - y2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i + x2i; a[j1 + 2] = y0r - y2r; a[j1 + 3] = y0i + y2i; x0r = x1r + x3i; x0i = x1i + x3r; a[j2] = wk1r * x0r - wk1i * x0i; a[j2 + 1] = wk1r * x0i + wk1i * x0r; x0r = y1r + y3i; x0i = y1i + y3r; a[j2 + 2] = wd1r * x0r - wd1i * x0i; a[j2 + 3] = wd1r * x0i + wd1i * x0r; x0r = x1r - x3i; x0i = x1i - x3r; a[j3] = wk3r * x0r + wk3i * x0i; a[j3 + 1] = wk3r * x0i - wk3i * x0r; x0r = y1r - y3i; x0i = y1i - y3r; a[j3 + 2] = wd3r * x0r + wd3i * x0i; a[j3 + 3] = wd3r * x0i - wd3i * x0r; j0 = m - j; j1 = j0 + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j0] + a[j2]; x0i = -a[j0 + 1] - a[j2 + 1]; x1r = a[j0] - a[j2]; x1i = -a[j0 + 1] + a[j2 + 1]; y0r = a[j0 - 2] + a[j2 - 2]; y0i = -a[j0 - 1] - a[j2 - 1]; y1r = a[j0 - 2] - a[j2 - 2]; y1i = -a[j0 - 1] + a[j2 - 1]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; y2r = a[j1 - 2] + a[j3 - 2]; y2i = a[j1 - 1] + a[j3 - 1]; y3r = a[j1 - 2] - a[j3 - 2]; y3i = a[j1 - 1] - a[j3 - 1]; a[j0] = x0r + x2r; a[j0 + 1] = x0i - x2i; a[j0 - 2] = y0r + y2r; a[j0 - 1] = y0i - y2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i + x2i; a[j1 - 2] = y0r - y2r; a[j1 - 1] = y0i + y2i; x0r = x1r + x3i; x0i = x1i + x3r; a[j2] = wk1i * x0r - wk1r * x0i; a[j2 + 1] = wk1i * x0i + wk1r * x0r; x0r = y1r + y3i; x0i = y1i + y3r; a[j2 - 2] = wd1i * x0r - wd1r * x0i; a[j2 - 1] = wd1i * x0i + wd1r * x0r; x0r = x1r - x3i; x0i = x1i - x3r; a[j3] = wk3i * x0r + wk3r * x0i; a[j3 + 1] = wk3i * x0i - wk3r * x0r; x0r = y1r - y3i; x0i = y1i - y3r; a[j3 - 2] = wd3i * x0r + wd3r * x0i; a[j3 - 1] = wd3i * x0i - wd3r * x0r; } wk1r = csc1 * (wd1r + wn4r); wk1i = csc1 * (wd1i + wn4r); wk3r = csc3 * (wd3r - wn4r); wk3i = csc3 * (wd3i - wn4r); j0 = mh; j1 = j0 + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j0 - 2] + a[j2 - 2]; x0i = -a[j0 - 1] - a[j2 - 1]; x1r = a[j0 - 2] - a[j2 - 2]; x1i = -a[j0 - 1] + a[j2 - 1]; x2r = a[j1 - 2] + a[j3 - 2]; x2i = a[j1 - 1] + a[j3 - 1]; x3r = a[j1 - 2] - a[j3 - 2]; x3i = a[j1 - 1] - a[j3 - 1]; a[j0 - 2] = x0r + x2r; a[j0 - 1] = x0i - x2i; a[j1 - 2] = x0r - x2r; a[j1 - 1] = x0i + x2i; x0r = x1r + x3i; x0i = x1i + x3r; a[j2 - 2] = wk1r * x0r - wk1i * x0i; a[j2 - 1] = wk1r * x0i + wk1i * x0r; x0r = x1r - x3i; x0i = x1i - x3r; a[j3 - 2] = wk3r * x0r + wk3i * x0i; a[j3 - 1] = wk3r * x0i - wk3i * x0r; x0r = a[j0] + a[j2]; x0i = -a[j0 + 1] - a[j2 + 1]; x1r = a[j0] - a[j2]; x1i = -a[j0 + 1] + a[j2 + 1]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; a[j0] = x0r + x2r; a[j0 + 1] = x0i - x2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i + x2i; x0r = x1r + x3i; x0i = x1i + x3r; a[j2] = wn4r * (x0r - x0i); a[j2 + 1] = wn4r * (x0i + x0r); x0r = x1r - x3i; x0i = x1i - x3r; a[j3] = -wn4r * (x0r + x0i); a[j3 + 1] = -wn4r * (x0i - x0r); x0r = a[j0 + 2] + a[j2 + 2]; x0i = -a[j0 + 3] - a[j2 + 3]; x1r = a[j0 + 2] - a[j2 + 2]; x1i = -a[j0 + 3] + a[j2 + 3]; x2r = a[j1 + 2] + a[j3 + 2]; x2i = a[j1 + 3] + a[j3 + 3]; x3r = a[j1 + 2] - a[j3 + 2]; x3i = a[j1 + 3] - a[j3 + 3]; a[j0 + 2] = x0r + x2r; a[j0 + 3] = x0i - x2i; a[j1 + 2] = x0r - x2r; a[j1 + 3] = x0i + x2i; x0r = x1r + x3i; x0i = x1i + x3r; a[j2 + 2] = wk1i * x0r - wk1r * x0i; a[j2 + 3] = wk1i * x0i + wk1r * x0r; x0r = x1r - x3i; x0i = x1i - x3r; a[j3 + 2] = wk3i * x0r + wk3r * x0i; a[j3 + 3] = wk3i * x0i - wk3r * x0r;}
void cftrec4(int n, double *a, int nw, double *w) { int cfttree(int n, int j, int k, double *a, int nw, double *w); void cftleaf(int n, int isplt, double *a, int nw, double *w); void cftmdl1(int n, double *a, double *w); int isplt, j, k, m;
m = n; while (m > 512) { m >>= 2; cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]); } cftleaf(m, 1, &a[n - m], nw, w); k = 0; for (j = n - m; j > 0; j -= m) { k++; isplt = cfttree(m, j, k, a, nw, w); cftleaf(m, isplt, &a[j - m], nw, w); }}
int cfttree(int n, int j, int k, double *a, int nw, double *w) { void cftmdl1(int n, double *a, double *w); void cftmdl2(int n, double *a, double *w); int i, isplt, m;
if ((k & 3) != 0) { isplt = k & 1; if (isplt != 0) { cftmdl1(n, &a[j - n], &w[nw - (n >> 1)]); } else { cftmdl2(n, &a[j - n], &w[nw - n]); } } else { m = n; for (i = k; (i & 3) == 0; i >>= 2) { m <<= 2; } isplt = i & 1; if (isplt != 0) { while (m > 128) { cftmdl1(m, &a[j - m], &w[nw - (m >> 1)]); m >>= 2; } } else { while (m > 128) { cftmdl2(m, &a[j - m], &w[nw - m]); m >>= 2; } } } return isplt;}
void cftleaf(int n, int isplt, double *a, int nw, double *w) { void cftmdl1(int n, double *a, double *w); void cftmdl2(int n, double *a, double *w); void cftf161(double *a, double *w); void cftf162(double *a, double *w); void cftf081(double *a, double *w); void cftf082(double *a, double *w);
if (n == 512) { cftmdl1(128, a, &w[nw - 64]); cftf161(a, &w[nw - 8]); cftf162(&a[32], &w[nw - 32]); cftf161(&a[64], &w[nw - 8]); cftf161(&a[96], &w[nw - 8]); cftmdl2(128, &a[128], &w[nw - 128]); cftf161(&a[128], &w[nw - 8]); cftf162(&a[160], &w[nw - 32]); cftf161(&a[192], &w[nw - 8]); cftf162(&a[224], &w[nw - 32]); cftmdl1(128, &a[256], &w[nw - 64]); cftf161(&a[256], &w[nw - 8]); cftf162(&a[288], &w[nw - 32]); cftf161(&a[320], &w[nw - 8]); cftf161(&a[352], &w[nw - 8]); if (isplt != 0) { cftmdl1(128, &a[384], &w[nw - 64]); cftf161(&a[480], &w[nw - 8]); } else { cftmdl2(128, &a[384], &w[nw - 128]); cftf162(&a[480], &w[nw - 32]); } cftf161(&a[384], &w[nw - 8]); cftf162(&a[416], &w[nw - 32]); cftf161(&a[448], &w[nw - 8]); } else { cftmdl1(64, a, &w[nw - 32]); cftf081(a, &w[nw - 8]); cftf082(&a[16], &w[nw - 8]); cftf081(&a[32], &w[nw - 8]); cftf081(&a[48], &w[nw - 8]); cftmdl2(64, &a[64], &w[nw - 64]); cftf081(&a[64], &w[nw - 8]); cftf082(&a[80], &w[nw - 8]); cftf081(&a[96], &w[nw - 8]); cftf082(&a[112], &w[nw - 8]); cftmdl1(64, &a[128], &w[nw - 32]); cftf081(&a[128], &w[nw - 8]); cftf082(&a[144], &w[nw - 8]); cftf081(&a[160], &w[nw - 8]); cftf081(&a[176], &w[nw - 8]); if (isplt != 0) { cftmdl1(64, &a[192], &w[nw - 32]); cftf081(&a[240], &w[nw - 8]); } else { cftmdl2(64, &a[192], &w[nw - 64]); cftf082(&a[240], &w[nw - 8]); } cftf081(&a[192], &w[nw - 8]); cftf082(&a[208], &w[nw - 8]); cftf081(&a[224], &w[nw - 8]); }}
void cftmdl1(int n, double *a, double *w) { int j, j0, j1, j2, j3, k, m, mh; double wn4r, wk1r, wk1i, wk3r, wk3i; double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
mh = n >> 3; m = 2 * mh; j1 = m; j2 = j1 + m; j3 = j2 + m; x0r = a[0] + a[j2]; x0i = a[1] + a[j2 + 1]; x1r = a[0] - a[j2]; x1i = a[1] - a[j2 + 1]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; a[0] = x0r + x2r; a[1] = x0i + x2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i - x2i; a[j2] = x1r - x3i; a[j2 + 1] = x1i + x3r; a[j3] = x1r + x3i; a[j3 + 1] = x1i - x3r; wn4r = w[1]; k = 0; for (j = 2; j < mh; j += 2) { k += 4; wk1r = w[k]; wk1i = w[k + 1]; wk3r = w[k + 2]; wk3i = w[k + 3]; j1 = j + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j] + a[j2]; x0i = a[j + 1] + a[j2 + 1]; x1r = a[j] - a[j2]; x1i = a[j + 1] - a[j2 + 1]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; a[j] = x0r + x2r; a[j + 1] = x0i + x2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i - x2i; x0r = x1r - x3i; x0i = x1i + x3r; a[j2] = wk1r * x0r - wk1i * x0i; a[j2 + 1] = wk1r * x0i + wk1i * x0r; x0r = x1r + x3i; x0i = x1i - x3r; a[j3] = wk3r * x0r + wk3i * x0i; a[j3 + 1] = wk3r * x0i - wk3i * x0r; j0 = m - j; j1 = j0 + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j0] + a[j2]; x0i = a[j0 + 1] + a[j2 + 1]; x1r = a[j0] - a[j2]; x1i = a[j0 + 1] - a[j2 + 1]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; a[j0] = x0r + x2r; a[j0 + 1] = x0i + x2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i - x2i; x0r = x1r - x3i; x0i = x1i + x3r; a[j2] = wk1i * x0r - wk1r * x0i; a[j2 + 1] = wk1i * x0i + wk1r * x0r; x0r = x1r + x3i; x0i = x1i - x3r; a[j3] = wk3i * x0r + wk3r * x0i; a[j3 + 1] = wk3i * x0i - wk3r * x0r; } j0 = mh; j1 = j0 + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j0] + a[j2]; x0i = a[j0 + 1] + a[j2 + 1]; x1r = a[j0] - a[j2]; x1i = a[j0 + 1] - a[j2 + 1]; x2r = a[j1] + a[j3]; x2i = a[j1 + 1] + a[j3 + 1]; x3r = a[j1] - a[j3]; x3i = a[j1 + 1] - a[j3 + 1]; a[j0] = x0r + x2r; a[j0 + 1] = x0i + x2i; a[j1] = x0r - x2r; a[j1 + 1] = x0i - x2i; x0r = x1r - x3i; x0i = x1i + x3r; a[j2] = wn4r * (x0r - x0i); a[j2 + 1] = wn4r * (x0i + x0r); x0r = x1r + x3i; x0i = x1i - x3r; a[j3] = -wn4r * (x0r + x0i); a[j3 + 1] = -wn4r * (x0i - x0r);}
void cftmdl2(int n, double *a, double *w) { int j, j0, j1, j2, j3, k, kr, m, mh; double wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i; double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
mh = n >> 3; m = 2 * mh; wn4r = w[1]; j1 = m; j2 = j1 + m; j3 = j2 + m; x0r = a[0] - a[j2 + 1]; x0i = a[1] + a[j2]; x1r = a[0] + a[j2 + 1]; x1i = a[1] - a[j2]; x2r = a[j1] - a[j3 + 1]; x2i = a[j1 + 1] + a[j3]; x3r = a[j1] + a[j3 + 1]; x3i = a[j1 + 1] - a[j3]; y0r = wn4r * (x2r - x2i); y0i = wn4r * (x2i + x2r); a[0] = x0r + y0r; a[1] = x0i + y0i; a[j1] = x0r - y0r; a[j1 + 1] = x0i - y0i; y0r = wn4r * (x3r - x3i); y0i = wn4r * (x3i + x3r); a[j2] = x1r - y0i; a[j2 + 1] = x1i + y0r; a[j3] = x1r + y0i; a[j3 + 1] = x1i - y0r; k = 0; kr = 2 * m; for (j = 2; j < mh; j += 2) { k += 4; wk1r = w[k]; wk1i = w[k + 1]; wk3r = w[k + 2]; wk3i = w[k + 3]; kr -= 4; wd1i = w[kr]; wd1r = w[kr + 1]; wd3i = w[kr + 2]; wd3r = w[kr + 3]; j1 = j + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j] - a[j2 + 1]; x0i = a[j + 1] + a[j2]; x1r = a[j] + a[j2 + 1]; x1i = a[j + 1] - a[j2]; x2r = a[j1] - a[j3 + 1]; x2i = a[j1 + 1] + a[j3]; x3r = a[j1] + a[j3 + 1]; x3i = a[j1 + 1] - a[j3]; y0r = wk1r * x0r - wk1i * x0i; y0i = wk1r * x0i + wk1i * x0r; y2r = wd1r * x2r - wd1i * x2i; y2i = wd1r * x2i + wd1i * x2r; a[j] = y0r + y2r; a[j + 1] = y0i + y2i; a[j1] = y0r - y2r; a[j1 + 1] = y0i - y2i; y0r = wk3r * x1r + wk3i * x1i; y0i = wk3r * x1i - wk3i * x1r; y2r = wd3r * x3r + wd3i * x3i; y2i = wd3r * x3i - wd3i * x3r; a[j2] = y0r + y2r; a[j2 + 1] = y0i + y2i; a[j3] = y0r - y2r; a[j3 + 1] = y0i - y2i; j0 = m - j; j1 = j0 + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j0] - a[j2 + 1]; x0i = a[j0 + 1] + a[j2]; x1r = a[j0] + a[j2 + 1]; x1i = a[j0 + 1] - a[j2]; x2r = a[j1] - a[j3 + 1]; x2i = a[j1 + 1] + a[j3]; x3r = a[j1] + a[j3 + 1]; x3i = a[j1 + 1] - a[j3]; y0r = wd1i * x0r - wd1r * x0i; y0i = wd1i * x0i + wd1r * x0r; y2r = wk1i * x2r - wk1r * x2i; y2i = wk1i * x2i + wk1r * x2r; a[j0] = y0r + y2r; a[j0 + 1] = y0i + y2i; a[j1] = y0r - y2r; a[j1 + 1] = y0i - y2i; y0r = wd3i * x1r + wd3r * x1i; y0i = wd3i * x1i - wd3r * x1r; y2r = wk3i * x3r + wk3r * x3i; y2i = wk3i * x3i - wk3r * x3r; a[j2] = y0r + y2r; a[j2 + 1] = y0i + y2i; a[j3] = y0r - y2r; a[j3 + 1] = y0i - y2i; } wk1r = w[m]; wk1i = w[m + 1]; j0 = mh; j1 = j0 + m; j2 = j1 + m; j3 = j2 + m; x0r = a[j0] - a[j2 + 1]; x0i = a[j0 + 1] + a[j2]; x1r = a[j0] + a[j2 + 1]; x1i = a[j0 + 1] - a[j2]; x2r = a[j1] - a[j3 + 1]; x2i = a[j1 + 1] + a[j3]; x3r = a[j1] + a[j3 + 1]; x3i = a[j1 + 1] - a[j3]; y0r = wk1r * x0r - wk1i * x0i; y0i = wk1r * x0i + wk1i * x0r; y2r = wk1i * x2r - wk1r * x2i; y2i = wk1i * x2i + wk1r * x2r; a[j0] = y0r + y2r; a[j0 + 1] = y0i + y2i; a[j1] = y0r - y2r; a[j1 + 1] = y0i - y2i; y0r = wk1i * x1r - wk1r * x1i; y0i = wk1i * x1i + wk1r * x1r; y2r = wk1r * x3r - wk1i * x3i; y2i = wk1r * x3i + wk1i * x3r; a[j2] = y0r - y2r; a[j2 + 1] = y0i - y2i; a[j3] = y0r + y2r; a[j3 + 1] = y0i + y2i;}
void cftfx41(int n, double *a, int nw, double *w) { void cftf161(double *a, double *w); void cftf162(double *a, double *w); void cftf081(double *a, double *w); void cftf082(double *a, double *w);
if (n == 128) { cftf161(a, &w[nw - 8]); cftf162(&a[32], &w[nw - 32]); cftf161(&a[64], &w[nw - 8]); cftf161(&a[96], &w[nw - 8]); } else { cftf081(a, &w[nw - 8]); cftf082(&a[16], &w[nw - 8]); cftf081(&a[32], &w[nw - 8]); cftf081(&a[48], &w[nw - 8]); }}
void cftf161(double *a, double *w) { double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
wn4r = w[1]; wk1r = w[2]; wk1i = w[3]; x0r = a[0] + a[16]; x0i = a[1] + a[17]; x1r = a[0] - a[16]; x1i = a[1] - a[17]; x2r = a[8] + a[24]; x2i = a[9] + a[25]; x3r = a[8] - a[24]; x3i = a[9] - a[25]; y0r = x0r + x2r; y0i = x0i + x2i; y4r = x0r - x2r; y4i = x0i - x2i; y8r = x1r - x3i; y8i = x1i + x3r; y12r = x1r + x3i; y12i = x1i - x3r; x0r = a[2] + a[18]; x0i = a[3] + a[19]; x1r = a[2] - a[18]; x1i = a[3] - a[19]; x2r = a[10] + a[26]; x2i = a[11] + a[27]; x3r = a[10] - a[26]; x3i = a[11] - a[27]; y1r = x0r + x2r; y1i = x0i + x2i; y5r = x0r - x2r; y5i = x0i - x2i; x0r = x1r - x3i; x0i = x1i + x3r; y9r = wk1r * x0r - wk1i * x0i; y9i = wk1r * x0i + wk1i * x0r; x0r = x1r + x3i; x0i = x1i - x3r; y13r = wk1i * x0r - wk1r * x0i; y13i = wk1i * x0i + wk1r * x0r; x0r = a[4] + a[20]; x0i = a[5] + a[21]; x1r = a[4] - a[20]; x1i = a[5] - a[21]; x2r = a[12] + a[28]; x2i = a[13] + a[29]; x3r = a[12] - a[28]; x3i = a[13] - a[29]; y2r = x0r + x2r; y2i = x0i + x2i; y6r = x0r - x2r; y6i = x0i - x2i; x0r = x1r - x3i; x0i = x1i + x3r; y10r = wn4r * (x0r - x0i); y10i = wn4r * (x0i + x0r); x0r = x1r + x3i; x0i = x1i - x3r; y14r = wn4r * (x0r + x0i); y14i = wn4r * (x0i - x0r); x0r = a[6] + a[22]; x0i = a[7] + a[23]; x1r = a[6] - a[22]; x1i = a[7] - a[23]; x2r = a[14] + a[30]; x2i = a[15] + a[31]; x3r = a[14] - a[30]; x3i = a[15] - a[31]; y3r = x0r + x2r; y3i = x0i + x2i; y7r = x0r - x2r; y7i = x0i - x2i; x0r = x1r - x3i; x0i = x1i + x3r; y11r = wk1i * x0r - wk1r * x0i; y11i = wk1i * x0i + wk1r * x0r; x0r = x1r + x3i; x0i = x1i - x3r; y15r = wk1r * x0r - wk1i * x0i; y15i = wk1r * x0i + wk1i * x0r; x0r = y12r - y14r; x0i = y12i - y14i; x1r = y12r + y14r; x1i = y12i + y14i; x2r = y13r - y15r; x2i = y13i - y15i; x3r = y13r + y15r; x3i = y13i + y15i; a[24] = x0r + x2r; a[25] = x0i + x2i; a[26] = x0r - x2r; a[27] = x0i - x2i; a[28] = x1r - x3i; a[29] = x1i + x3r; a[30] = x1r + x3i; a[31] = x1i - x3r; x0r = y8r + y10r; x0i = y8i + y10i; x1r = y8r - y10r; x1i = y8i - y10i; x2r = y9r + y11r; x2i = y9i + y11i; x3r = y9r - y11r; x3i = y9i - y11i; a[16] = x0r + x2r; a[17] = x0i + x2i; a[18] = x0r - x2r; a[19] = x0i - x2i; a[20] = x1r - x3i; a[21] = x1i + x3r; a[22] = x1r + x3i; a[23] = x1i - x3r; x0r = y5r - y7i; x0i = y5i + y7r; x2r = wn4r * (x0r - x0i); x2i = wn4r * (x0i + x0r); x0r = y5r + y7i; x0i = y5i - y7r; x3r = wn4r * (x0r - x0i); x3i = wn4r * (x0i + x0r); x0r = y4r - y6i; x0i = y4i + y6r; x1r = y4r + y6i; x1i = y4i - y6r; a[8] = x0r + x2r; a[9] = x0i + x2i; a[10] = x0r - x2r; a[11] = x0i - x2i; a[12] = x1r - x3i; a[13] = x1i + x3r; a[14] = x1r + x3i; a[15] = x1i - x3r; x0r = y0r + y2r; x0i = y0i + y2i; x1r = y0r - y2r; x1i = y0i - y2i; x2r = y1r + y3r; x2i = y1i + y3i; x3r = y1r - y3r; x3i = y1i - y3i; a[0] = x0r + x2r; a[1] = x0i + x2i; a[2] = x0r - x2r; a[3] = x0i - x2i; a[4] = x1r - x3i; a[5] = x1i + x3r; a[6] = x1r + x3i; a[7] = x1i - x3r;}
void cftf162(double *a, double *w) { double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, x0r, x0i, x1r, x1i, x2r, x2i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
wn4r = w[1]; wk1r = w[4]; wk1i = w[5]; wk3r = w[6]; wk3i = -w[7]; wk2r = w[8]; wk2i = w[9]; x1r = a[0] - a[17]; x1i = a[1] + a[16]; x0r = a[8] - a[25]; x0i = a[9] + a[24]; x2r = wn4r * (x0r - x0i); x2i = wn4r * (x0i + x0r); y0r = x1r + x2r; y0i = x1i + x2i; y4r = x1r - x2r; y4i = x1i - x2i; x1r = a[0] + a[17]; x1i = a[1] - a[16]; x0r = a[8] + a[25]; x0i = a[9] - a[24]; x2r = wn4r * (x0r - x0i); x2i = wn4r * (x0i + x0r); y8r = x1r - x2i; y8i = x1i + x2r; y12r = x1r + x2i; y12i = x1i - x2r; x0r = a[2] - a[19]; x0i = a[3] + a[18]; x1r = wk1r * x0r - wk1i * x0i; x1i = wk1r * x0i + wk1i * x0r; x0r = a[10] - a[27]; x0i = a[11] + a[26]; x2r = wk3i * x0r - wk3r * x0i; x2i = wk3i * x0i + wk3r * x0r; y1r = x1r + x2r; y1i = x1i + x2i; y5r = x1r - x2r; y5i = x1i - x2i; x0r = a[2] + a[19]; x0i = a[3] - a[18]; x1r = wk3r * x0r - wk3i * x0i; x1i = wk3r * x0i + wk3i * x0r; x0r = a[10] + a[27]; x0i = a[11] - a[26]; x2r = wk1r * x0r + wk1i * x0i; x2i = wk1r * x0i - wk1i * x0r; y9r = x1r - x2r; y9i = x1i - x2i; y13r = x1r + x2r; y13i = x1i + x2i; x0r = a[4] - a[21]; x0i = a[5] + a[20]; x1r = wk2r * x0r - wk2i * x0i; x1i = wk2r * x0i + wk2i * x0r; x0r = a[12] - a[29]; x0i = a[13] + a[28]; x2r = wk2i * x0r - wk2r * x0i; x2i = wk2i * x0i + wk2r * x0r; y2r = x1r + x2r; y2i = x1i + x2i; y6r = x1r - x2r; y6i = x1i - x2i; x0r = a[4] + a[21]; x0i = a[5] - a[20]; x1r = wk2i * x0r - wk2r * x0i; x1i = wk2i * x0i + wk2r * x0r; x0r = a[12] + a[29]; x0i = a[13] - a[28]; x2r = wk2r * x0r - wk2i * x0i; x2i = wk2r * x0i + wk2i * x0r; y10r = x1r - x2r; y10i = x1i - x2i; y14r = x1r + x2r; y14i = x1i + x2i; x0r = a[6] - a[23]; x0i = a[7] + a[22]; x1r = wk3r * x0r - wk3i * x0i; x1i = wk3r * x0i + wk3i * x0r; x0r = a[14] - a[31]; x0i = a[15] + a[30]; x2r = wk1i * x0r - wk1r * x0i; x2i = wk1i * x0i + wk1r * x0r; y3r = x1r + x2r; y3i = x1i + x2i; y7r = x1r - x2r; y7i = x1i - x2i; x0r = a[6] + a[23]; x0i = a[7] - a[22]; x1r = wk1i * x0r + wk1r * x0i; x1i = wk1i * x0i - wk1r * x0r; x0r = a[14] + a[31]; x0i = a[15] - a[30]; x2r = wk3i * x0r - wk3r * x0i; x2i = wk3i * x0i + wk3r * x0r; y11r = x1r + x2r; y11i = x1i + x2i; y15r = x1r - x2r; y15i = x1i - x2i; x1r = y0r + y2r; x1i = y0i + y2i; x2r = y1r + y3r; x2i = y1i + y3i; a[0] = x1r + x2r; a[1] = x1i + x2i; a[2] = x1r - x2r; a[3] = x1i - x2i; x1r = y0r - y2r; x1i = y0i - y2i; x2r = y1r - y3r; x2i = y1i - y3i; a[4] = x1r - x2i; a[5] = x1i + x2r; a[6] = x1r + x2i; a[7] = x1i - x2r; x1r = y4r - y6i; x1i = y4i + y6r; x0r = y5r - y7i; x0i = y5i + y7r; x2r = wn4r * (x0r - x0i); x2i = wn4r * (x0i + x0r); a[8] = x1r + x2r; a[9] = x1i + x2i; a[10] = x1r - x2r; a[11] = x1i - x2i; x1r = y4r + y6i; x1i = y4i - y6r; x0r = y5r + y7i; x0i = y5i - y7r; x2r = wn4r * (x0r - x0i); x2i = wn4r * (x0i + x0r); a[12] = x1r - x2i; a[13] = x1i + x2r; a[14] = x1r + x2i; a[15] = x1i - x2r; x1r = y8r + y10r; x1i = y8i + y10i; x2r = y9r - y11r; x2i = y9i - y11i; a[16] = x1r + x2r; a[17] = x1i + x2i; a[18] = x1r - x2r; a[19] = x1i - x2i; x1r = y8r - y10r; x1i = y8i - y10i; x2r = y9r + y11r; x2i = y9i + y11i; a[20] = x1r - x2i; a[21] = x1i + x2r; a[22] = x1r + x2i; a[23] = x1i - x2r; x1r = y12r - y14i; x1i = y12i + y14r; x0r = y13r + y15i; x0i = y13i - y15r; x2r = wn4r * (x0r - x0i); x2i = wn4r * (x0i + x0r); a[24] = x1r + x2r; a[25] = x1i + x2i; a[26] = x1r - x2r; a[27] = x1i - x2i; x1r = y12r + y14i; x1i = y12i - y14r; x0r = y13r - y15i; x0i = y13i + y15r; x2r = wn4r * (x0r - x0i); x2i = wn4r * (x0i + x0r); a[28] = x1r - x2i; a[29] = x1i + x2r; a[30] = x1r + x2i; a[31] = x1i - x2r;}
void cftf081(double *a, double *w) { double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
wn4r = w[1]; x0r = a[0] + a[8]; x0i = a[1] + a[9]; x1r = a[0] - a[8]; x1i = a[1] - a[9]; x2r = a[4] + a[12]; x2i = a[5] + a[13]; x3r = a[4] - a[12]; x3i = a[5] - a[13]; y0r = x0r + x2r; y0i = x0i + x2i; y2r = x0r - x2r; y2i = x0i - x2i; y1r = x1r - x3i; y1i = x1i + x3r; y3r = x1r + x3i; y3i = x1i - x3r; x0r = a[2] + a[10]; x0i = a[3] + a[11]; x1r = a[2] - a[10]; x1i = a[3] - a[11]; x2r = a[6] + a[14]; x2i = a[7] + a[15]; x3r = a[6] - a[14]; x3i = a[7] - a[15]; y4r = x0r + x2r; y4i = x0i + x2i; y6r = x0r - x2r; y6i = x0i - x2i; x0r = x1r - x3i; x0i = x1i + x3r; x2r = x1r + x3i; x2i = x1i - x3r; y5r = wn4r * (x0r - x0i); y5i = wn4r * (x0r + x0i); y7r = wn4r * (x2r - x2i); y7i = wn4r * (x2r + x2i); a[8] = y1r + y5r; a[9] = y1i + y5i; a[10] = y1r - y5r; a[11] = y1i - y5i; a[12] = y3r - y7i; a[13] = y3i + y7r; a[14] = y3r + y7i; a[15] = y3i - y7r; a[0] = y0r + y4r; a[1] = y0i + y4i; a[2] = y0r - y4r; a[3] = y0i - y4i; a[4] = y2r - y6i; a[5] = y2i + y6r; a[6] = y2r + y6i; a[7] = y2i - y6r;}
void cftf082(double *a, double *w) { double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
wn4r = w[1]; wk1r = w[2]; wk1i = w[3]; y0r = a[0] - a[9]; y0i = a[1] + a[8]; y1r = a[0] + a[9]; y1i = a[1] - a[8]; x0r = a[4] - a[13]; x0i = a[5] + a[12]; y2r = wn4r * (x0r - x0i); y2i = wn4r * (x0i + x0r); x0r = a[4] + a[13]; x0i = a[5] - a[12]; y3r = wn4r * (x0r - x0i); y3i = wn4r * (x0i + x0r); x0r = a[2] - a[11]; x0i = a[3] + a[10]; y4r = wk1r * x0r - wk1i * x0i; y4i = wk1r * x0i + wk1i * x0r; x0r = a[2] + a[11]; x0i = a[3] - a[10]; y5r = wk1i * x0r - wk1r * x0i; y5i = wk1i * x0i + wk1r * x0r; x0r = a[6] - a[15]; x0i = a[7] + a[14]; y6r = wk1i * x0r - wk1r * x0i; y6i = wk1i * x0i + wk1r * x0r; x0r = a[6] + a[15]; x0i = a[7] - a[14]; y7r = wk1r * x0r - wk1i * x0i; y7i = wk1r * x0i + wk1i * x0r; x0r = y0r + y2r; x0i = y0i + y2i; x1r = y4r + y6r; x1i = y4i + y6i; a[0] = x0r + x1r; a[1] = x0i + x1i; a[2] = x0r - x1r; a[3] = x0i - x1i; x0r = y0r - y2r; x0i = y0i - y2i; x1r = y4r - y6r; x1i = y4i - y6i; a[4] = x0r - x1i; a[5] = x0i + x1r; a[6] = x0r + x1i; a[7] = x0i - x1r; x0r = y1r - y3i; x0i = y1i + y3r; x1r = y5r - y7r; x1i = y5i - y7i; a[8] = x0r + x1r; a[9] = x0i + x1i; a[10] = x0r - x1r; a[11] = x0i - x1i; x0r = y1r + y3i; x0i = y1i - y3r; x1r = y5r + y7r; x1i = y5i + y7i; a[12] = x0r - x1i; a[13] = x0i + x1r; a[14] = x0r + x1i; a[15] = x0i - x1r;}
void cftf040(double *a) { double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
x0r = a[0] + a[4]; x0i = a[1] + a[5]; x1r = a[0] - a[4]; x1i = a[1] - a[5]; x2r = a[2] + a[6]; x2i = a[3] + a[7]; x3r = a[2] - a[6]; x3i = a[3] - a[7]; a[0] = x0r + x2r; a[1] = x0i + x2i; a[2] = x1r - x3i; a[3] = x1i + x3r; a[4] = x0r - x2r; a[5] = x0i - x2i; a[6] = x1r + x3i; a[7] = x1i - x3r;}
void cftb040(double *a) { double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
x0r = a[0] + a[4]; x0i = a[1] + a[5]; x1r = a[0] - a[4]; x1i = a[1] - a[5]; x2r = a[2] + a[6]; x2i = a[3] + a[7]; x3r = a[2] - a[6]; x3i = a[3] - a[7]; a[0] = x0r + x2r; a[1] = x0i + x2i; a[2] = x1r + x3i; a[3] = x1i - x3r; a[4] = x0r - x2r; a[5] = x0i - x2i; a[6] = x1r - x3i; a[7] = x1i + x3r;}
void cftx020(double *a) { double x0r, x0i;
x0r = a[0] - a[2]; x0i = a[1] - a[3]; a[0] += a[2]; a[1] += a[3]; a[2] = x0r; a[3] = x0i;}
void rftfsub(int n, double *a, int nc, double *c) { int j, k, kk, ks, m; double wkr, wki, xr, xi, yr, yi;
m = n >> 1; ks = 2 * nc / m; kk = 0; for (j = 2; j < m; j += 2) { k = n - j; kk += ks; wkr = 0.5 - c[nc - kk]; wki = c[kk]; xr = a[j] - a[k]; xi = a[j + 1] + a[k + 1]; yr = wkr * xr - wki * xi; yi = wkr * xi + wki * xr; a[j] -= yr; a[j + 1] -= yi; a[k] += yr; a[k + 1] -= yi; }}
void rftbsub(int n, double *a, int nc, double *c) { int j, k, kk, ks, m; double wkr, wki, xr, xi, yr, yi;
m = n >> 1; ks = 2 * nc / m; kk = 0; for (j = 2; j < m; j += 2) { k = n - j; kk += ks; wkr = 0.5 - c[nc - kk]; wki = c[kk]; xr = a[j] - a[k]; xi = a[j + 1] + a[k + 1]; yr = wkr * xr + wki * xi; yi = wkr * xi - wki * xr; a[j] -= yr; a[j + 1] -= yi; a[k] += yr; a[k + 1] -= yi; }}
void dctsub(int n, double *a, int nc, double *c) { int j, k, kk, ks, m; double wkr, wki, xr;
m = n >> 1; ks = nc / n; kk = 0; for (j = 1; j < m; j++) { k = n - j; kk += ks; wkr = c[kk] - c[nc - kk]; wki = c[kk] + c[nc - kk]; xr = wki * a[j] - wkr * a[k]; a[j] = wkr * a[j] + wki * a[k]; a[k] = xr; } a[m] *= c[0];}
void dstsub(int n, double *a, int nc, double *c) { int j, k, kk, ks, m; double wkr, wki, xr;
m = n >> 1; ks = nc / n; kk = 0; for (j = 1; j < m; j++) { k = n - j; kk += ks; wkr = c[kk] - c[nc - kk]; wki = c[kk] + c[nc - kk]; xr = wki * a[k] - wkr * a[j]; a[k] = wkr * a[k] + wki * a[j]; a[j] = xr; } a[m] *= c[0];}