Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[10744] trunk/octave-forge/extra/lssa
|
From: <be...@us...> - 2012-07-17 00:18:12
|
Revision: 10744
http://octave.svn.sourceforge.net/octave/?rev=10744&view=rev
Author: benjf5
Date: 2012-07-17 00:18:06 +0000 (Tue, 17 Jul 2012)
Log Message:
-----------
fastlscomplex has been mostly rewritten, and has problems with rounding errors (but generates some correct results!).
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/SampleScriptWithVostokData
trunk/octave-forge/extra/lssa/fastlscomplex.cc
Modified: trunk/octave-forge/extra/lssa/SampleScriptWithVostokData
===================================================================
--- trunk/octave-forge/extra/lssa/SampleScriptWithVostokData 2012-07-15 20:13:22 UTC (rev 10743)
+++ trunk/octave-forge/extra/lssa/SampleScriptWithVostokData 2012-07-17 00:18:06 UTC (rev 10744)
@@ -35,17 +35,190 @@
ls_complex_co2_gas_age = lscomplex(co2(:,3),co2(:,4),1,100,20);
ls_real_co2_ice_age = lsreal(co2(:,2),co2(:,4),1,100,20);
ls_real_co2_gas_age = lsreal(co2(:,3),co2(:,4),1,100,20);
+ls_complex_ch4_ice_age = lscomplex(ch4(:,2),ch4(:,4),1,100,20);
+ls_complex_ch4_gas_age = lscomplex(ch4(:,3),ch4(:,4),1,100,20);
+ls_real_ch4_ice_age = lsreal(ch4(:,2),ch4(:,4),1,100,20);
+ls_real_ch4_gas_age = lsreal(ch4(:,3),ch4(:,4),1,100,20);
+ls_complex_o18_ice_age = lscomplex(o18(:,2),o18(:,4),1,100,20);
+ls_complex_o18_gas_age = lscomplex(o18(:,3),o18(:,4),1,100,20);
+ls_real_o18_ice_age = lsreal(o18(:,2),o18(:,4),1,100,20);
+ls_real_o18_gas_age = lsreal(o18(:,3),o18(:,4),1,100,20);
+ls_complex_deut = lscomplex(deut(:,2),deut(:,3),1,100,20);
+ls_real_deut = lsreal(deut(:,2),deut(:,3),1,100,20);
+ls_complex_dust = lscomplex(dust(:,2),dust(:,3),1,100,20);
+ls_real_dust = lsreal(dust(:,2),dust(:,3),1,100,20);
+x_data_axis_vector = [ -430000, 0 ];
+## Useful because all of the data extends over 430 000 years up to the
+## present.
+
+## Setting up the CO2 plots:
figure(co2_fig);
subplot(4,2,1);
+axis(x_data_axis_vector);
plot(-(co2(:,2)),co2(:,4));
title("Gas levels over ice age");
subplot(4,2,2);
+axis(x_data_axis_vector);
plot(-(co2(:,3),co2(:,4));
title("Gas levels over gas age");
subplot(4,2,3);
+plot(real(ls_complex_co2_ice_age));
+hold on;
+plot(imag(ls_complex_co2_ice_age),'r');
+title("Complex L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,4);
+plot(real(ls_complex_co2_gas_age));
+hold on;
+plot(imag(ls_complex_co2_gas_age),'r');
+title("Complex L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+subplot(4,2,5);
+plot(real(ls_real_co2_ice_age));
+hold on;
+plot(imag(ls_real_co2_ice_age),'r');
+title("Real L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,6);
+plot(real(ls_real_co2_gas_age));
+hold on;
+plot(imag(ls_real_co2_gas_age));
+title("Real L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+## At this point, we have transforms of both datasets, real and complex,
+## and just need to figure out what cool thing to do with the remaining slot.
+## Setting up the CH4 plots
+figure(ch4_fig);
+subplot(4,2,1);
+axis(x_data_axis_vector);
+plot(-(ch4(:,2)),ch4(:,4));
+title("Gas levels over ice age");
+subplot(4,2,2);
+axis(x_data_axis_vector);
+plot(-(ch4(:,3),ch4(:,4));
+title("Gas levels over gas age");
+subplot(4,2,3);
+plot(real(ls_complex_ch4_ice_age));
+hold on;
+plot(imag(ls_complex_ch4_ice_age),'r');
+title("Complex L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,4);
+plot(real(ls_complex_ch4_gas_age));
+hold on;
+plot(imag(ls_complex_ch4_gas_age),'r');
+title("Complex L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+subplot(4,2,5);
+plot(real(ls_real_ch4_ice_age));
+hold on;
+plot(imag(ls_real_ch4_ice_age),'r');
+title("Real L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,6);
+plot(real(ls_real_ch4_gas_age));
+hold on;
+plot(imag(ls_real_ch4_gas_age));
+title("Real L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+## Setting up the O18 plots:
+figure(o18_fig);
+subplot(4,2,1);
+axis(x_data_axis_vector);
+plot(-(o18(:,2)),o18(:,4));
+title("Gas levels over ice age");
+subplot(4,2,2);
+axis(x_data_axis_vector);
+plot(-(o18(:,3),o18(:,4));
+title("Gas levels over gas age");
+subplot(4,2,3);
+plot(real(ls_complex_o18_ice_age));
+hold on;
+plot(imag(ls_complex_o18_ice_age),'r');
+title("Complex L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,4);
+plot(real(ls_complex_o18_gas_age));
+hold on;
+plot(imag(ls_complex_o18_gas_age),'r');
+title("Complex L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+subplot(4,2,5);
+plot(real(ls_real_o18_ice_age));
+hold on;
+plot(imag(ls_real_o18_ice_age),'r');
+title("Real L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,6);
+plot(real(ls_real_o18_gas_age));
+hold on;
+plot(imag(ls_real_o18_gas_age));
+title("Real L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+
+## Setting up Dust plots:
+figure(dust_fig);
+subplot(4,1,1);
+axis(x_data_axis_vector);
+plot(-(dust(:,2)),dust(:,3));
+title("Dust levels over ice age");
+subplot(4,1,2);
+plot(real(ls_complex_dust_ice_age));
+hold on;
+plot(imag(ls_complex_dust_ice_age),'r');
+title("Complex L-S transform of Dust/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,1,3);
+plot(real(ls_real_dust_ice_age));
+hold on;
+plot(imag(ls_real_dust_ice_age),'r');
+title("Real L-S transform of Dust/ice age data");
+legend("Real part","Imaginary part");
+
+## Setting up Deuterium plots:
+figure(deut_fig);
+subplot(4,1,1);
+axis(x_data_axis_vector);
+plot(-(deut(:,2)),deut(:,3));
+title("Deuterium levels over ice age");
+subplot(4,1,2);
+plot(real(ls_complex_deut_ice_age));
+hold on;
+plot(imag(ls_complex_deut_ice_age),'r');
+title("Complex L-S transform of Deuterium/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,1,3);
+plot(real(ls_real_deut_ice_age));
+hold on;
+plot(imag(ls_real_deut_ice_age),'r');
+title("Real L-S transform of Deuterium/ice age data");
+legend("Real part","Imaginary part");
+
+co2_ch4_comparison_figure = figure("visible","off","name","CO2/CH4
+comparison");
+subplot(4,1,1);
+axes(x_data_axis_vector);
+plot(-(co2(:,2)),co2(:,4));
+hold on;
+plot(-(ch4(:,2)),ch4(:,4),'g');
+title("CO2 and CH4 data");
+legend("CO2","CH4");
+
+subplot(4,1,2);
+plot(abs(ls_complex_co2_ice_age));
+hold on;
+plot(abs(ls_complex_ch4_gas_age),'g');
+title("Abs. values of CO2 and CH4 L-S complex transforms");
+legend("CO2,CH4");
+
+
+
+
+
+
## to implement:
## - displays of all the data and flaws in trying to model with just
## using L-S data
Modified: trunk/octave-forge/extra/lssa/fastlscomplex.cc
===================================================================
--- trunk/octave-forge/extra/lssa/fastlscomplex.cc 2012-07-15 20:13:22 UTC (rev 10743)
+++ trunk/octave-forge/extra/lssa/fastlscomplex.cc 2012-07-17 00:18:06 UTC (rev 10744)
@@ -11,6 +11,8 @@
#include <iostream>
#include <exception>
+
+
ComplexRowVector flscomplex( RowVector tvec , ComplexRowVector xvec ,
double maxfreq , int octaves , int coefficients);
@@ -67,14 +69,13 @@
ComplexRowVector results = ComplexRowVector (coefficients * octaves );
- double tau, delta_tau, tau_0, tau_h, n_inv, mu,
+ double tau, delta_tau, tau_0, tau_h, n_inv, mu, te,
omega_oct, omega_multiplier, octavemax, omega_working,
- loop_tau_0, loop_delta_tau;
+ loop_tau_0, loop_delta_tau, on_1, n_1, o;
double length = ( tvec((tvec.numel()-1)) - tvec( octave_idx_type (0)));
int octave_iter, coeff_iter;
- std::complex<double> zeta, z_accumulator, exp_term, exp_multiplier, alpha;
+ std::complex<double> zeta, zz, z_accumulator, exp_term, exp_multiplier, alpha, h, *tpra, *temp_ptr_alpha, temp_alpha[12], *tprb, *temp_ptr_beta, temp_beta[12], temp_array[12], *p, x;
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 );
}
@@ -83,7 +84,7 @@
for ( int i = 1 ; i < 12 ; i++ ) {
factorial_array[i] = factorial_array[i-1] * i;
}
- n_inv = 1.0 / n;
+ n_1 = n_inv = 1.0 / n;
mu = (0.5 * M_PI)/length; // Per the article; this is in place to improve numerical accuracy if desired.
/* Viz. the paper, in which Dtau = c / omega_max, and c is stated as pi/2 for floating point processors,
* In the case of this computation, I'll go by the recommendation.
@@ -91,64 +92,118 @@
delta_tau = M_PI / ( 2 * maxfreq );
tau_0 = tvec(0) + delta_tau;
tau_h = tau_0;
- size_t precomp_subset_count = (size_t) ceil( ( tvec(tvec.numel()-1) - tvec(0) ) / ( 2 * delta_tau ) );
- // I've used size_t because it will work for my purposes without threatening undefined behaviour.
- const std::complex<double> im = std::complex<double> ( 0 , 1 );
+ te = tau_h + delta_tau;
octave_idx_type k ( 0 ); // Iterator for accessing xvec, tvec.
Precomputation_Record * precomp_records_head, *record_current, *record_tail, *record_ref, *record_next;
record_current = precomp_records_head = new Precomputation_Record;
- for ( size_t p_r_iter = 1 ; p_r_iter < precomp_subset_count ; p_r_iter++ ) {
+ for ( te = tvec(k) + (2 * delta_tau) ; ; ) {
+ x = std::complex<double>(xvec(k));
+ {
+ double t = mu*(tvec(k)-tau_h), tt;
+ p = record_current->power_series;
+ // p = 0
+ *p++ = std::complex<double>(x);
+ // p = 1
+ tt = -t;
+ h = x * tt;
+ *p++ = std::complex<double>(-h.imag(),h.real());
+ // p = 2
+ tt *= t*(1.0/2.0);
+ *p++ = x*tt;
+ // p = 3
+ tt *= t*(-1.0/3.0);
+ h = x * tt;
+ *p++ = std::complex<double>(-h.imag(),h.real());
+ // p = 4
+ tt *= t*(1.0/4.0);
+ *p++ = x*tt;
+ // p = 5
+ tt *= t*(-1.0/5.0);
+ h = x * tt;
+ *p++ = std::complex<double>(-h.imag(),h.real());
+ // p = 6
+ tt *= t*(1.0/6.0);
+ *p++ = x*tt;
+ // p = 7
+ tt *= t*(-1.0/7.0);
+ h = x * tt;
+ *p++ = std::complex<double>(-h.imag(),h.real());
+ // p = 8
+ tt *= t*(1.0/8.0);
+ *p++ = x*tt;
+ // p = 9
+ tt *= t*(-1.0/9.0);
+ h = x * tt;
+ *p++ = std::complex<double>(-h.imag(),h.real());
+ // p = 10
+ tt *= t*(1.0/10.0);
+ *p++ = x*tt;
+ // p = 11
+ tt *= t*(-1.0/11.0);
+ h = x * tt;
+ *p++ = std::complex<double>(-h.imag(),h.real());
+ }
+ record_current->stored_data = true;
+ for(k++; ( k < n ) && tvec(k) < te ; k++ ) {
+ x = std::complex<double>(xvec(k));
+ {
+ double t = mu*(tvec(k)-tau_h), tt;
+ p = record_current->power_series;
+ // p = 0
+ *p++ += std::complex<double>(x);
+ // p = 1
+ tt = -t;
+ h = x * tt;
+ *p++ += std::complex<double>(- h.imag(), h.real());
+ // p = 2
+ tt *= t*(1.0/2.0);
+ *p++ += x*tt;
+ // p = 3
+ tt *= t*(-1.0/3.0);
+ h = x * tt;
+ *p++ += std::complex<double>(-h.imag(),h.real());
+ // p = 4
+ tt *= t*(1.0/4.0);
+ *p++ += x*tt;
+ // p = 5
+ tt *= t*(-1.0/5.0);
+ h = x * tt;
+ *p++ += std::complex<double>(-h.imag(),h.real());
+ // p = 6
+ tt *= t*(1.0/6.0);
+ *p++ += x*tt;
+ // p = 7
+ tt *= t*(-1.0/7.0);
+ h = x * tt;
+ *p++ += std::complex<double>(-h.imag(),h.real());
+ // p = 8
+ tt *= t*(1.0/8.0);
+ *p++ += x*tt;
+ // p = 9
+ tt *= t*(-1.0/9.0);
+ h = x * tt;
+ *p++ += std::complex<double>(-h.imag(),h.real());
+ // p = 10
+ tt *= t*(1.0/10.0);
+ *p++ += x*tt;
+ // p = 11
+ tt *= t*(-1.0/11.0);
+ h = x * tt;
+ *p++ += std::complex<double>(-h.imag(),h.real());
+ }
+ record_current->stored_data = true;
+ }
+ 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;
- /* 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( ; record_current != 0 ; record_current = record_current->next ) {
- for ( int j = 0 ; j < 12 ; j++ ) {
- record_current->power_series[j] = std::complex<double> ( 0 , 0 );
- } // To avoid any trouble down the line, although it is an annoyance.
- // Error is traced this far. Time to see if it's in this loop.
- while ( (k < n) && (abs(tvec(k)-tau_h) <= delta_tau) ) {
- double p;
- // alpha = std::complex<double> ( 0 , 0 );
- for ( int j = 0 ; j < 12 ; j++ ) {
- alpha.real() = xvec(k).real();
- alpha.imag() = xvec(k).imag();
- // alpha = std::complex<double> ( xvec(k).real() , xvec(k).imag() );
- // alpha *= ( -1 * im * mu * ( tvec(k) - tau_h ) ) * p;
- 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];
- }
- }
- record_current->power_series[j].real() += alpha.real();
- 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.
- record_current->stored_data = true;
- k++;
- }
- tau_h += ( 2 * delta_tau );
- }
- // At this point all precomputation records have been exhausted; short-circuit is abused to avoid overflow errors.
/* 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)
@@ -162,40 +217,37 @@
octave_idx_type iter ( 0 );
- double real_part = 0, imag_part = 0, real_part_accumulator = 0, 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;
- exp_term = std::complex<double> ( cos( - omega_working * loop_tau_0 ) ,
- sin ( - omega_working * loop_tau_0 ) );
- 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){
- real_part_accumulator = 0;
- imag_part_accumulator = 0;
- real_part = 0;
- imag_part = 0;
- for ( record_current = precomp_records_head ; record_current ;
+ 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_multiplier = std::complex<double> ( cos ( - 2 * omega_working * loop_delta_tau ) ,
+ sin ( - 2 * omega_working * loop_delta_tau ) );
+ for ( zeta = 0, record_current = precomp_records_head ; record_current ;
record_current = record_current->next, exp_term *= exp_multiplier ) {
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- z_accumulator = ( pow(omega_working,array_iter) * record_current->power_series[array_iter] );
- real_part_accumulator += z_accumulator.real();
- imag_part_accumulator += z_accumulator.imag();
+ if ( record_current->stored_data ) {
+ int p;
+ for ( zz = 0 , p = 0, on_1 = n_1 ; p < 12 ; p++ ) {
+ zz += record_current->power_series[p] * on_1 ;
+ on_1 *= o;
+ }
+ zeta += exp_term * zz;
}
- real_part = real_part + ( exp_term.real() * real_part_accumulator - ( exp_term.imag() * imag_part_accumulator ) );
- imag_part = imag_part + ( exp_term.imag() * real_part_accumulator + exp_term.real() * imag_part_accumulator );
}
- results(iter) = std::complex<double> ( n_inv * real_part , n_inv * imag_part );
+ results(iter) = std::complex<double> (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 ) );
@@ -204,61 +256,92 @@
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 ) {
- 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) * record_current->power_series[p - ( 2 * q )];
+ 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++) = std::complex<double>(record_current->power_series[p]);
+ 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) * record_current->power_series[p - ( 2 * q ) - 1];
- }
}
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- record_current->power_series[array_iter].real() = temp[array_iter].real();
- record_current->power_series[array_iter].imag() = temp[array_iter].imag();
- }
if ( ! record_ref ) break; // Last record was reached
}
else {
record_next = record_ref;
if ( 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) * ( record_current->power_series[p - ( 2 * q )] - record_next->power_series[p - (2*q)] );
+ 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) * ( record_current->power_series[p - ( 2 * q ) - 1] - record_next->power_series[p - ( 2 * q ) - 1 ] );
- }
}
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- record_current->power_series[array_iter].real() = temp[array_iter].real();
- 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) * 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)) * record_next->power_series[p - ( 2 * q ) - 1];
- }
}
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- record_current->power_series[array_iter].real() = temp[array_iter].real();
- record_current->power_series[array_iter].imag() = temp[array_iter].imag();
- }
}
record_current->stored_data = true;
record_ref = record_next;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
