Mercurial > forge
changeset 10677:08dd8ccec62f octave-forge
Removed looped lsreal, lscomplex; detabified fastlsreal.cc, cleaned up spaces in fastlscomplex.
| author | benjf5 |
|---|---|
| date | Tue, 14 Aug 2012 18:55:38 +0000 |
| parents | 82432e2a8d82 |
| children | 4cf5fa6c82d8 |
| files | extra/lssa/inst/lscomplex.m.old extra/lssa/inst/lsreal.m.old extra/lssa/src/fastlscomplex.cc extra/lssa/src/fastlsreal.cc |
| diffstat | 4 files changed, 257 insertions(+), 406 deletions(-) [+] |
line wrap: on
line diff
--- a/extra/lssa/inst/lscomplex.m.old Tue Aug 14 16:51:37 2012 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,68 +0,0 @@ -## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com> -## -## This program is free software; you can redistribute it and/or modify it under -## the terms of the GNU General Public License as published by the Free Software -## Foundation; either version 3 of the License, or (at your option) any later -## version. -## -## This program is distributed in the hope that it will be useful, but WITHOUT -## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -## details. -## -## You should have received a copy of the GNU General Public License along with -## this program; if not, see <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{t} =} lscomplex (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{numoctaves}) -## -## Return the complex least-squares transform of the (@var{time},@var{mag}) -## series, considering frequencies up to @var{maxfreq}, over @var{numoctaves} -## octaves and @var{numcoeff} coefficients. -## -## @seealso{lsreal} -## @end deftypefn - - -function transform = lscomplex (t, x, omegamax, ncoeff, noctave) - - ## VECTOR ONLY, and since t and x have the same number of - ## entries, there's no problem. - n = length (t); - - - transform = zeros (1, ncoeff * noctave); - - o = omegamax; - - omul = 2 ^ (- 1 / ncoeff); - - for iter = 1:ncoeff * noctave - - ot = o .* t; - - ## See the paper for the expression below - transform(iter) = sum ((cos (ot) - (sin (ot) .* i)) .* x) / n; - - - ## Advance the transform to the next coefficient in the octave - o *= omul; - - endfor - -endfunction - -%!test -%! maxfreq = 4 / ( 2 * pi ); -%! t = [0:0.008:8]; -%! x = ( 2 .* sin (maxfreq .* t) + -%! 3 .* sin ( (3 / 4) * maxfreq .* t)- -%! 0.5 .* sin ((1/4) * maxfreq .* t) - -%! 0.2 .* cos (maxfreq .* t) + -%! cos ((1/4) * maxfreq .* t)); -%! o = [ maxfreq , 3 / 4 * maxfreq , 1 / 4 * maxfreq ]; -%! assert (lscomplex (t, x, maxfreq, 2, 2), -%! [(-0.400924546169395 - 2.371555305867469i), ... -%! (1.218065147708429 - 2.256125004156890i), ... -%! (1.935428592212907 - 1.539488163739336i), ... -%! (2.136692292751917 - 0.980532175174563i)], 5e-10);
--- a/extra/lssa/inst/lsreal.m.old Tue Aug 14 16:51:37 2012 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,79 +0,0 @@ -## Copyright (C) 2012 Benjamin Lewis -## -## This program is free software; you can redistribute it and/or modify it under -## the terms of the GNU General Public License as published by the Free Software -## Foundation; either version 3 of the License, or (at your option) any later -## version. -## -## This program is distributed in the hope that it will be useful, but WITHOUT -## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -## details. -## -## You should have received a copy of the GNU General Public License along with -## this program; if not, see <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{transform} =} lsreal (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{numoctaves}) -## -## Return the real least-squares transform of the time series -## defined, based on the maximal frequency @var{maxfreq}, the -## number of coefficients @var{numcoeff}, and the number of -## octaves @var{numoctaves}. Each complex-valued result is the -## pair (c_o, s_o) defining the coefficients which best fit the -## function y = c_o * cos(ot) + s_o * sin(ot) to the (@var{time}, @var{mag}) data. -## -## @seealso{lscomplex} -## @end deftypefn - -function transform = lsreal (t, x, omegamax, ncoeff, noctave) - - ## FIXME : THIS IS VECTOR-ONLY. I'd need to add another bit of code to - ## make it array-safe, and that's not knowing right now what else - ## will be necessary. - k = n = length (t); - - transform = zeros (1, (noctave * ncoeff)); - - od = 2 ^ (- 1 / ncoeff); - o = omegamax; - n1 = 1 / n; - - ncoeffp = ncoeff * noctave; - - for iter = 1:ncoeffp - ## This method is an application of Eq. 8 on - ## page 6 of the text, as well as Eq. 7 - ot = o .* t; - - zeta = n1 * sum ((cos (ot) - i * sin (ot)) .* x); - - ot *= 2; - - iota = n1 * sum (cos (ot) - i * sin (ot)); - - - transform(iter) = (2 * (conj (zeta) - (conj (iota) * zeta)) / - (1 - (real (iota) ^ 2) - (imag (iota) ^ 2))); - - o *= od; - endfor - -endfunction - -%!test -%! maxfreq = 4 / ( 2 * pi ); -%! t = linspace(0,8); -%! x = ( 2 .* sin ( maxfreq .* t ) + -%! 3 .* sin ( (3/4) * maxfreq .* t ) - -%! 0.5 .* sin ( (1/4) * maxfreq .* t ) - -%! 0.2 .* cos ( maxfreq .* t ) + -%! cos ( (1/4) * maxfreq .* t ) ); -%! # In the assert here, I've got an error bound large enough to catch -%! # individual system errors which would present no real issue. -%! assert (lsreal (t,x,maxfreq,2,2), -%! [(-1.68275915310663 + 4.70126183846743i), ... -%! (1.93821553170889 + 4.95660209883437i), ... -%! (4.38145452686697 + 2.14403733658600i), ... -%! (5.27425332281147 - 0.73933440226597i)], -%! 5e-10)
--- a/extra/lssa/src/fastlscomplex.cc Tue Aug 14 16:51:37 2012 +0000 +++ b/extra/lssa/src/fastlscomplex.cc Tue Aug 14 18:55:38 2012 +0000 @@ -43,12 +43,12 @@ \n\ @seealso{lscomplex, fastlsreal}\n\ \n\ -@end deftypefn") +@end deftypefn") { - octave_value_list retval; + octave_value_list retval; - if (args.length() != 5) + if (args.length() != 5) print_usage(); else { @@ -60,38 +60,37 @@ int ncoeff = args(4).int_value (); if (tvals.numel () != xvals.numel ()) - if (tvals.numel () > xvals.numel ()) + if (tvals.numel () > xvals.numel ()) error ("fastlscomplex: More time values than magnitude values"); - else + else error ("fastlscomplex: More magnitude values than time values"); - - if (ncoeff == 0) + if (ncoeff == 0) error ("fastlscomplex: No coefficients to compute"); - if (noctaves == 0) + if (noctaves == 0) error ("fastlscomplex: No octaves to compute over"); - if (omegamax == 0) + if (omegamax == 0) error ("fastlscomplex: No difference between minimal and maximal frequency"); - if (! error_state) + if (! error_state) { ComplexRowVector results; if (flscomplex (tvals, xvals, omegamax, noctaves, ncoeff, results)) retval(0) = octave_value (results); else error ("fastlscomplex: error in the underlying flscomplex function"); - } - + } + } return retval; } bool flscomplex (const RowVector & tvec, const ComplexRowVector & xvec, - double maxfreq, int octaves, int coefficients, - ComplexRowVector & results) + double maxfreq, int octaves, int coefficients, + ComplexRowVector & results) { - - struct Precomputation_Record + + struct Precomputation_Record { Precomputation_Record *next; std::complex<double> power_series[12]; // I'm using 12 as a matter of compatibility, only. @@ -112,19 +111,19 @@ int octave_iter, coeff_iter; - std::complex<double> zeta, zz, z_accumulator, exp_term, exp_multiplier, - alpha, h, *tpra, *temp_ptr_alpha, temp_alpha[12], *tprb, *temp_ptr_beta, + 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 (); - for (int array_iter = 0; array_iter < 12; array_iter++) + for (int array_iter = 0; array_iter < 12; array_iter++) temp_array[array_iter] = std::complex<double> (0 , 0); - + int factorial_array[12]; factorial_array[0] = 1; - for (int i = 1; i < 12; i++) + for (int i = 1; i < 12; i++) factorial_array[i] = factorial_array[i-1] * i; n_1 = n_inv = 1.0 / n; @@ -244,7 +243,7 @@ record_current->stored_data = true; } - if (k >= n) + if (k >= n) break; tau_h = te + delta_tau; @@ -257,7 +256,7 @@ 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) */ @@ -267,18 +266,18 @@ octavemax = maxfreq; loop_tau_0 = tau_0; loop_delta_tau = delta_tau; - + octave_idx_type iter = 0; - + // Loops need to first travel over octaves, then coefficients; - for (octave_iter = octaves; ; + for (octave_iter = octaves; ; omega_oct *= 0.5, octavemax *= 0.5, loop_tau_0 += loop_delta_tau, loop_delta_tau *= 2) { o = omega_oct; omega_working = octavemax; - for (coeff_iter = 0; - coeff_iter < coefficients; - coeff_iter++, o *= omega_multiplier, + 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), @@ -287,14 +286,14 @@ 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; + for (zeta = 0, record_current = precomp_records_head; record_current; record_current = record_current->next, exp_term *= exp_multiplier ) { - if (record_current->stored_data) + if (record_current->stored_data) { int p; - for (zz = 0, p = 0, on_1 = n_1; p < 12; p++) + for (zz = 0, p = 0, on_1 = n_1; p < 12; p++) { zz += record_current->power_series[p] * on_1; on_1 *= o; @@ -302,40 +301,39 @@ zeta += exp_term * zz; } } - results(iter) = std::complex<double> (zeta); iter++; } - if (! (--octave_iter)) + if (! (--octave_iter)) break; - + /* If we've already reached the lowest value, stop. * Otherwise, merge with the next computation range. */ 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++) + + 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 (record_current = precomp_records_head; + { + for (record_current = precomp_records_head; record_current; - record_current = record_current->next) + record_current = record_current->next) { - if (! (record_ref = record_current->next ) - || ! record_ref->stored_data ) + 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) + 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++) + p < 12; + p++ , exp_pse_ptr++) { tpra = temp_ptr_alpha; *(temp_ptr_alpha++) = std::complex<double>(record_current->power_series[p]); @@ -348,29 +346,29 @@ * actually be accessed for the first iterations, it's easier. */ - if (++exp_ptr >= exp_pse_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 ) + + if (++exp_ptr >= exp_pse_ptr ) break; --tpra; record_current->power_series[p] -= *tpra * *exp_ptr; - if (++exp_ptr >= exp_pse_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) + + if (++exp_ptr >= exp_pse_ptr) break; --tpra; @@ -378,112 +376,112 @@ } } - if ( ! record_ref ) + if ( ! record_ref ) break; // Last record was reached - } - else + } + else { record_next = record_ref; - if ( record_current->stored_data ) + 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++) + 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 ) + 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 ) + + if ( ++exp_ptr >= exp_pse_ptr ) break; - + tpra -= 2; record_current->power_series[p] -= *tpra * *exp_ptr; - - if (++exp_ptr >= exp_pse_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) + + if (++exp_ptr >= exp_pse_ptr) break; - + tpra -= 2; record_current->power_series[p] += *tpra * *exp_ptr; } } - } + } else { 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++) + 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, + + for (exp_ptr = exp_power_series_elements, record_next->power_series[q] = *tpra * *exp_ptr; ;) { - if (++exp_ptr >= exp_pse_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) + + if (++exp_ptr >= exp_pse_ptr) break; - + --tpra; record_next->power_series[q] -= *tpra * *exp_ptr; - - if ( ++exp_ptr >= exp_pse_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) + + if (++exp_ptr >= exp_pse_ptr) break; - + --tpra; record_next->power_series[q] += *tpra * *exp_ptr; } } } - + record_current->stored_data = true; record_ref = record_next; record_current->next = record_ref->next; delete record_ref; - + } } } - } - catch (std::exception & e) + } + catch (std::exception & e) {//This section was part of my debugging, and may be removed. std::cout << "Exception thrown: " << e.what() << std::endl; return (false);
--- a/extra/lssa/src/fastlsreal.cc Tue Aug 14 16:51:37 2012 +0000 +++ b/extra/lssa/src/fastlsreal.cc Tue Aug 14 18:55:38 2012 +0000 @@ -25,11 +25,11 @@ #include <exception> ComplexRowVector flsreal( RowVector tvec , ComplexRowVector xvec , - double maxfreq , int octaves , int coefficients); + double maxfreq , int octaves , int coefficients); DEFUN_DLD(fastlsreal,args,nargout, - "-*- texinfo -*-\n\ + "-*- texinfo -*-\n\ @deftypefn {Function File} { C = } fastlsreal(@var{time},@var{magnitude},@var{maximum_frequency},@var{octaves},@var{coefficients})\n\ \n\ Return the real least-sqaures spectral fit to the (@var{time},@var{magnitude})\n\ @@ -68,7 +68,7 @@ } ComplexRowVector flsreal( RowVector tvec , RowVector xvec , - double maxfreq, int octaves, int coefficients ) { + double maxfreq, int octaves, int coefficients ) { struct XTElem { double x, t; }; @@ -169,55 +169,55 @@ 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; + 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; } record_current->stored_data = true; } @@ -250,33 +250,33 @@ 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 ) ); + 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 ) ); + 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, + 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; - } + 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; + } } results(iter) = 2 / ( 1 - ( iota.real() * iota.real() ) - (iota.imag() * - iota.imag() ) - * ( conj(zeta) - conj(iota) * zeta ); + iota.imag() ) + * ( conj(zeta) - conj(iota) * zeta ); iter++; } if ( !(--octave_iter) ) break; @@ -288,104 +288,104 @@ 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 ) ); + * ( mu * loop_delta_tau) * ( 1.0 / ( (double) r_iter ) ); } try{ 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 , temp_ptr_alpha = temp_alpha ; ; ) { - 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_pse_ptr++, record_current->power_series[p]->x = tpra->x * *exp_ptr, record_current->power_series[p]->t = tpra->t * *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_power_series_elements ) break; - ++tpra; - record_current->power_series[p]->x -= tpra->x * *exp_ptr; - record_current->power_series[p]->t -= tpra->t * *exp_ptr; - if ( --exp_ptr < exp_power_series_elements ) break; - ++tpra; - record_current->power_series[p]->x += tpra->x * *exp_ptr; - record_current->power_series[p]->t += tpra->x * *exp_ptr; - } - if ( ++p >= 12 ) break; - 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( tprb = temp_beta, exp_ptr = exp_pse_ptr++, record_current->power_series[p]->t = tprb->t * *exp_ptr; exp_ptr > exp_power_series_elements ; ) { - ++tprb; - --exp_ptr; - record_current->power_series[p]->t += tprb->t * *exp_ptr; - } - if ( ++p >= 12 ) break; - } - } - if ( ! record_ref ) break; // Last record was reached - } - else { - record_next = record_ref; - if ( record_current->stored_data ) { - int p = 0; - for( exp_pse_ptr = exp_power_series_elements, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; ; ) { - temp_ptr_alpha->x = record_current->power_series[p]->x + record_next->power_series[p]->x; - (temp_ptr_alpha++)->t = record_current->power_series[p]->t + record_next->power_series[p]->t; - temp_ptr_beta->x = record_ref->power_series[p]->x - record_current->power_series[p]->x; - (temp_ptr_beta++)->t = record_ref->power_series[p]->t - record_current->power_series[p]->t; - for( tpra = temp_alpha, exp_ptr = exp_pse_ptr++, record_current->power_series[p]->x = tpra->x * *exp_ptr, record_current->power_series[p]->t = tpra->x * *exp_ptr; ; ) { - if ( --exp_ptr < exp_pse_ptr ) break; - ++tpra; - record_current->power_series[p]->x -= tpra->x * *exp_ptr; - record_current->power_series[p]->t -= tpra->t * *exp_ptr; - if ( --exp_ptr < exp_pse_ptr ) break; - ++tpra; - record_current->power_series[p]->x += tpra->x * *exp_ptr; - record_current->power_series[p]->t += tpra->t * *exp_ptr; - } - if ( ++p >= 12 ) break; - temp_ptr_alpha->x = record_next->power_series[p]->x - record_current->power_series[p]->x; - (temp_ptr_alpha++)->t = record_next->power_series[p]->t - record_current->power_series[p]->t; - temp_ptr_beta->x = record_current->power_series[p]->x + record_next->power_series[p]->x; - (temp_ptr_beta++)->t = record_current->power_series[p]->t + record_next->power_series[p]->t; - for(tprb = temp_beta, exp_ptr = exp_pse_ptr++, record_current->power_series[p]->x = tprb->x * *exp_ptr, record_current->power_series[p]->t = tprb->x * *exp_ptr; exp_ptr > exp_power_series_elements; ) { - ++tprb; - --exp_ptr; - record_current->power_series[p]->x += tprb->x * *exp_ptr; - record_current->power_series[p]->t += tprb->t * *exp_ptr; - } - if ( ++p >= 12 ) break; - } - } else { - int q = 0; - for( exp_pse_ptr = exp_power_series_elements, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; ; ) { - temp_ptr_alpha->x = record_next->power_series[q]->x; - temp_ptr_alpha->t = record_next->power_series[q]->t; - for(tpra = temp_alpha, exp_ptr = exp_pse_ptr++, record_next->power_series[q]->x = tpra->x * *exp_ptr, record_next->power_series[q]->t = tpra->t * *exp_ptr; exp_ptr > exp_power_series_elements; ) { - ++tpra; - --exp_ptr; - record_next->power_series[q]->x += tpra->x * *exp_ptr; - record_next->power_series[q]->t += tpra->t * *exp_ptr; - } - if ( ++q >= 12 ) break; - } - record_current->stored_data = true; - record_ref = record_next; - record_current->next = record_ref->next; - record_next = 0; - delete record_ref; - } + 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 , temp_ptr_alpha = temp_alpha ; ; ) { + 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_pse_ptr++, record_current->power_series[p]->x = tpra->x * *exp_ptr, record_current->power_series[p]->t = tpra->t * *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_power_series_elements ) break; + ++tpra; + record_current->power_series[p]->x -= tpra->x * *exp_ptr; + record_current->power_series[p]->t -= tpra->t * *exp_ptr; + if ( --exp_ptr < exp_power_series_elements ) break; + ++tpra; + record_current->power_series[p]->x += tpra->x * *exp_ptr; + record_current->power_series[p]->t += tpra->x * *exp_ptr; + } + if ( ++p >= 12 ) break; + 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( tprb = temp_beta, exp_ptr = exp_pse_ptr++, record_current->power_series[p]->t = tprb->t * *exp_ptr; exp_ptr > exp_power_series_elements ; ) { + ++tprb; + --exp_ptr; + record_current->power_series[p]->t += tprb->t * *exp_ptr; + } + if ( ++p >= 12 ) break; + } + } + if ( ! record_ref ) break; // Last record was reached + } + else { + record_next = record_ref; + if ( record_current->stored_data ) { + int p = 0; + for( exp_pse_ptr = exp_power_series_elements, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; ; ) { + temp_ptr_alpha->x = record_current->power_series[p]->x + record_next->power_series[p]->x; + (temp_ptr_alpha++)->t = record_current->power_series[p]->t + record_next->power_series[p]->t; + temp_ptr_beta->x = record_ref->power_series[p]->x - record_current->power_series[p]->x; + (temp_ptr_beta++)->t = record_ref->power_series[p]->t - record_current->power_series[p]->t; + for( tpra = temp_alpha, exp_ptr = exp_pse_ptr++, record_current->power_series[p]->x = tpra->x * *exp_ptr, record_current->power_series[p]->t = tpra->x * *exp_ptr; ; ) { + if ( --exp_ptr < exp_pse_ptr ) break; + ++tpra; + record_current->power_series[p]->x -= tpra->x * *exp_ptr; + record_current->power_series[p]->t -= tpra->t * *exp_ptr; + if ( --exp_ptr < exp_pse_ptr ) break; + ++tpra; + record_current->power_series[p]->x += tpra->x * *exp_ptr; + record_current->power_series[p]->t += tpra->t * *exp_ptr; + } + if ( ++p >= 12 ) break; + temp_ptr_alpha->x = record_next->power_series[p]->x - record_current->power_series[p]->x; + (temp_ptr_alpha++)->t = record_next->power_series[p]->t - record_current->power_series[p]->t; + temp_ptr_beta->x = record_current->power_series[p]->x + record_next->power_series[p]->x; + (temp_ptr_beta++)->t = record_current->power_series[p]->t + record_next->power_series[p]->t; + for(tprb = temp_beta, exp_ptr = exp_pse_ptr++, record_current->power_series[p]->x = tprb->x * *exp_ptr, record_current->power_series[p]->t = tprb->x * *exp_ptr; exp_ptr > exp_power_series_elements; ) { + ++tprb; + --exp_ptr; + record_current->power_series[p]->x += tprb->x * *exp_ptr; + record_current->power_series[p]->t += tprb->t * *exp_ptr; + } + if ( ++p >= 12 ) break; + } + } else { + int q = 0; + for( exp_pse_ptr = exp_power_series_elements, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; ; ) { + temp_ptr_alpha->x = record_next->power_series[q]->x; + temp_ptr_alpha->t = record_next->power_series[q]->t; + for(tpra = temp_alpha, exp_ptr = exp_pse_ptr++, record_next->power_series[q]->x = tpra->x * *exp_ptr, record_next->power_series[q]->t = tpra->t * *exp_ptr; exp_ptr > exp_power_series_elements; ) { + ++tpra; + --exp_ptr; + record_next->power_series[q]->x += tpra->x * *exp_ptr; + record_next->power_series[q]->t += tpra->t * *exp_ptr; + } + if ( ++q >= 12 ) break; + } + record_current->stored_data = true; + record_ref = record_next; + record_current->next = record_ref->next; + record_next = 0; + delete record_ref; + } } } return results;