Re: octave-speed

[Paul Kienzle]

On Feb 5, 2005, at 1:54 AM, Paul Thomas wrote:

> (ii) Your sample programme can be rewritten as:
>
> tic;
> Ns=1e5;
> m=[1:Ns-1];
> T=sqrt(m.*(m+1));
> n=zeros(1,Ns);
> n(1)=8;
> for i=m
>  n(i+1)=sqrt(T(i)*n(i));
> end
> printf("time=%g\n",toc);
> printf("n=%g\n",n(1,Ns));
>
> which exposes the core iteration - n(i+1) = function ( weight(i) * 
> n(i) ). If this were linear, you could use the trick with filter, 
> given in (i).  As it is, as Paul Kienzle advised you, you will have to 
> write a C++ function to do the job.

Hmmm...

Expanding the recursion:

        n_{k+1} = sqrt(T_k n_k)
                = sqrt(T_k) sqrt(n_k)
                = sqrt(T_k) sqrt(sqrt(T_{k-1} n_{k-1}))
            = ...

I get:

        n_{k+1} = n_1^{1/2^k} \prod_{i=1}^k T_i^{1/2^{k-i+1}}

which in octave syntax is:

        n1^(2^-k) prod(T.^(2.^(-k:-1)))

This doesn't get you the intermediate values though.

If you need the intermediate values, the following will work in theory:

        z = cumprod (T.^(2.^(-k:-1))) .^ (2.^(k-1:-1:0)) .* n1 .^ (2 .^ -(1:k))

Looking at z(2)

        z2 = (T1^(2^(-k))*T2^(2^(-k+1))) ^ (2^(k-2)) * n1^(2^-2)
                = (T1^(2^(-k)))^(2^(k-2)) * (T2^(2^(-k+1)))^(2^(k-2)) * 
n1^(2^-2)
                = T1^(2^-2) * T2^(2^-1) * n1^(2^-2)
                = sqrt(sqrt(T1)) sqrt(T2) sqrt(sqrt(n1))

which is just n_3, so the entire sequence is:

        n = [n1, z]

In practice loss of precision will be an issue for any sequence long 
enough
that vectorizing is worthwhile.

- Paul





-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------