From help-request at octave dot org Sat Feb  5 15:30:34 2005
Subject: Re: octave-speed
From: Paul Kienzle <pkienzle at users dot sourceforge dot net>
To: "Paul Thomas" <paulthomas2 at wanadoo dot fr>
Cc: =?ISO-8859-1?Q?Rafael_Rodr=EDguez_Velilla?= <rrv66 at tid dot es>,         <help@octave.org>
Date: Sat, 5 Feb 2005 16:35:00 -0500

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
-------------------------------------------------------------