00001
00002
00003
00004
00005
00006
00007 #ifdef STORM_USEF
00008 #include "rtrsmrartrswt_.h"
00009 #include "drtrsmrartrswt_.h"
00010 #include "irtrsmrartrswt_.h"
00011 #include "artrsmrartrswt_.h"
00012 #include "dartrsmrartrswt_.h"
00013 #include "iartrsmrartrswt_.h"
00014 #endif
00015
00016 #include "stormdef.h"
00017 #define INVSQRT2 STORMDEF_MATHCONST_INVSQRT2
00018
00019 void rtrsmrartrswt(
00020 const int& Msize,
00021 const int& mmax,
00022 int& mrastep,
00023 float* mracoef,
00024 float* mracoefrs,
00025 const float& datavalue
00026 )
00027 {
00028 #ifdef STORM_USEF
00029 rtrsmrartrswt_(
00030 &Msize,
00031 &mmax,
00032 &mrastep, mracoef, mracoefrs,
00033 &datavalue
00034 );
00035 #else
00036 mrastep++;
00037 int mcycle = 0;
00038 while ( ( mrastep & (1<<mcycle) ) && ( mcycle < mmax ) ) { mcycle++; }
00039 int mcyclers = mcycle + 1;
00040 if ( ( mrastep < (1<<mcyclers) ) || ( mcyclers > mmax ) ) { mcyclers = 0; }
00041 float tmp = datavalue;
00042 int m=1;
00043 float tmpold,tmps,tmpw;
00044 while ( m <= mcycle ) {
00045 tmpold = mracoef[ 2*(m-1) ];
00046 tmps = (tmpold+tmp) * INVSQRT2;
00047 tmpw = (tmpold-tmp) * INVSQRT2;
00048 mracoef[ 2*(m-1) ] = tmp;
00049 mracoef[ 2*m + 1 ] = tmpw;
00050 tmp=tmps;
00051 m++;
00052 }
00053 if ( m == mcyclers ) {
00054 tmpold = mracoef[ 2*(m-1) ];
00055 tmps = (tmpold+tmp) * INVSQRT2;
00056 tmpw = (tmpold-tmp) * INVSQRT2;
00057 mracoefrs[ 2*m ] = tmps;
00058 mracoefrs[ 2*m + 1 ] = tmpw;
00059 }
00060 mracoef[ 2*(m-1) ] = tmp;
00061 #endif
00062
00063 }
00064
00065
00066 void rtrsmrartrswt(
00067 const int& Msize,
00068 const int& mmax,
00069 int& mrastep,
00070 double* mracoef,
00071 double* mracoefrs,
00072 const double& datavalue
00073 )
00074 {
00075 #ifdef STORM_USEF
00076 drtrsmrartrswt_(
00077 &Msize,
00078 &mmax,
00079 &mrastep, mracoef, mracoefrs,
00080 &datavalue
00081 );
00082 #else
00083 mrastep++;
00084 int mcycle = 0;
00085 while ( ( mrastep & (1<<mcycle) ) && ( mcycle < mmax ) ) { mcycle++; }
00086 int mcyclers = mcycle + 1;
00087 if ( ( mrastep < (1<<mcyclers) ) || ( mcyclers > mmax ) ) { mcyclers = 0; }
00088 double tmp = datavalue;
00089 int m=1;
00090 double tmpold,tmps,tmpw;
00091 while ( m <= mcycle ) {
00092 tmpold = mracoef[ 2*(m-1) ];
00093 tmps = (tmpold+tmp) * INVSQRT2;
00094 tmpw = (tmpold-tmp) * INVSQRT2;
00095 mracoef[ 2*(m-1) ] = tmp;
00096 mracoef[ 2*m + 1 ] = tmpw;
00097 tmp = tmps;
00098 m++;
00099 }
00100 if ( m == mcyclers ) {
00101 tmpold = mracoef[ 2*(m-1) ];
00102 tmps = (tmpold+tmp) * INVSQRT2;
00103 tmpw = (tmpold-tmp) * INVSQRT2;
00104 mracoefrs[ 2*m ] = tmps;
00105 mracoefrs[ 2*m + 1 ] = tmpw;
00106 }
00107 mracoef[ 2*(m-1) ]=tmp;
00108 #endif
00109 }
00110
00111
00112 void rtrsmrartrswt(
00113 const int& Msize,
00114 const int& mmax,
00115 int& mrastep,
00116 int* mracoef,
00117 int* mracoefrs,
00118 const int& datavalue
00119 )
00120 {
00121 #ifdef STORM_USEF
00122 irtrsmrartrswt_( &Msize, &mmax, &mrastep, mracoef, mracoefrs, &datavalue );
00123 #else
00124 mrastep++;
00125 int mcycle = 0;
00126 while ( ( mrastep & (1<<mcycle) ) && ( mcycle < mmax ) ) { mcycle++; }
00127 int mcyclers = mcycle + 1;
00128 if ( ( mrastep < (1<<mcyclers) ) || ( mcyclers > mmax ) ) { mcyclers = 0; }
00129 int tmp = datavalue;
00130 int m=1;
00131 int tmpold,tmprs;
00132 while ( m <= mcycle ) {
00133 tmpold = mracoef[ 2*(m-1) ];
00134 mracoef[ 2*(m-1) ] = tmp;
00135 mracoef[ 2*m + 1 ] = tmpold-tmp;
00136 tmp = tmpold+tmp;
00137 if ( (tmp >= 0) || ((tmp & 0x00000001) != 0 ) ) {
00138 tmp = tmp/2;
00139 } else {
00140 tmp = tmp/2-1;
00141 }
00142 m++;
00143 }
00144 if ( m == mcyclers ) {
00145 tmpold = mracoef[ 2*(m-1) ];
00146 mracoefrs[ 2*m + 1 ] = tmpold-tmp;
00147 tmprs = tmpold+tmp;
00148 if ( (tmprs >= 0) || ((tmprs & 0x00000001) != 0 ) ) {
00149 tmprs = tmprs/2;
00150 } else {
00151 tmprs = tmprs/2-1;
00152 }
00153 mracoefrs[ 2*m ] = tmprs;
00154 }
00155 mracoef[ 2*(m-1) ]=tmp;
00156 #endif
00157 }
00158
00159
00160 void rtrsmrartrswt(
00161 const int& Asize,
00162 const int& Msize,
00163 const int& amax,
00164 const int& mmax,
00165 int& mrastep,
00166 float* mracoef,
00167 float* mracoefrs,
00168 const float* datavalue
00169 )
00170 {
00171 #ifdef STORM_USEF
00172 artrsmrartrswt_(
00173 &Asize, &Msize,
00174 &amax, &mmax,
00175 &mrastep, mracoef, mracoefrs,
00176 datavalue
00177 );
00178 #else
00179 mrastep++;
00180 int mcycle = 0;
00181 while ( ( mrastep & (1<<mcycle) ) && ( mcycle < mmax ) ) { mcycle++; }
00182 int mcyclers = mcycle + 1;
00183 if ( ( mrastep < (1<<mcyclers) ) || ( mcyclers > mmax ) ) { mcyclers = 0; }
00184 for( int a=0; a <= amax; a++ ) {
00185 float tmp = datavalue[a];
00186 int m = 1;
00187 float tmpold,tmps,tmpw;
00188 while ( m <= mcycle ) {
00189 tmpold = mracoef[ Msize*2*a + 2*(m-1) ];
00190 tmps = (tmpold+tmp) * INVSQRT2;
00191 tmpw = (tmpold-tmp) * INVSQRT2;
00192 mracoef[ Msize*2*a + 2*(m-1) ] = tmp;
00193 mracoef[ Msize*2*a + 2*m + 1 ] = tmpw;
00194 tmp = tmps;
00195 m++;
00196 }
00197 if ( m == mcyclers ) {
00198 tmpold = mracoef[ Msize*2*a + 2*(m-1) ];
00199 tmps = (tmpold+tmp) * INVSQRT2;
00200 tmpw = (tmpold-tmp) * INVSQRT2;
00201 mracoefrs[ Msize*2*a + 2*m ] = tmps;
00202 mracoefrs[ Msize*2*a + 2*m + 1 ] = tmpw;
00203 }
00204 mracoef[ Msize*2*a + 2*(m-1) ] = tmp;
00205 }
00206 #endif
00207 }
00208
00209
00210 void rtrsmrartrswt(
00211 const int& Asize,
00212 const int& Msize,
00213 const int& amax,
00214 const int& mmax,
00215 int& mrastep,
00216 double* mracoef,
00217 double* mracoefrs,
00218 const double* datavalue
00219 )
00220 {
00221 #ifdef STORM_USEF
00222 dartrsmrartrswt_(
00223 &Asize, &Msize,
00224 &amax, &mmax,
00225 &mrastep, mracoef, mracoefrs,
00226 datavalue
00227 );
00228 #else
00229 mrastep++;
00230 int mcycle = 0;
00231 while ( ( mrastep & (1<<mcycle) ) && ( mcycle < mmax ) ) { mcycle++; }
00232 int mcyclers = mcycle+1;
00233 if ( ( mrastep < (1<<mcyclers) ) || ( mcyclers > mmax ) ) { mcyclers = 0; }
00234 for( int a = 0; a <= amax; a++) {
00235 double tmp = datavalue[a];
00236 int m = 1;
00237 double tmpold,tmps,tmpw;
00238 while ( m <= mcycle ) {
00239 tmpold = mracoef[ Msize*2*a + 2*(m-1) ];
00240 tmps = (tmpold+tmp) * INVSQRT2;
00241 tmpw = (tmpold-tmp) * INVSQRT2;
00242 mracoef[ Msize*2*a + 2*(m-1) ] = tmp;
00243 mracoef[ Msize*2*a + 2*m + 1 ] = tmpw;
00244 tmp = tmps;
00245 m++;
00246 }
00247 if ( m == mcyclers ) {
00248 tmpold = mracoef[ Msize*2*a + 2*(m-1) ];
00249 tmps = (tmpold+tmp) * INVSQRT2;
00250 tmpw = (tmpold-tmp) * INVSQRT2;
00251 mracoefrs[ Msize*2*a + 2*m ] = tmps;
00252 mracoefrs[ Msize*2*a + 2*m + 1 ] = tmpw;
00253 }
00254 mracoef[ Msize*2*a + 2*(m-1) ] = tmp;
00255 }
00256 #endif
00257 }
00258
00259
00260 void rtrsmrartrswt(
00261 const int& Asize,
00262 const int& Msize,
00263 const int& amax,
00264 const int& mmax,
00265 int& mrastep,
00266 int* mracoef,
00267 int* mracoefrs,
00268 const int* datavalue
00269 )
00270 {
00271 #ifdef STORM_USEF
00272 iartrsmrartrswt_(
00273 &Asize, &Msize,
00274 &amax, &mmax,
00275 &mrastep, mracoef, mracoefrs,
00276 datavalue
00277 );
00278 #else
00279 mrastep++;
00280 int mcycle = 0;
00281 while ( ( mrastep & (1<<mcycle) ) && ( mcycle < mmax ) ) { mcycle++; }
00282 int mcyclers = mcycle+1;
00283 if ( ( mrastep < (1<<mcyclers) ) || ( mcyclers > mmax ) ) { mcyclers = 0; }
00284 for( int a = 0; a <= amax; a++) {
00285 int tmp = datavalue[a];
00286 int m = 1;
00287 int tmpold,tmprs;
00288 while ( m <= mcycle ) {
00289 tmpold = mracoef[ Msize*2*a + 2*(m-1) ];
00290 mracoef[ Msize*2*a + 2*(m-1) ] = tmp;
00291 mracoef[ Msize*2*a + 2*m + 1 ] = tmpold-tmp;
00292 tmp = tmpold+tmp;
00293 if ( (tmp >= 0) || ((tmp & 0x00000001) != 0 ) ) {
00294 tmp = tmp/2;
00295 } else {
00296 tmp = tmp/2-1;
00297 }
00298 m++;
00299 }
00300 if ( m == mcyclers ) {
00301 tmpold = mracoef[ Msize*2*a + 2*(m-1) ];
00302 mracoefrs[ Msize*2*a + 2*m + 1 ] = tmpold-tmp;
00303 tmprs = tmpold+tmp;
00304 if ( (tmprs >= 0) || ((tmprs & 0x00000001) != 0 ) ) {
00305 tmprs = tmprs/2;
00306 } else {
00307 tmprs = tmprs/2-1;
00308 }
00309 mracoefrs[ Msize*2*a + 2*m ] = tmprs;
00310 }
00311 mracoef[ Msize*2*a + 2*(m-1) ] = tmp;
00312 }
00313 #endif
00314 }
00315
