Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[10785] trunk/octave-forge/extra/lssa/src
|
From: <be...@us...> - 2012-07-30 18:00:55
|
Revision: 10785
http://octave.svn.sourceforge.net/octave/?rev=10785&view=rev
Author: benjf5
Date: 2012-07-30 18:00:45 +0000 (Mon, 30 Jul 2012)
Log Message:
-----------
Some modifications made, still getting everything ready.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/src/fastlsreal.cc
Added Paths:
-----------
trunk/octave-forge/extra/lssa/src/Makefile
Added: trunk/octave-forge/extra/lssa/src/Makefile
===================================================================
--- trunk/octave-forge/extra/lssa/src/Makefile (rev 0)
+++ trunk/octave-forge/extra/lssa/src/Makefile 2012-07-30 18:00:45 UTC (rev 10785)
@@ -0,0 +1,14 @@
+MKOCTFILE ?= mkoctfile
+
+all: fastlscomplex.oct \
+ fastlsreal.oct
+
+fastlscomplex.oct: fastlscomplex.cc
+ $(MKOCTFILE) fastlscomplex.cc
+
+fastlsreal.oct: fastlsreal.cc
+ $(MKOCTFILE) fastlsreal.cc
+
+# helper function just in case
+clean:
+ rm *.o *.oct *~ octave-core
Modified: trunk/octave-forge/extra/lssa/src/fastlsreal.cc
===================================================================
--- trunk/octave-forge/extra/lssa/src/fastlsreal.cc 2012-07-29 21:51:22 UTC (rev 10784)
+++ trunk/octave-forge/extra/lssa/src/fastlsreal.cc 2012-07-30 18:00:45 UTC (rev 10785)
@@ -55,11 +55,14 @@
}
-ComplexRowVector flsreal( RowVector tvec , ComplexRowVector xvec ,
+ComplexRowVector flsreal( RowVector tvec , RowVector xvec ,
double maxfreq, int octaves, int coefficients ) {
+ struct XTElem {
+ double x, t;
+ };
struct Precomputation_Record {
Precomputation_Record *next;
- std::complex<double> power_series[12]; // I'm using 12 as a matter of compatibility, only.
+ XTElem power_series[12]; // I'm using 12 as a matter of compatibility, only.
bool stored_data;
};
@@ -67,16 +70,15 @@
double tau, delta_tau, tau_0, tau_h, n_inv, mu,
omega_oct, omega_multiplier, octavemax, omega_working,
- loop_tau_0, loop_delta_tau;
+ loop_tau_0, loop_delta_tau, x;
double length = ( tvec((tvec.numel()-1)) - tvec( octave_idx_type (0)));
int octave_iter, coeff_iter;
std::complex<double> zeta, z_accumulator, zeta_exp_term, zeta_exp_multiplier, alpha,
- iota, i_accumulator, iota_exp_term, iota_exp_multiplier;
+ iota, i_accumulator, iota_exp_term, iota_exp_multiplier, exp_squared, exp_squared_multiplier;
octave_idx_type n = tvec.numel();
- std::complex<double> temp_array[12];
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- temp_array[array_iter] = std::complex<double> ( 0 , 0 );
- }
+ XTElem *tpra, *temp_ptr_alpha, temp_alpha[12], *tprb, *temp_ptr_beta, temp_beta[12], temp_array[12];
+
+
int factorial_array[12];
factorial_array[0] = 1;
for ( int i = 1 ; i < 12 ; i++ ) {
@@ -95,104 +97,127 @@
const std::complex<double> im = std::complex<double> ( 0 , 1 ); //I seriously prefer C99's complex numbers.
octave_idx_type k ( 0 ); // Iterator for accessing xvec, tvec.
-
- Precomputation_Record * complex_precomp_records_head, *complex_record_current,
- *complex_record_tail, *complex_record_ref, *complex_record_next, *iota_precomp_records_head,
- *iota_record_current, *iota_record_tail, *iota_record_ref, *iota_record_next;
- complex_record_current = complex_precomp_records_head = new Precomputation_Record;
- iota_record_current = iota_precomp_records_head = new Precomputation_Record;
- for ( size_t p_r_iter = 1 ; p_r_iter < precomp_subset_count ; p_r_iter++ ) {
- complex_record_current->next = new Precomputation_Record;
- iota_record_current->next = new Precomputation_Record;
- complex_record_current = complex_record_current->next;
- iota_record_current = iota_record_current->next;
- }
- // Error's past this point
- complex_record_tail = complex_record_current;
- iota_record_tail = iota_record_current;
- complex_record_current = complex_precomp_records_head;
- iota_record_current = iota_precomp_records_head;
- complex_record_tail->next = 0;
- iota_record_tail->next = 0;
- /* A test needs to be included for if there was a failure, but since
- * precomp_subset_count is of type size_t, it should be okay. */
- for( ; complex_record_current != 0 ; complex_record_current = complex_record_current->next ) {
- for ( int j = 0 ; j < 12 ; j++ ) {
- complex_record_current->power_series[j] = std::complex<double> ( 0 , 0 );
- } // To avoid any trouble down the line, although it is an annoyance.
- while ( (k < n) && (abs(tvec(k)-tau_h) <= delta_tau) ) {
- double p;
- for ( int j = 0 ; j < 12 ; j++ ) {
- alpha.real() = xvec(k).real();
- alpha.imag() = xvec(k).imag();
- // Change to switches for easier manipulation, fewer tests. This is two tests where one will do.
- if ( !( j % 2 ) ) {
- if ( ! ( j % 4 ) ) {
- alpha.real() = xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- } else {
- alpha.real() = -1 * xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = -1 * xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- }
- } else {
- if ( ! ( j % 3 ) ) {
- alpha.real() = -1 * xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = -1 * xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- } else {
- alpha.real() = xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- }
- }
- complex_record_current->power_series[j].real() += alpha.real();
- complex_record_current->power_series[j].imag() += alpha.imag();
- }
- // Computes each next step of the power series for the given power series element.
- // j was reused since it was a handy inner-loop variable, even though I used it twice here.
- complex_record_current->stored_data = true;
- k++;
+
+ Precomputation_Record * precomp_records_head, *record_current, *record_tail, *record_ref, *record_next;
+ record_current = precomp_records_head = new Precomputation_Record;
+ for ( te = tvec(k) + (2 * delta_tau) ; ; ) {
+ x = xvec(k);
+ {
+ double t = mu*(tvec(k)-tau_h), tt;
+ p = record_current->power_series;
+ // p = 0
+ p->x = x;
+ (p++)->t = 1;
+ // p = 1
+ tt = -t;
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 2
+ tt *= t*(1.0/2.0);
+ p->x = x*tt;
+ (p++)->t = tt;
+ // p = 3
+ tt *= t*(-1.0/3.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 4
+ tt *= t*(1.0/4.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 5
+ tt *= t*(-1.0/5.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 6
+ tt *= t*(1.0/6.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 7
+ tt *= t*(-1.0/7.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 8
+ tt *= t*(1.0/8.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 9
+ tt *= t*(-1.0/9.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 10
+ tt *= t*(1.0/10.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 11
+ tt *= t*(-1.0/11.0);
+ p->x = x * tt;
+ (p++)->t = tt;
}
- tau_h += ( 2 * delta_tau );
- }
- // At this point all precomputation records have been
- // exhausted for complex records. Short-circuit is abused
- // to avoid overflow errors.
- // Reset k, tau_h to reset the process. I may rewrite
- // these loops to be one, since running twice as long to
- // do the same thing is painful. May also move to a switch
- // in the previous section too.
- k = 0;
- tau_h = tau_0;
- for( ; iota_record_current != 0 ; iota_record_current = iota_record_current->next ) {
- for ( int j = 0 ; j < 12 ; j++ ) {
- iota_record_current->power_series[j] = std::complex<double> ( 0 , 0 );
- }
- while( ( k < n ) && (abs(tvec(k)-tau_h) <= delta_tau) ) {
- double comps[12];
- iota_record_current->power_series[0].real() = 1;
- comps[0] = 1;
- for ( int j = 1 ; j < 12 ; j++ ) {
- comps[j] = comps[j-1] * mu * (tvec(k)-tau_h);
- switch ( j % 4 ) {
- case 0 :
- iota_record_current->power_series[j].real() += comps[j] / factorial_array[j] ;
- break;
- case 1:
- iota_record_current->power_series[j].imag() += comps[j] / factorial_array[j] ;
- break;
- case 2:
- iota_record_current->power_series[j].real() -= comps[j] / factorial_array[j] ;
- break;
- case 3:
- iota_record_current->power_series[j].imag() -= comps[j] / factorial_array[j] ;
- break;
- }
+ record_current->stored_data = true;
+ for(k++; ( k < n ) && tvec(k) < te ; k++ ) {
+ x = xvec(k);
+ {
+ double t = mu*(tvec(k)-tau_h), tt;
+ p = record_current->power_series;
+ // p = 0
+ p->x += x;
+ (p++)->t += 1;
+ // p = 1
+ tt = -t;
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 2
+ tt *= t*(1.0/2.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 3
+ tt *= t*(-1.0/3.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 4
+ tt *= t*(1.0/4.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 5
+ tt *= t*(-1.0/5.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 6
+ tt *= t*(1.0/6.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 7
+ tt *= t*(-1.0/7.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 8
+ tt *= t*(1.0/8.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 9
+ tt *= t*(-1.0/9.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 10
+ tt *= t*(1.0/10.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 11
+ tt *= t*(-1.0/11.0);
+ p->x += x * tt;
+ (p++)->t += tt;
}
- iota_record_current->stored_data = true;
- k++;
+ record_current->stored_data = true;
}
- tau_h += ( 2 * delta_tau );
+ if( k >= n ) break;
+ tau_h = te + delta_tau;
+ te = tau_h + delta_tau;
+ record_current->next = new Precomputation_Record;
+ record_current = record_current->next;
}
-
+ record_tail = record_current;
+ record_current = precomp_records_head;
+ record_tail->next = 0;
/* Summation of coefficients for each frequency. As we have ncoeffs * noctaves elements,
* it makes sense to work from the top down, as we have omega_max by default (maxfreq)
@@ -206,142 +231,155 @@
octave_idx_type iter ( 0 );
- double zeta_real_part = 0, zeta_imag_part = 0, zeta_real_part_accumulator = 0, zeta_imag_part_accumulator = 0,
- iota_real_part = 0, iota_imag_part = 0, iota_real_part_accumulator = 0, iota_imag_part_accumulator = 0;
-
// Loops need to first travel over octaves, then coefficients;
for ( octave_iter = octaves ; ; omega_oct *= 0.5 , octavemax *= 0.5 , loop_tau_0 += loop_delta_tau , loop_delta_tau *= 2 ) {
- omega_working = omega_oct;
- zeta_exp_term = std::complex<double> ( cos ( - omega_working * loop_tau_0 ) ,
- sin ( - omega_working * loop_tau_0 ) );
- zeta_exp_multiplier = std::complex<double> ( cos ( - 2 * omega_working * loop_delta_tau ) ,
- sin ( - 2 * omega_working * loop_delta_tau ) );
- iota_exp_term = std::complex<double> ( cos ( - 2 * omega_working * loop_tau_0 ) ,
- sin ( - 2 * omega_working * loop_tau_0 ) );
- iota_exp_multiplier = std::complex<double> ( cos ( - 2 * omega_working * loop_delta_tau ) ,
- sin ( - 2 * omega_working * loop_delta_tau ) );
- for ( coeff_iter = 0 ; coeff_iter < coefficients ; coeff_iter++, omega_working *= omega_multiplier){
- zeta_real_part_accumulator = 0;
- zeta_imag_part_accumulator = 0;
- zeta_real_part = 0;
- zeta_imag_part = 0;
- for ( complex_record_current = complex_precomp_records_head ; complex_record_current ;
- complex_record_current = complex_record_current->next, zeta_exp_term *= zeta_exp_multiplier ) {
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- z_accumulator = ( pow(omega_working,array_iter) * complex_record_current->power_series[array_iter] );
- zeta_real_part_accumulator += z_accumulator.real();
- zeta_imag_part_accumulator += z_accumulator.imag();
+ o = omega_oct;
+ omega_working = octavemax;
+ for ( coeff_iter = 0 ; coeff_iter < coefficients ; coeff_iter++, o *= omega_multiplier, omega_working *= omega_multiplier){
+ exp_term = std::complex<double> ( cos( - omega_working * loop_tau_0 ) ,
+ sin ( - omega_working * loop_tau_0 ) );
+ exp_squared = exp_term * exp_term;
+ exp_multiplier = std::complex<double> ( cos ( - 2 * omega_working * loop_delta_tau ) ,
+ sin ( - 2 * omega_working * loop_delta_tau ) );
+ exp_squared_multiplier = exp_multiplier * exp_multiplier;
+ for ( zeta = iota = 0, record_current = precomp_records_head ; record_current ;
+ record_current = record_current->next, exp_term *= exp_multiplier,
+ exp_squared *= exp_squared_multiplier ) {
+ if ( record_current->stored_data ) {
+ int p;
+ for ( zz = ii = 0 , p = 0, on_1 = n_1 ; p < 12 ; ) {
+ zz.real() += record_current->power_series[p]->x * on_1;
+ ii.real() += record_current->power_series[p++]-> t * o2n_1;
+ on_1 *= o;
+ o2n_1 *= o2;
+ zz.imag() += record_current->power_series[p]->x * on_1;
+ ii.imag() += record_current->power_series[p++]-> t * o2n_1;
+ on_1 *= o;
+ o2n_1 *= o2;
+ }
+ zeta += exp_term * zz;
+ iota += exp_squared * ii;
}
- zeta_real_part = zeta_real_part + ( zeta_exp_term.real() * zeta_real_part_accumulator - zeta_exp_term.imag() * zeta_imag_part_accumulator );
- zeta_imag_part = zeta_imag_part + ( zeta_exp_term.imag() * zeta_real_part_accumulator + zeta_exp_term.real() * zeta_imag_part_accumulator );
}
- for ( iota_record_current = iota_precomp_records_head; iota_record_current ;
- iota_record_current = iota_record_current->next, iota_exp_term *= iota_exp_multiplier ) {
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- i_accumulator = ( pow(omega_working,array_iter) * iota_record_current->power_series[array_iter] );
- iota_real_part_accumulator += i_accumulator.real();
- iota_imag_part_accumulator += i_accumulator.imag();
- }
- iota_real_part = iota_real_part + ( iota_exp_term.real() * iota_real_part_accumulator - iota_exp_term.imag() * iota_imag_part_accumulator );
- iota_imag_part = iota_imag_part + ( iota_exp_term.imag() * iota_real_part_accumulator + iota_exp_term.real() * iota_imag_part_accumulator );
- }
- // c + i s = 2 * ( zeta-omega-conjugate - iota-2omega-conjuage * zeta-omega ) / ( 1 - abs(iota-2omega) ^ 2 )
- // (is what the results will be.)
- zeta_real_part *= n_inv;
- zeta_imag_part *= n_inv;
- iota_real_part *= n_inv;
- iota_imag_part *= n_inv;
- double result_real_part = 2 * ( zeta_real_part - iota_real_part * zeta_real_part - iota_imag_part * zeta_imag_part )
- / ( 1 - iota_real_part * iota_real_part - iota_imag_part * iota_imag_part );
- double result_imag_part = 2 * ( - zeta_imag_part + iota_imag_part * zeta_real_part - iota_real_part * zeta_imag_part )
- / ( 1 - iota_real_part * iota_real_part - iota_imag_part * iota_imag_part );
- results(iter) = std::complex<double> ( result_real_part , result_imag_part );
+ results(iter) = 2 / ( 1 - ( iota.real() * iota.real() ) - (iota.imag() *
+ iota.imag() )
+ * ( conj(zeta) - conj(iota) * zeta );
iter++;
}
if ( !(--octave_iter) ) break;
/* If we've already reached the lowest value, stop.
* Otherwise, merge with the next computation range.
*/
- double exp_power_series_elements[12];
+ double *exp_pse_ptr, *exp_ptr, exp_power_series_elements[12];
exp_power_series_elements[0] = 1;
+ exp_pse_ptr = exp_ptr = exp_power_series_elements;
for ( int r_iter = 1 ; r_iter < 12 ; r_iter++ ) {
exp_power_series_elements[r_iter] = exp_power_series_elements[r_iter-1]
* ( mu * loop_delta_tau) * ( 1.0 / ( (double) r_iter ) );
}
try{
- for ( complex_record_current = complex_precomp_records_head ; complex_record_current ;
- complex_record_current = complex_record_current->next ) {
- if ( ! ( complex_record_ref = complex_record_current->next ) || ! complex_record_ref->stored_data ) {
- if ( complex_record_current->stored_data ) {
- std::complex<double> temp[12];
- for( int array_init = 0 ; array_init < 12 ; array_init++ ) { temp[array_init] = std::complex<double>(0,0); }
- for( int p = 0 ; p < 12 ; p ++ ) {
- double step_floor_r = floor( ( (double) p ) / 2.0 );
- double step_floor_i = floor( ( (double) ( p - 1 ) ) / 2.0 );
- for( int q = 0 ; q < step_floor_r ; q++ ) {
- temp[p] += exp_power_series_elements[2*q] * pow((double)-1,q) * complex_record_current->power_series[p - ( 2 * q )];
+ for ( record_current = precomp_records_head ; record_current ;
+ record_current = record_current->next ) {
+ if ( ! ( record_ref = record_current->next ) || ! record_ref->stored_data ) {
+ // In this case, there is no next record, but this record has data.
+ if ( record_current->stored_data ) {
+ int p = 0;
+ for( exp_pse_ptr = exp_power_series_elements + 1 , temp_ptr_alpha = temp_alpha ; p < 12 ; p++ , exp_pse_ptr++ ) {
+ tpra = temp_ptr_alpha;
+ temp_ptr_alpha->x = record_current->power_series[p]->x;
+ (temp_ptr_alpha++)->t = record_current->power_series[p]->t;
+ temp_ptr_beta->x = -record_current->power_series[p]->x;
+ (temp_ptr_beta++)->t = -record_current->power_series[p]->t;
+ for( exp_ptr = exp_power_series_elements, record_current->power_series[p] = *temp_ptr_alpha * *exp_ptr; ; ) {
+ /* This next block is from Mathias' code, and it does a few
+ * ... unsavoury things. First off, it uses conditionals with
+ * break in order to avoid potentially accessing null regions
+ * of memory, and then it does ... painful things with a few
+ * numbers. However, remembering that most of these will not
+ * actually be accessed for the first iterations, it's easier.
+ */
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ h = *tpra * *exp_ptr;
+ record_current->power_series[p].real() -= h.imag();
+ record_current->power_series[p].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ record_current->power_series[p] -= *tpra * *exp_ptr;
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ h = -*tpra * *exp_ptr;
+ record_current->power_series[p].real() -= h.imag();
+ record_current->power_series[p].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ record_current->power_series[p] += *tpra * *exp_ptr;
}
- for( int q = 0 ; q <= step_floor_i ; q++ ) {
- temp[p] += im * exp_power_series_elements[2 * q + 1] * pow((double)-1,q) * complex_record_current->power_series[p - ( 2 * q ) - 1];
- }
}
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- complex_record_current->power_series[array_iter].real() = temp[array_iter].real();
- complex_record_current->power_series[array_iter].imag() = temp[array_iter].imag();
- }
- if ( ! complex_record_ref ) break; // Last record was reached
+ if ( ! record_ref ) break; // Last record was reached
}
else {
- complex_record_next = complex_record_ref;
- if ( complex_record_current->stored_data ) {
- std::complex<double> temp[12];
- for( int array_init = 0 ; array_init < 12 ; array_init++ ) { temp[array_init] = std::complex<double>(0,0); }
- for( int p = 0 ; p < 12 ; p ++ ) {
- double step_floor_r = floor( ( (double) p ) / 2.0 );
- double step_floor_i = floor( ( (double) ( p - 1 ) ) / 2.0 );
- for( int q = 0 ; q < step_floor_r ; q++ ) {
- temp[p] += exp_power_series_elements[2*q] * pow((double)-1,q) * ( complex_record_current->power_series[p - ( 2 * q )] - complex_record_next->power_series[p - (2*q)] );
+ record_next = record_ref;
+ if ( record_current->stored_data ) {
+ int p = 0, q = 0;
+ for( exp_pse_ptr = exp_power_series_elements + 1, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; p < 12 ; p++, q++, exp_pse_ptr++ ) {
+ tpra = temp_ptr_alpha;
+ *temp_ptr_alpha++ = record_current->power_series[p] + record_next->power_series[q];
+ *temp_ptr_beta++ = record_current->power_series[p] - record_next->power_series[1];
+ tprb = temp_ptr_beta;
+ for( exp_ptr = exp_power_series_elements, record_current->power_series[p] = *tpra * *exp_ptr; ; ) {
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ tprb -= 2;
+ h = *tprb * *exp_ptr;
+ record_current->power_series[p].real() -= h.imag();
+ record_current->power_series[p].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ tpra -= 2;
+ record_current->power_series[p] -= *tpra * *exp_ptr;
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ tprb -= 2;
+ h = - *tprb * *exp_ptr;
+ record_current->power_series[p].real() -= h.imag();
+ record_current->power_series[p].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ tpra -= 2;
+ record_current->power_series[p] += *tpra * *exp_ptr;
}
- for( int q = 0 ; q <= step_floor_i ; q++ ) {
- temp[p] += im * exp_power_series_elements[2 * q + 1] * pow((double)-1,q) * ( complex_record_current->power_series[p - ( 2 * q ) - 1] - complex_record_next->power_series[p - ( 2 * q ) - 1 ] );
- }
}
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- complex_record_current->power_series[array_iter].real() = temp[array_iter].real();
- complex_record_current->power_series[array_iter].imag() = temp[array_iter].imag();
- }
} else {
- std::complex<double> temp[12];
- for( int array_init = 0 ; array_init < 12 ; array_init++ ) { temp[array_init] = std::complex<double>(0,0); }
- for( int p = 0 ; p < 12 ; p ++ ) {
- double step_floor_r = floor( ( (double) p ) / 2.0 );
- double step_floor_i = floor( ( (double) ( p - 1 ) ) / 2.0 );
- for( int q = 0 ; q < step_floor_r ; q++ ) {
- temp[p] += exp_power_series_elements[2*q] * pow((double)-1,q) * complex_record_next->power_series[p - ( 2 * q )];
+ int q = 0;
+ for( exp_pse_ptr = exp_power_series_elements + 1, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; q < 12; q++, exp_pse_ptr++ ) {
+ tpra = temp_ptr_alpha;
+ *temp_ptr_alpha++ = std::complex<double>(record_next->power_series[q]);
+ for ( exp_ptr = exp_power_series_elements, record_next->power_series[q] = *tpra * *exp_ptr; ; ) {
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ h = *tpra * *exp_ptr;
+ record_next->power_series[q].real() -= h.imag();
+ record_next->power_series[q].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ record_next->power_series[q] -= *tpra * *exp_ptr;
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ h = -*tpra * *exp_ptr;
+ record_next->power_series[q].real() -= h.imag();
+ record_next->power_series[q].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ record_next->power_series[q] += *tpra * *exp_ptr;
}
- for( int q = 0 ; q <= step_floor_i ; q++ ) {
- temp[p] += im * exp_power_series_elements[2 * q + 1] * pow((double)-1,(q+1)) * complex_record_next->power_series[p - ( 2 * q ) - 1];
- }
}
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- complex_record_current->power_series[array_iter].real() = temp[array_iter].real();
- complex_record_current->power_series[array_iter].imag() = temp[array_iter].imag();
- }
}
- complex_record_current->stored_data = true;
- complex_record_ref = complex_record_next;
- complex_record_current->next = complex_record_ref->next;
- delete complex_record_ref;
+ record_current->stored_data = true;
+ record_ref = record_next;
+ record_current->next = record_ref->next;
+ delete record_ref;
}
}
}
- } catch (std::exception& e) { //This section was part of my debugging, and may be removed.
- std::cout << "Exception thrown: " << e.what() << std::endl;
- ComplexRowVector exception_result (1);
- exception_result(0) = std::complex<double> ( 0,0);
- return exception_result;
- }
- }
+
+
return results;
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
