mm-421
===
Subject: Re: 2D FFT for fast polynomial multiplication?>
Another question about this: how is the fourier-transformed
function > represented? Is it unnecessary to expand the range
of sampling points > (which would require an iFFT )? After
all, more and more complex > functions (polynomials of higher
degree) are supposed to require more > and more frequency
sampling points, aren't they? So if I multiply two >
fourier-transformed Polynomials of the same m sampling points
each, > wouldn't I need to evaluate new sampling points, which
might require an > iFFT? To multiply (convolve) two polynomials
(arrays) of degree N and M, the usual method is to zero-pad
both arrays to degree at least N+M-1, do two FFTs of the same
order, multiply the Fourier amplitudes, then IFFT with the
same order to get your resulting polynomial of degree
N+M-1.(If N > M or vice-versa, you can save some time by
using overlap-add or similar tricks to express as a sequence
of convolutions of order 2*min(M,N).)(Polynomial
multiplication is just a linear convolution, for which there
are zillions of references online and elsewhere. Try Numerical
Recipes, to pick one at random.)
===Subject: Re: 2D FFT for fast
polynomial multiplication?> To multiply (convolve) two
polynomials (arrays) of degree N and M, the> usual method is
to zero-pad both arrays to degree at least N+M-1, do two> FFTs
of the same order, multiply the Fourier amplitudes, then IFFT
with> the same order to get your resulting polynomial of
degree N+M-1.But don't you need longer zero padding parts if
you multiply more than twofactors?
===Subject: Re: 2D FFT for
fast polynomial multiplication?>To multiply (convolve) two
polynomials (arrays) of degree N and M, the>usual method is to
zero-pad both arrays to degree at least N+M-1, do two>FFTs of
the same order, multiply the Fourier amplitudes, then IFFT
with>the same order to get your resulting polynomial of degree
N+M-1.> But don't you need longer zero padding parts if you
multiply more than two> factors?Yes, of course. That's what he
said inAlways pad to n, the size of the resulting polynomial.
===Subject: Re: 2D FFT for fast polynomial multiplication?>
Another question about this: how is the fourier-transformed
function > represented? Is it unnecessary to expand the range
of sampling points > (which would require an iFFT )? After
all, more and more complex > functions (polynomials of higher
degree) are supposed to require more > and more frequency
sampling points, aren't they? So if I multiply two >
fourier-transformed Polynomials of the same m sampling points
each, > wouldn't I need to evaluate new sampling points, which
might require > an iFFT?> To multiply (convolve) two
polynomials (arrays) of degree N and M, the > usual method is
to zero-pad both arrays to degree at least N+M-1, do two >
FFTs of the same order, multiply the Fourier amplitudes, then
IFFT with > the same order to get your resulting polynomial of
degree N+M-1.Ok. So this rules out the O(n) multiplication of
fourier-transformed representations because for padding, the
polynomial must be in power-coefficient-form (or however you'd
call it), and the transform costs O(n log n)> (If N > M or
vice-versa, you can save some time by using overlap-add or >
similar tricks to express as a sequence of convolutions of
order > 2*min(M,N).)But the runtime is still O(size of result
* log(size of result)), isn't it?> (Polynomial multiplication
is just a linear convolution, for which there > are zillions
of references online and elsewhere. Try Numerical Recipes, >
to pick one at random.)Great, I'll have a look at Numerical
Recipes.Best regards,
===Subject: Re: 2D FFT for fast
polynomial multiplication?> (If N > M or vice-versa, you can
save some time by using overlap-add > or similar tricks to
express as a sequence of convolutions of order > 2*min(M,N).)>
But the runtime is still O(size of result * log(size of
result)), isn't it?Not quite; it's O((M+N) * log min(M,N)) in
that case. So, there are some savings, but only in the
log.
===Subject: Re: 2D FFT for fast polynomial
multiplication?> c) If I have not only two, but k polynomials
to be multiplied, what> is the time complexity of determining
the product? O(k n log n). > Are you sure about the factor k?
Keep in mind that n is the number of > coefficients in the
*result* - and if there are more factors to start > with, the
result is also bigger so the k is already contained in n, so >
to speak.Right, sorry; I was thinking of cyclic convolutions
(i.e. a fixed degree n, mod z^n - 1). This also affects the
answer to question (b), which is basically subsumed by (or at
least similar to) your question (c).Suppose you have k
polynomials of degree m, with n = km.If you transform them all
at once to multiply together in the Fourier domain, you need to
zero-pad them all to size n, so it is O(k n log n).If you
multiply them pairwise in the obvious way, you get degree 2m,
then degree 3m, then degree 4m, and so on until degree n. This
has complexity: O(2m log 2m + 3m log 3m + ... + km log km) =
O(m log m * (2 + 3 + ... + k)) + O(m * (2 log 2 + 3 log 3 +
... + k log k)) = O(k n log (n/k)) + O(k n log k)I didn't
bother working out the second term in detail, so that term is
only an upper bound (as O(..) technically means anyway).
Besides, you can use overlap-add tricks to do the convolution
of a degree-m polynomial with a degree-Nm polynomial in O(Nm
log m) time, and this gets rid of the second term entirely.On
the other hand, if you multiply them in equal-degree pairs in
a tree with depth lg k, the complexity is: O(k/2 * 2m log 2m +
k/4 * 4m log 4m + ... + 1 * km log km) = O(km * (log 2m + log
4m + ... + log km)) = O(km * ([log m] * [lg k] + [log 2] * [1
+ 2 + 3 + ... lg k]) = O(km * ([log m] * [lg k] + [log 2] *
[lg k]^2) = O(km * (log m * k) * (log k)) = O(n (log n) (log
k))which is better than the previous approach. I'm just
working this out on the fly, so hopefully I'm not making any
algebra errors, but you should be able to work it out
yourself> Or, which would suffice: > which literature could I
use as reference to justify a statement of the > product of
multiple multivariate Polynomials can be determined in O(n >
log n) time, where n is the size of the result?theory behind
the fast-multiplication-via-FFT. There are lots of references
on the web for this, so you should read up on it. Once you
understand it, the multivariate case is the same: it's just a
multivariate convolution, and the proof of the convolution
theorem is the same (just change the indices into vectors of
indices).People use the FFT for 2d convolutions all the time
in image processing, I'm sure you can find some links there.>
Err..the above thing *was* the third question, wasn't it? Why
should it > be redundant?(You had three question marks in
question (c). The third question mark was redundant with the
first.)
===Subject: Re: learning how to use Matlab> Hello
matlab users. I am trying to figure out how to write a Matlab>
program to evaluate the function y=ln(1/1-x)> for any user
specified value of x, where x is a number <1(note that ln> is
the natural logarithm to the base e) Use an if structure to>
varify that the value passed to the program is legal. If the
value of> x is legal, calculate y(x). If not wright a suitable
error message and> quit.> Please give me some pointers on how
to write this program or even> better if you can solve it
email me, it would be better. Try something like this (I
didn't check it):function y = f(x)if (x >= 1.0) disp('bad
argument: x >= 1.0'); return;endy = log(1.0/(1.0 -
x);end
===Subject: How to calculate the total coverage area of
a few circles?In my simulation, N circles with the same radius
r are randomlyplaced. Let P_i denote the center of circle i.
For any i, p_i lieswithin the coverage range of at least one
other clicle, i.e. at leastone other circle contains p_i. How
to calculate the total coveragearea of the N overlaped circle?
The method should be easy to beimplemented by programming for
simulation.Any comments is welcome. Thank you very much.Leng
Supeng
===Subject: Pointers to Solve Modified Transport Problem
using Matlab I was taking a look at Vivekanand's problem
(message copied below)and came across the answer provided by
Erwin (copied below). Wascurious to know whether we could
solve the problem simply byconsidering a bounded
transportation problem. If instead of having 1 source with
tiered costs for transportation,i.e. X monetary units for the
first 1000 units, Y monetary units forthe next 1000 units,
couldn't we formulate it as following: We split the source
into 2 seperate sources A and B. Consider thecost of source A
as X monetary units with a 1000 units upperbound. Andthe cost
of source B as Y monetary units with a 1000 units
lowerboundand a 2000units upperbound, and so forth. This would
translate the original transportation problem intoanother
bounded transportation problem, and we could solve it using
aregular bounded transportation problem algorithm. Wouldn't
this work? Erwin mentioned non-convexities, can anyone
plzexplain. And wouldn't the method just described take that
intoconsideration and simplify the
problem?-Auro--------------- Original Message
---------------------I believe this will lead to
non-convexities. You couldformulate this as a mixed-integer
programming problemand it could even be formulated using
so-calledSOS2 variables. For a small example
see:http://www.gams.com/~erwin/interpolation.pdf--------------
--------------------------------------------------erwin@
gams.com,
http://www.gams.com/~erwin------------------------------------
----------------------------> I am trying to solve a typical
transportation problem with a small> twist. The cost of
transportation is not a fixed value but has a> tiered values.
Which means that instead of a fixed cost of> transporting, the
cost of transporting first 1000 units is some value> and the
cost of transporting next 1000 is another value less than the>
previous value (due to discounts and others). > Help me
understand if there is any standard algorithm to solve this
in> one pass. I can think of an algorithm that takes multiple
passes for> each tier and solves this.
===Subject: Counter
Intuitive Results: Optimising an FFT routine for Speedme some
good leads.I present justification for a lot of the comments
that drew(constructive) criticism below. Firstly, let me
summarise the feedbackthat has pointed me in the right
direction :-> Steven> 20.13: What's the best way of making my
program efficient?> A: By picking good algorithms and
implementing them carefully. Count me in for that,
there was no intention or suggestionof micro-optimisation or
use of inline assembler in my post. What I amtrying to find
out is what algorithm is effective and what C syntax ismore
efficient, albeit only on the small Turbo-C development
systemthat I have. The biggest speed gains in an FFT (and
most algorithms) come from much higher-> level
transformations. For example, our FFTW (www.fftw.org) ...I see
that I previously downloaded FFTW 3.0.1 is the latest official
version of FFTWa while back, I am overwhelmed with volume :458
*.c files (2.9MB)42 *.h filesCould take me a while to
assimilate all of this.> you should use a different routine
than clock()I thought I did quite well on the Athlon /Win98,
synchronising on the20ms ticks and counting in between I could
improve resolution to600ns. Still dependent on the system time
sharing on the CPU though.> The best resolution comes from
reading the CPU cycle counter directly (see >
www.fftw.org/download.html for code to do this on many
CPUs).Yes, I knew there was something called the software
stopwatch,hopefully I will come right with CYCLE.H.I can't go
and peep at this URL at the mo' as my C-FAQ download isstill
stepping through each answer page. CYCLE.H has quite a bit
ofcode to it, I suppose I just chuck out most of the CPU
specific codeand keep_( Pentium cycle counter)typedef unsigned
long long ticks;static __inline__ ticks getticks(void) { ticks
ret; __asm__ __volatile__(rdtsc: =A (ret));return ret; }>
TristanI followed my usual method of QUERY GOOGLE -> READ ->
FILTER -> POSTand using this I selected the four newsgroups
that had hits on C Speed Optimisation FFT>
comp.lang.oberon
===Subject: Re: Optimizing array bounds
checking> comp.lang.functional
===Subject: Re: Academia isn't
interested.......But what you *can* do is generate low-level
code from high-levelspecifications. This has already been
well-demonstrated in e.g. theFFTW(written in OCaml), which
generates extremely efficient C-algorithmsfor the FFT from
specifications and specialises them for the givenarchitecture.
This, however, is high-level programming again, notlow-level
programming. There is no escape from imperativeness if youwork
with von-Neumann architectures on a low level.> Eric> You'd
have to examine the actual generated code to have a hope> of
figuring out what's really happening in any particular
case.Don't want to do that, anyway my solution has to be
fairly portable,from me (design engineer) it goes to software
engineer, rebuilt underVisual C++, and then targetted to
TMS320 (or equiv) DSP.> Steven> You realize, of course, that
your first mistake was starting with the> radix-2 Numerical
Recipes in C code and spending your time attempting>
micro-optimizations. The biggest speed gains in an FFT (and
most> algorithms) come from much higher-level
transformations.Actually I started with someone else's
Fortran->C implementation of analgorithm that looks like it
came from the Cooley-Tukey paper. Then Istepped to the
Numerical Recipes (C) algorithm (Danielson - Lanczos?)(which
by the way still carries the legacy of Fortran array indexing)
This took -72% off the time.I have refined this in about five
major steps, so that at the point ofposting I was -23% off the
time of the NumRep routine as published (only highlevel
tweaks).I will perform my last trial of converting the
original code to doublewithout changing the indexing, just to
see what I would have got,before ...Now, however, I will trawl
through all the source in FFTW and see whatI can find.> If you
just want a fast FFT and are not doing this as a learning
experience, you'll> be much better off downloading an existing
optimized code.Even though I will probably do better with an
FFTW routine, thelessons learned will be valuable for other
pieces of time criticalcode.> A decent compiler should
transform array references a[j] in a loop> into separate
incremented pointers for you, if it's advantageous to do> so.
By being clever you can easily just confuse the compiler's>
optimizer unless you know what you're doing.The two sets of
indices do not change sequentially, which is why Ihoped that
me forcing the use of pointers would be faster.> Kurt> I've
done the above many a time. The only way I've been able to
cope is> by learning the instruction set and reading the
(equivalent) assembly code.> It's pretty easy to tell when
you've confused the optimizer.This would be interesting given
unlimited time, I do know that the DSPcode was definitely hand
tweaked in the past as portability andmaintainability were far
lower priority than performance.Paul
===Subject: Re: Counter
Intuitive Results: Optimising an FFT routine for
Speedpostmaster@pwi.mailshell.com says...> The best resolution
comes from reading the CPU cycle counter directly (see >
www.fftw.org/download.html for code to do this on many CPUs).>
Yes, I knew there was something called the software stopwatch,>
hopefully I will come right with CYCLE.H.> I can't go and peep
at this URL at the mo' as my C-FAQ download is> still stepping
through each answer page. CYCLE.H has quite a bit of> code to
it, I suppose I just chuck out most of the CPU specific code>
and keep> _( Pentium cycle counter)> typedef unsigned long
long ticks;> static __inline__ ticks getticks(void) {> ticks
ret;> __asm__ __volatile__(rdtsc: =A (ret));> return ret;
}This, and most of this thread is OT here, but since you've
alreadyposted this, I don't want others to be bit by this in
the future.Danger: The above is TOTALLY UNRELIABLE on SMP
systems. Gettingdependent upon this will bite you, as even
notebook computersare now available with dual CPUs. There is
no synchronizationbetween the TSCs on the individual CPUs. So,
you are as likelyas not to get bogus time values if your
process gets bouncedbetween processors.
===
===Subject: Re:
Counter Intuitive Results: Optimising an FFT routine for
Speed[...]> http://www.eskimo.com/~scs/C-faq/top.html> through
all the URLs at this site.The other versions link,
,points you to
a compressed plain-text copy of the whole thing
at. You can
uncompress itwith uncompress or gunzip. Apparently it's the
most up-to-dateversion available.
===Keith Thompson
(The_Other_Keith) kst@cts.com San
Diego Supercomputer Center <*>
Schroedinger does Shakespeare: To be
*and* not to be
===Subject: Re: Counter Intuitive Results:
Optimising an FFT routine for Speed> Don't want to do that,
anyway my solution has to be fairly portable,> from me (design
engineer) it goes to software engineer, rebuilt under> Visual
C++, and then targetted to TMS320 (or equiv) DSP.Whoa, whoa,
hold on, you want this to be fast on a DSP chip? But you're
planning to benchmark and optimize it on a general-purpose
CPU?Sigh...
===Subject: Re: Counter Intuitive Results:
Optimising an FFT routine for SpeedNNTP-Proxy-Relay:
library2.airnews.net> tempr=wr*data[ixj]-wi*data[ixj+1];>
tempi=wr*data[ixj+1]+wi*data[ixj];>
data[ixj]=data[ixData]-tempr;>
data[ixj+1]=data[ixData+1]-tempi;> data[ixData] += tempr;>
data[ixData+1] += tempi;It has been quite a few years since I
last looked at FFT code, and that wasa casual look during some
self-education time. This is the fundamentalradix-2 butterfly,
isn't it?
===Subject: Re: Counter Intuitive Results: Optimising
an FFT routine for Speed> tempr=wr*data[ixj]-wi*data[ixj+1];>
tempi=wr*data[ixj+1]+wi*data[ixj];>
data[ixj]=data[ixData]-tempr;>
data[ixj+1]=data[ixData+1]-tempi;> data[ixData] += tempr;>
data[ixData+1] += tempi;> It has been quite a few years since
I last looked at FFT code, and that was> a casual look during
some self-education time. This is the fundamental> radix-2
butterfly, isn't it?Looks like it, and programmed in a way
that is guaranteed to run slow by preventing
optimisations.
===Subject: Re: Counter Intuitive Results:
Optimising an FFT routine for Speed> Why should accessing
vectorised data using pointer references be so> much slower
than indexing - surely Turbo-C cannot have pointer>
de-referencing so badly wrong?x86 has a shortage of registers.
By forcing the compiler to storeintermediate results in
pointers, you might be running it out of registersand causing
thrashing of the register file.Just a guess - you'd need to
examine the assembly to see where it is reallygoing wrong.If
you need to optimise code, try not to use an old compiler for
a newtarget. It can't take advantage of any changes to the
architecture sincewhen it was written.
===Subject: Re: Counter
Intuitive Results: Optimising an FFT routine for Speed
<3f82fb30$0$574$b45e6eb0@senator-bedfellow.mit.eduDon't want
to do that, anyway my solution has to be fairly portable,>
from me (design engineer) it goes to software engineer,
rebuilt> under Visual C++, and then targetted to TMS320 (or
equiv) DSP.> Whoa, whoa, hold on, you want this to be fast on
a DSP chip? But> you're planning to benchmark and optimize it
on a general-purpose CPU?> Sigh...Optimising for the TMS320
processors (they vary a *lot* BTW) is *way* OTfor this group
since on the ones I've used there are all sorts of trickyou
can play in assembler that you can't (and the compiler I used
didnot) use.The OP should find a group or mailing list
specific to the processor hewants to use and probably obtain a
library written in assembler to makefull use of the facilities
the processor has. Also read the processormanual, ISTR the
TMS320C1x/2x/5x books had example FFT code...
===Mark
GordonPaid to be a Geek & a Senior Software DeveloperAlthough
my email address says spamtrap, it is real and I read
it.
===Subject: Counter Intuitive Results: Optimising an FFT
routine for SpeedOK - I feel somewhat exonerated.I tried
various combinations of options with the high iteration
codesegment. Indeed the pointerized version DID RUN
SIGNIFICANTLY FASTER(as in 21%) when I revert to 32 bit
floating point here instead ofdouble (for those entering the
thread here I repeat the salient code).This is strange but
credible, I understood that all internalarithmetic was carried
out with doubles, and casts to and from floatperformed as
necessary to match the variables data types. This doesnot
appear to be the case with Turbo-CNow I can put my present
trials to rest and start hunting throughFFTW.BTW - I could not
get CYCLE.h to work, Turbo-C does not support__INLINE__ or
__ASM__. Good thing I spent the time to get my****PORTABLE****
timer to work./*_1 curReal = wr *data[ ixj ] -wi *data[
ixj+1];_1 curImag = wr *data[ ixj+1 ] +wi *data[ ixj ];_1
data[ ixj] = data[ ixData ] -curReal;_1 data[ ixj+1] = data[
ixData+1] -curImag;_1 data[ ixData] += curReal;_1 data[
ixData+1] += curImag;*/#if defined(FFT_FLOAT)# define srcPtr0_
floatPtr_1__# define srcPtr1_ floatPtr_2__# define tgtPtr0_
floatPtr_3__# define tgtPtr1_ floatPtr_4__#endif srcPtr1_ =
srcPtr0_ = &( data[ixj]); ++srcPtr1_; tgtPtr1_ = tgtPtr0_ = &(
data[ixData]); ++tgtPtr1_; curReal = wr **srcPtr0_ -wi
**srcPtr1_; curImag= wr **srcPtr1_ +wi **srcPtr0_; *srcPtr0_ =
*tgtPtr0_ -curReal; *srcPtr1_ = *tgtPtr1_ -curImag; *tgtPtr0_
+= curReal; *tgtPtr1_ += curImag;# undef srcPtr0_ # undef
srcPtr1_ # undef tgtPtr0_ # undef tgtPtr1_ This is a typical
execution time with full double variablesFFT_NRC starting to
execute 5000 loopsFFT_NRC time = 923647 88(3.5s) nsHere is the
comparison between 32 bit floats, FFT_NRC_prev is the _1code
commented out above and FFT_NRC is the pointerized code
:FFT_NRC starting to execute 5000 loopsFFT_NRC time = 519175
88(3.5s) nsFFT_NRC_prev starting to execute 5000
loopsFFT_NRC_prev time = 658008 88(3.5s) nsFFT_NRC time change
:FFT_NRC_Prev = -138.833 0.124(3.5s)
[Micro]s(-21.10%)FFT_NRC_prev starting to execute 1000
loopsFFT_NRC_prev time = 654509 451(3.5s) nsFFT_NRC starting
to execute 1000 loopsFFT_NRC time = 515313 451(3.5s)
ns
===Subject: Re: Counter Intuitive Results: Optimising an FFT
routine for Speed> This is strange but credible, I understood
that all internal> arithmetic was carried out with doubles,
and casts to and from float> performed as necessary to match
the variables data types. This does> not appear to be the case
with Turbo-CThe x86 registers are extended precision (or at
least double, depending upon the mode), so your compiler
doesn't have much choice in the matter. Nor does this affect
the speed.Single precision is almost always faster than double
precision, though, at least for significant transform sizes,
simply because of cache usage.(Double precision variables also
need to be 8-byte aligned on x86 or the speed dies horribly and
unpredictably; not all compilers do this.)> BTW - I could not
get CYCLE.h to work, Turbo-C does not support> __INLINE__ or
__ASM__. Good thing I spent the time to get my>
****PORTABLE**** timer to work.You asked for something with
more resolution than clock(), so I told you the best
resolution you can get; if you're happy with clock(), why did
you ask? Unfortunately, accessing the cycle counter is
necessarily compiler/CPU-specific. (You also inexplicably
grabbed the code from cycle.h labeled as gcc-specific.)BTW,
clock() may be portable, but any assumptions about its
resolution are not....Not that any of this is relevant to the
DSP chip you eventually plan to run on.
===Subject: Re: 3-d
integration routine is neededThank you for all your
suggestions, I am going to work on these suggestionsand try to
come up with a solution. By the way my integration limits are
allscalar which makes the problem even simpler.> google but I
could not find anything. Do you have any suggestions? Any>
reliable routine which involves c/c++/fortran is good enough
for me.
===Subject: Re: 3-d integration routine is needed>
google but I could not find anything. Do you have any
suggestions? Any> reliable routine which involves
c/c++/fortran is good enough for
me.http://www.mat.univie.ac.at/~neum/software.html#
integraloffers a number of choices. [~ is a tilde]
===Subject:
Re: 3-d integration routine is needed> google but I could not
find anything. Do you have any suggestions? Any> reliable
routine which involves c/c++/fortran is good enough for me.in
case your integral is a convolution, you may apply FFTs
threetimesfor fast integrations.Axel
===Subject: Re: 3-d
integration routine is neededQuick and dirty solution:Use any
one dimensional integration routine recursively. Look at
netlibfor a reliable one.Mike> google but I could not find
anything. Do you have any suggestions? Any> reliable routine
which involves c/c++/fortran is good enough for me.
===Subject:
LSODARHas anyone had any success trying to solve two
independent sets ofODE's within a single program using
LSODAR?Reading the comments in the LSODAR I can see that with
this kind ofproblem the user must call the DSRCAR routine to
restore and saveLSODAR's common block before and after the
LSODAR call for each ODEset (ie each problem requires the
preservation of its own LSODARcommon block)I have implemented
the calls to DSRCAR in my program but could not getsensible
results. So as a test case I am trying to solve the samesimple
set of ODE's twice, and should therefore get two sets
ofidentical results. At the moment only the first set gives me
thecorrect answer. If I only solve for one or the other ODE
sets, I getsensible results. It just wont work with two
sets.The only way I have managed to get round this problem is
to set istate=1 and t =0 on every call and force lsodar to
re-initialise itselfeach time it is called. ie its now solving
many problems for time t =0 to tout, but using updated the y
array values from the provioussolution. I am then incrementing
time externally from lsodar. This isobviously not very
efficient.Does anyone have any ideas as to why I might not be
able to solve theproblem by using the DSRCAR subroutine?I am
using the 8th of August 2001 release of DLSODAR.I am currently
NOT providing a jacobian routine. This might beimportant as the
two sets of results only start to differ once thestiff method
is used (I think!).
===Subject: Re: LSODAR> Has anyone had any
success trying to solve two independent sets of> ODE's within
a single program using LSODAR?> Reading the comments in the
LSODAR I can see that with this kind of> problem the user must
call the DSRCAR routine to restore and save> LSODAR's common
block before and after the LSODAR call for each ODE> set (ie
each problem requires the preservation of its own LSODAR>
common block)> I have implemented the calls to DSRCAR in my
program but could notget> sensible results. So as a test case
I am trying to solve the same> simple set of ODE's twice, and
should therefore get two sets of> identical results. At the
moment only the first set gives me the> correct answer. If I
only solve for one or the other ODE sets, I get> sensible
results. It just wont work with two sets.> The only way I have
managed to get round this problem is to setistate> =1 and t =0
on every call and force lsodar to re-initialise itself> each
time it is called. ie its now solving many problems for time
t=> 0 to tout, but using updated the y array values from the
provious> solution. I am then incrementing time externally
from lsodar. Thisis> obviously not very efficient.> Does
anyone have any ideas as to why I might not be able to
solvethe> problem by using the DSRCAR subroutine?> I am using
the 8th of August 2001 release of DLSODAR.> I am currently NOT
providing a jacobian routine. This might be> important as the
two sets of results only start to differ once the> stiff
method is used (I think!).> TryRSAV = 0.d0ISAV = 0just before
calling SRCAR.
===
===Subject: Re: LSODAR> Has anyone had any
success trying to solve two independent sets of> ODE's within
a single program using LSODAR?> Reading the comments in the
LSODAR I can see that with this kind of> problem the user must
call the DSRCAR routine to restore and save> LSODAR's common
block before and after the LSODAR call for each ODE> set (ie
each problem requires the preservation of its own LSODAR>
common block)> I have implemented the calls to DSRCAR in my
program but could not> get> sensible results. So as a test
case I am trying to solve the same> simple set of ODE's twice,
and should therefore get two sets of> identical results. At the
moment only the first set gives me the> correct answer. If I
only solve for one or the other ODE sets, I get> sensible
results. It just wont work with two sets.> The only way I have
managed to get round this problem is to set> istate> =1 and t
=0 on every call and force lsodar to re-initialise itself>
each time it is called. ie its now solving many problems for
time t> => 0 to tout, but using updated the y array values
from the provious> solution. I am then incrementing time
externally from lsodar. This> is> obviously not very
efficient.> Does anyone have any ideas as to why I might not
be able to solve> the> problem by using the DSRCAR
subroutine?> I am using the 8th of August 2001 release of
DLSODAR.> I am currently NOT providing a jacobian routine.
This might be> important as the two sets of results only start
to differ once the> stiff method is used (I think!).Try> RSAV =
0.d0> ISAV = 0> just before calling SRCAR.Gerry,Could you
please expand on what you mean by this. RSAV and ISAV
arearrays. Do you mean set their contents to zero before the
first callto SRCAR? If so why ? - In my first call of SRCAR,
JOB = 1 .Anythingin ISAV and RSAV will be overwritten. Setting
them to zero at anyother time will surely result in the loss of
the common block data.Have you successfully used the SRCAR
routine before? I was assumingthe simulation wasnt working
because there is some additional commondata which this routine
does not backup.Simon
===Subject: Re: LSODAR> Has anyone had any
success trying to solve two independent setsof> ODE's within a
single program using LSODAR?> Reading the comments in the
LSODAR I can see that with this kindof> problem the user must
call the DSRCAR routine to restore andsave> LSODAR's common
block before and after the LSODAR call for eachODE> set (ie
each problem requires the preservation of its own LSODAR>
common block)> I have implemented the calls to DSRCAR in my
program but couldnot> get> sensible results. So as a test case
I am trying to solve thesame> simple set of ODE's twice, and
should therefore get two sets of> identical results. At the
moment only the first set gives me the> correct answer. If I
only solve for one or the other ODE sets, Iget> sensible
results. It just wont work with two sets.> The only way I have
managed to get round this problem is to set> istate> =1 and t
=0 on every call and force lsodar to re-initialiseitself> each
time it is called. ie its now solving many problems fortime t>
=> 0 to tout, but using updated the y array values from
theprovious> solution. I am then incrementing time externally
from lsodar.This> is> obviously not very efficient.> Does
anyone have any ideas as to why I might not be able tosolve>
the> problem by using the DSRCAR subroutine?> I am using the
8th of August 2001 release of DLSODAR.> I am currently NOT
providing a jacobian routine. This might be> important as the
two sets of results only start to differ oncethe> stiff method
is used (I think!).Try> RSAV = 0.d0> ISAV = 0> just before
calling SRCAR.> Gerry,> Could you please expand on what you
mean by this. RSAV and ISAV are> arrays. Do you mean set their
contents to zero before the first call> to SRCAR? If so why ? -
In my first call of SRCAR, JOB = 1 .Anything> in ISAV and RSAV
will be overwritten. Setting them to zero at any> other time
will surely result in the loss of the common block data.> Have
you successfully used the SRCAR routine before? I was assuming>
the simulation wasnt working because there is some additional
common> data which this routine does not backup.Simon,My
previous post was confusing as it related to ODEPACK prior to
theJune/2001 revision which had no SRCAR and JOB and one had
to 'RESTORE,FROM RSAV AND ISAV, THE CONTENTS OF THE INTERNAL
COMMON BLOCKS USED BYLSODA. ' Now it's easier.The only advise
I can offer is that:if IOPT<>0 , reinitialize rwork and iwork,
ie,rwork = 0.d0 !F95 syntax, actually only 5-7 need to be
zeroed.iwork = 0 !F95 syntax, actually only 5-9 need to be
zeroed.If rwork/iwork are not reinitialized for each ivp you
have noassurance that they are zero and, consequently, you
might not getLSODA(R)'s default values for inputs other than
the ones you elect tooverride (eg, iwork(5)=1 for verbose
messages or iwork(6) to changethe default number of steps the
solver can do per step). If IOPT=1, avalue of zero for any of
these optional inputs will cause the defaultvalue to be
used;and for successive ivp's, reinitialize t and istate, ie,t
= 0.d0istate = 1Let me know if this was of any
help.
===
===Subject: A simple
question?Content-transfer-encoding: 8bitWhile studying
Carmichael numbers I came across this curiosity:Why is it that
for any prime P, P > 29 that P^2 = K (mod (29^2-1))where K
results always in one of (1, 11^2, 13^2, 17^2, 19^2,
23^2).This is true for all tested primes, but not for all
composites. It appears that aprox. 13% of carmichael numbers
fail on this simpletest up to
10^12.P^
===#############Imagination is more important than
knowledge - A. Einstein
===Subject: Re: A simple question?>
While studying Carmichael numbers I came across this
curiosity:> Why is it that for any prime P, P > 29 that P^2 =
K (mod (29^2-1))> where K results always in one of (1, 11^2,
13^2, 17^2, 19^2, 23^2).> This is true for all tested primes,
but not for all composites. > It appears that aprox. 13% of
carmichael numbers fail on this simple> test up to 10^12.I'm
not sure why this is in sci.math.num-analysis, but:The same
thing happens for any number that has no commonfactor with
29^2-1. You can verify this by checking justthe numbers from 0
to 29^2-2. Of course any prime > 29satisfies that condition
:-). So it appears that approximately13% of Carmichael numbers
up to 10^12 have 2, 3, 5, or 7as factors.So why are there so
few squares mod 29^2-1? Well, it's2^3.3.5.7; mod those factors
there are, respectively,1,1,2,3 nonzero squares. So there are
only 6 squaresmod 29^2-1 (other than the non-coprime ones),
and it'snot all that surprising that they're represented by1
and 5 squares of small primes.
===Gareth McCaughan.sig under
construc
===Subject: who to invert 3*3 singular matrix by
support1.mathforum.org (8.11.6/8.11.6/The Math Forum,
$Revision: 1.9 primary) id h973a5P15817;I will be very
appreciated if you could invert the following 3*3matrix;fist
row entries are (2.25,.433,0) 2nd row entries are (.433,.75,0)
3rd entries are (0,0,0)thanks alot
===Subject: Re: who to invert
3*3 singular matrixX-Enigmail-Version:
0.76.1.0X-Enigmail-Supports: pgp-inline, pgp-mime> I will be
very appreciated if you could invert the following 3*3>
matrix;> fist row entries are (2.25,.433,0)> 2nd row entries
are (.433,.75,0)> 3rd entries are (0,0,0)alotThis one is a
possibility, the Moore-Penrose pseudo-inverse:(0.5, -0.2889,
0)(-0.2889, 1.5, 0)(0, 0, 0) Dan
===Subject: Re: who to invert
3*3 singular matrix> I will be very appreciated if you could
invert the following 3*3> matrix;> fist row entries are
(2.25,.433,0)> 2nd row entries are (.433,.75,0)> 3rd entries
are (0,0,0)alotWhy do you think it is called singular?
===
^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/
===
Subject: Re: who to invert 3*3 singular matrixI will be very
appreciated if you could invert the following 3*3 >matrix;
>fist row entries are (2.25,.433,0) > 2nd row entries are
(.433,.75,0) > 3rd entries are (0,0,0) >thanks alot >never try
the impossible. look in your textbook.
===Subject: Need Your
Expertise Please! by support1.mathforum.org (8.11.6/8.11.6/The
Math Forum, $Revision: 1.9 primary) id h97LQCi23500;Hi
Everyone!I have an odd request, here goes...I purchased an
Expedition with the keypad entry - no one knows whatthe
combination is and I don't have time or money to take it to
thedealer...I'd like to try combinations at my leisure until I
find the one thatworks. What I need is to know every possible
combination using thenumbers 1,3,5,7,9. Numbers may repeat, so
that needs to be taken intoconsideration as well...My email is
bblocksmith@earthlink.netSincerely,Jeanne
===Subject: Re: Need
Your Expertise Please!> Hi Everyone!> I have an odd request,
here goes...> I purchased an Expedition with the keypad entry
- no one knows what> the combination is and I don't have time
or money to take it to the> dealer...> I'd like to try
combinations at my leisure until I find the one that> works.
What I need is to know every possible combination using the>
numbers 1,3,5,7,9. Numbers may repeat, so that needs to be
taken into> consideration as well...> My email is
bblocksmith@earthlink.net> Sincerely,> JeanneHow many numbers
long is the combination? If you have a setof 5 numbers
allowing repeats, the number of possible combosis 5^n where n
is the length of the combination sequence.However it is hard
to believe someone who can afford an Expeditioncan't afford to
have a dealer read out the combination.
===
^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/
===
Subject: Re: Open source Conjugate GradientThis solves the
normal equation AT.A.x =
AT.b./*-----------------------------------------------------!
Conjugate Gradient Method for a Sparse Linear
System!------------------------------------------------------!
Ref.: Numerical Recipes, Cambridge University Press! 1986
(chapter 2.10).!! C++ demo1 version by J-P Moreau,
Paris.!-------------------------------------------------------
!! Solve sparse linear system (size=10):!! 2x1 -x10 = 0! 2x2
-x9 -x10 = 0! 2x3 -x8 -x9 = 0! 2x4 -x7 -x6 = 0! 2x5 -x6 -x7 =
0! -x5 +2x6 = 11! -x4 -x5 +2x7 = 0! -x3 -x4 +2x8 = 0! -x2 -x3
+2x9 = 0! -x1 -x2 +2x10 = 0!!! Structure of sparse matrix A:!!
X X! X XX! X XX! X XX! XXX! XX! XX X! XX X! XX X! XX X!! Any
key to continue...!! Number of iterations: 10!! Solution
vector:!! 1.0000 3.0000 5.0000 7.0000 9.0000! 10.0000 8.0000
6.0000 4.0000 2.0000!! RSQ = 1.386462e-018!!! Product A.X:!!
0.0000 0.0000 0.0000 0.0000 0.0000! 11.0000 0.0000 0.0000
0.0000
0.0000!!-----------------------------------------------------*
/#include #include #define SIZE 11#define EPS
1e-6//global variablesdouble A[SIZE][SIZE];double B[SIZE],
B1[SIZE], X[SIZE];double RSQ;int i,j,N;void ASUB(double
*,double *);void ATSUB(double *,double *);void SPARSE(double
*B,int M,double *X,double *RSQ)
{/*-----------------------------------------------------------
-------!Solves the linear system A.x = b for the vector X of
length N,!given the right-hand vector B, and given two
subroutines,!ASUB(XIN,XOUT) and ATSUB(XIN,XOUT), which
respectively calculate!A.x and AT.x (AT for Transpose of A)
for x given as first argument,!returning the result in their
second argument.!These subroutines should take every advantage
of the sparseness!of the matrix A. On input, X should be set to
a first guess of the!desired solution (all zero components is
fine). On output, X is!the solution vector, and RSQ is the sum
of the squares of the!components of the residual vector A.x -
b. If this is not small,!then the matrix is numerically
singular and the solution represents!a least-squares
approximation.
!-------------------------------------------------------------
----*///LABEL: L1 double G[SIZE],H[SIZE],XI[SIZE],XJ[SIZE];
int IRST,ITER,J; EPS2=N*EPS*EPS; IRST=0;L1:IRST++; ASUB(X,XI);
RP=0.0; BSQ=0.0; for (J=1; J<=M; J++) { BSQ=BSQ+B[J]*B[J];
XI[J]=XI[J]-B[J]; RP=RP+XI[J]*XI[J]; } ATSUB(XI,G); for (J=1;
J<=M; J++) { G[J]=-G[J]; H[J]=G[J]; } for (ITER=1; ITER<=10*M;
ITER++) { //main loop ASUB(H,XI); ANUM=0.0; ADEN=0.0; for (J=1;
J<=M; J++) { ANUM=ANUM+G[J]*H[J]; ADEN=ADEN+XI[J]*XI[J]; } if
(ADEN==0.0) printf( Very singular matrix.n); ANUM=ANUM/ADEN;
for (J=1; J<=M; J++) { XI[J]=X[J]; X[J]+=ANUM*H[J]; }
ASUB(X,XJ); *RSQ=0.0; for (J=1; J<=M; J++) { XJ[J]=XJ[J]-B[J];
*RSQ=*RSQ+XJ[J]*XJ[J]; } if (*RSQ==RP || *RSQ<=BSQ*EPS2) {
printf( Number of iterations: %dn,ITER); return; //normal
return } if (*RSQ > RP) { for (J=1; J<=M; J++) X[J]=XI[J]; if
(IRST>=3) return; //return if too many restarts goto L1; }
RP=*RSQ; ATSUB(XJ,XI); //compute gradient for next iteration
GG=0.0; DGG=0.0; for (J=1; J<=M; J++) { GG=GG+G[J]*G[J];
DGG=DGG+(XI[J]+G[J])*XI[J]; } if (GG==0.0) { printf( Number of
iterations: %dn, ITER); return; //rare but normal return } for
(J=1; J<=M; J++) { G[J]=-XI[J]; } } //main loop printf( Too
many iterations.nn);} /*In these versions of ASUB1 and ATSUB1,
we do not take any advantage of sparseness! On the contrary,
ASUB and ATSUB take into account the sparseness of the matrix
A. N.B. For big values of N, the difference in computation
time may be important. */void ASUB1(double *X,double *V) { int
I,J; for (I=1; I<=N; I++) { V[I]=0.0; for (J=1; J<=N; J++)
V[I]=V[I]+A[I][J]*X[J]; }}void ASUB(double *X,double *V) { int
I; V[1]=A[1][1]*X[1]+A[1][10]*X[10]; for (I=2; I<6; I++)
V[I]=A[I][I]*X[I]+A[I][N-I+1]*X[N-I+1]+A[I][N-I+2]*X[N-I+2];
V[6]=A[6][6]*X[6]+A[6][5]*X[5]; for (I=7; I<11; I++)
V[I]=A[I][I]*X[I]+A[I][N-I+1]*X[N-I+1]+A[I][N-I+2]*X[N-I+2];}
void ATSUB1(double *X,double *V) { int I,J; for (I=1; I<=N;
I++) { V[I]=0.0; for (J=1; J<=N; J++) V[I]=V[I]+A[J][I]*X[J];
}}void ATSUB(double *X,double *V) { int I;
V[1]=A[1][1]*X[1]+A[10][1]*X[10]; for (I=2; I<6; I++)
V[I]=A[I][I]*X[I]+A[N-I+1][I]*X[N-I+1]+A[N-I+2][I]*X[N-I+2];
V[6]=A[6][6]*X[6]+A[5][6]*X[5]; for (I=7; I<11; I++)
V[I]=A[I][I]*X[I]+A[N-I+1][I]*X[N-I+1]+A[N-I+2][I]*X[N-I+2];}/
/ show structure of sparse matrix Avoid ShowMat() { int i,j;
printf(n); for (i=1; i<=N; i++) { printf( ); for (j=1; j<=N;
j++) if (A[i][j]!=0.0) printf(X); else printf( ); printf(n); }
printf(n);}void main() {// matrix A N=10; for (i=1; i<=N; i++)
for (j=1; j<=N; j++) A[i][j]=0.0; A[1][1]=2.0; A[1][10]=-1.0;
A[2][2]=2.0; A[2][9]=-1.0; A[2][10]=-1.0; A[3][3]=2.0;
A[3][8]=-1.0; A[3][9]=-1.0; A[4][4]=2.0; A[4][7]=-1.0;
A[4][8]=-1.0; A[5][5]=2.0; A[5][6]=-1.0; A[5][7]=-1.0;
A[6][5]=-1.0; A[6][6]=2.0; A[7][4]=-1.0; A[7][5]=-1.0;
A[7][7]=2.0; A[8][3]=-1.0; A[8][4]=-1.0; A[8][8]=2.0;
A[9][2]=-1.0; A[9][3]=-1.0; A[9][9]=2.0; A[10][1]=-1.0;
A[10][2]=-1.0; A[10][10]=2.0;// vector B for (i=1; i<=N; i++)
B[i]=0.0; B[6]=11.0;// initial guess for X for (i=1; i<=N;
i++) X[i]=0.0; printf(n Structure of sparse matrix A:n);
ShowMat(); printf( Any key to continue... ); getchar();
printf(n); SPARSE(B,N,X,&RSQ); printf(n Solution vector:nn);
for (i=1; i<6; i++) printf(%8.4f,X[i]); printf(n); for (i=6;
i<=N; i++) printf(%8.4f,X[i]); printf(n); printf(n RSQ = %enn,
RSQ);// verification A.x = B ASUB(X,B1); printf(n Product
A.X:nn); for (i=1; i<6; i++) printf(%8.4f,B1[i]); printf(n);
for (i=6; i<=N; i++) printf(%8.4f,B1[i]); printf(nn);} // end
of file tsparse1.cpp
===Subject: MATLAB code to solve MCFPs
Have been looking for the past couple of days for some robust
andefficient MATLAB routines to solve Minimum Cost Flow
Problems (of theorder of 1000+ nodes). Would really appreciate
it if somebody couldhelp me find the same.-Auro
===Subject: need
help with a basic problemNeed to maximize NPV(net present
value) given budget constaint.Amt available for period 1 is
$50 & for period 2 is $20Cashflows and NPV required each
period given below :Project NPV Cash Outflow Cash Outflow
Period1 Period21 $14.00 $12.00 $3.002 17 54 73 17 6 64 15 6 25
40 30 356 12 6 67 14 48 48 10 36 39 12 18 3Set this up as a LP
problem and solve.
===Subject: Counter Intuitive Results:
Optimising an FFT routine for Speedratio is now approaching
the noise level asymptote, though there wasat least one good
point :> compressed plain-text copy of the whole thing at>
X-Last-Modified:
February 7, 1999Oh and by the way, since I have been
cross-examined on cross-posting :>
alt.comp.lang.learn.c-c++,comp.answers,alt.answers,
news.answersis the general reader of
alt.answers,news.answersgoing to be interested in the finer
points of C?> Whoa, whoa, hold on, you want this to be fast on
a DSP chip? But you're> planning to benchmark and optimize it
on a general-purpose CPU?Absolutely right - I am doing high
level system engineeting here, Iwill pass on the algorithm
implementation and source code (probablyFFTW derived, but may
still be Numerical Recipes version enhanced onseveral counts)
to the S/W engineer, he implements a version in VisualC++ for
the client (PC) workstation, and the embedded version
getsported to whatever DSP is chosen for the new hardware
platform. TheDSP speed is the most critical, but speed on the
PC client is alsoimportant. The lessons learned probably
remain with me.> The OP should find a group or mailing list
specific to the processor .....The reality of the situation is
that this won't happen.> The x86 registers are extended
precision (or at least double, depending> upon the mode), so
your compiler doesn't have much choice in the matter.> Nor
does this affect the speed.Surely the implication of this is
that float variable arithmeticshould be slower than double
arithmetic, because of the extra casts?> Single precision is
almost always faster than double precision, though,> at least
for significant transform sizes, simply because of cache
usage.should be the same speed as for float (there should be
cast overheadsthough to slow the float version down).> You
asked for something with more resolution than clock(), so I
told you> the best resolution you can get; if you're happy
with clock(), why did> you ask?> You also inexplicably grabbed
the code from cycle.h labeled as gcc-specific.toolbox, when I
have something apart from Turbo-C I am sure I will useit.2.
What resolution is this, are we talking of about 1ns, 10ns or
100ns(and I do remember that you or someone else said it is
for relativecomparisons anyway, but there will be a precise
answer for a fixed CPUand clock speed)?3. I am not happy with
clock, I am distinctly unhappy with same. Notonly is 18.2ms an
age for most fast applications, it only actuallycan get a
resolution of 666ns (900MHz Athlon, Win-98)and worst case
1145[Micro]s (998MHz Celeron, Win-2k).4. GCC specific? - I
grabbed the code for /* * Pentium cycle counter */#if
(defined(__GNUC__) || defined(__ICC)) && defined(__i386__) ...
__asm__ __volatile__(rdtsc: =A (ret));I have also tried /* *
X86-64 cycle counter */#if defined(__GNUC__) &&
defined(__x86_64__) &&!defined(HAVE_TICK_COUNTER) asm
volatile(rdtsc : =a (a), =d (d)); This at least generates the
compiler error that in-line assembly isnot allowed, the
__ASM__ keyword not being recognised at all in thePentium
example.Any chance of being able to get the assembled hex
codes for theseassembler statements and patching these into
the executable. Iimagine inserting a unique sequence of
straight C statements that asearch and replace executable
could replace with the opcodes for rdtsc: =A (ret)Seems
bizarre I know, but pragmatism was never thus
handicapped.__________________________________________________
__________________Another great Dane has made freeWith a
question of Be or Not be.Now might Schr.9adinger's puss,In
descending by Schuss,Leave one track on each side of a
tree?_________________________________________________________
___________Quantum mechanics, hmmm. You put a cat in a box,
along with a hammerandsome poison and a radioactive isotope
... I forget exactly how thisgoes.Anyway, keep some bandages
on hand, because I guarantee the cat
won'tbehappy._________________________________________________
___________________
===Subject: Re: Counter Intuitive Results:
Optimising an FFT routine for Speed> ratio is now approaching
the noise level asymptote, though there was> at least one good
point :> compressed plain-text copy of the whole thing at>
Well its plain text, and
it is contiguous, which is a lot better than> my 331 html
files totalling 982kB, but up to date it is not :>
X-Last-Modified: February 7, 1999I said it seems to be the
most up-to-date-version available; as far asI can tell, it
is.
===Keith Thompson (The_Other_Keith) kst@cts.com
San Diego Supercomputer Center <*>
Schroedinger does Shakespeare: To be
*and* not to be
===Subject: Re: Counter Intuitive Results:
Optimising an FFT routine for Speed>Whoa, whoa, hold on, you
want this to be fast on a DSP chip? But you're>planning to
benchmark and optimize it on a general-purpose CPU?>
Absolutely right - I am doing high level system engineeting
hereThe point is that optimizing for a DSP chip (with all
sorts of special-purpose instructions etc.) is completely
different from optimizing for a general-purpose CPU, and the
same code, or even the same *style* of code, is unlikely to be
fast on both kinds of architecture.> Surely the implication of
this is that float variable arithmetic> should be slower than
double arithmetic, because of the extra casts?There isn't a
one-to-one mapping from C operators to hardware operations.
The cast from float to double (or actually, to 80-bit extended
precision) is done costlessly in hardware as part of the
load/store.> 4. GCC specific? - I grabbed the code for > #if
(defined(__GNUC__) || defined(__ICC)) && defined(__i386__)
...Nota bene: __GNUC__. (If you scroll down in the file,
you'll see alternate code for Visual C++.) Anyway, I think
that this is the least of your worries.Good luck to you. I
hereby declare this thread dead.
===Subject: Make movie out of
sequential print screens?I have a complicated display that I
can save as an image usingctrl+printscrn. The problem is that
the display updates over a longperiod of time and I would like
to save each frame without having topersonally perform
ctrl+printscrn. After each frame is saved, I wouldlike to
stitch them together to make a AVI or MPEG file.Is there a way
do this?
===Subject: Re: Make movie out of sequential print
screens?This question is way off topic.The easiest format to
get to might be Motion JPEG.>I have a complicated display that
I can save as an image using>ctrl+printscrn. The problem is
that the display updates over a long>period of time and I
would like to save each frame without having to>personally
perform ctrl+printscrn. After each frame is saved, I
would>like to stitch them together to make a AVI or MPEG
file.>Is there a way do this?
===Subject: Numerical
optimization problem: smallest hypersphere?Hi all,given a
euclidian vector space V (dimension < 10) and a (possibly
ratherlarge) set of vectors x_i within V, I am looking for the
center y andradius r of the smallest hypersphere containing
these points.Therefore, I need to calculate r = min_y max_i |
x_i - y |efficiently, as I need to do many of these
optimizations.As far as I can tell, the objective function
f(y) = max_i | x_i - y | isconvex and Lipschitz continuous (if
I recall the definition correctly),thus it should possess a
unique minimum. However, because of the maxit is not
differentiable at the optimum. So my question is:
whichalgorithm for unconstrained optimization should I
consider?Playing around with a Nelder-Mead type optimizer gave
somewhatdiscouraging results. The tested algorithm proved to be
highlysensitive to the initial guess and converged to
non-optimal pointswhere it appeared to stick. This is
definitely not good. Since it isalso sometimes said that there
is no firm mathematical basis for thismethod, I am looking for
something else.Can anyone point to algorithms better suited to
the above problem?Freely available descriptions/write-ups would
be great, as well as areference implementation (any
higher-level language).
===-ha
===Subject: Re: Numerical
optimization problem: smallest hypersphere? >Hi all,given a
euclidian vector space V (dimension < 10) and a (possibly
rather >large) set of vectors x_i within V, I am looking for
the center y and >radius r of the smallest hypersphere
containing these points.Therefore, I need to calculate r =
min_y max_i | x_i - y |efficiently, as I need to do many of
these optimizations.As far as I can tell, the objective
function f(y) = max_i | x_i - y | is >convex and Lipschitz
continuous (if I recall the definition correctly), >thus it
should possess a unique minimum. However, because of the max
>it is not differentiable at the optimum. So my question is:
which >algorithm for unconstrained optimization should I
consider?Playing around with a Nelder-Mead type optimizer gave
somewhat >discouraging results. The tested algorithm proved to
be highly >sensitive to the initial guess and converged to
non-optimal points >where it appeared to stick. This is
definitely not good. Since it is >also sometimes said that
there is no firm mathematical basis for this >method, I am
looking for something else.Can anyone point to algorithms
better suited to the above problem? >Freely available
descriptions/write-ups would be great, as well as a >reference
implementation (any higher-level language). >-haminimizing the
max is the same as minimizing the max squared.name the value
of the max M. then norm(x_i - y )^2 <= M for all i and M is
the smallest such number hence minimize the function f(y,M)
def = M subject to the constraints (x_i-y)'(x_i-y) -M <= 0 for
all i happily enough this is a convex optimization problem with
a linear objective function and quadratic constraints , the
number of unkowns being quite modest (dim(y)+M ). since you
can provide also a reasonable initial and feasible guess
y0=center of mass M=the maximum radius any reasonable
nonlinear optimization code should do. for problems of this
type second order cone programming codes should be good. see
http://plato.la.asu.edu/topics/problems/nlores.html for codes
if for some reason you do not want to use one of these
packages then here even a simple augmented Lagrangian code
should do quite well since the problem has such a nice
structure
===Subject: SVD of a large matrix or different
solution.X-ID:
SAj69yZ1oeg5jGI1U3aYlua3yoMv+tEQP-8U+0I8LnTwpMJZz1m5wT Im
trying to fit a 4-6 dimensional spline to scattered
datapoints.So, essentially, I have to solve Ax=b.The
corresponding matrix A gets quite big, say, 10kx70k , with
about1k*70k of non zero entries. Since the matrix is most
probably rank deficient, I think a SVD is theright way to go.
Am I correct?The matrix is quite large, so a out-of-core
Algorithm seems in order. Does anyone have any pointers to a
library capable of decomposing sucha matrix?Or some
references/suggestions? ----Jan C. Bernauer
===Subject: Re: SVD
of a large matrix or different solution. Im trying to fit a 4-6
dimensional spline to scattered datapoints. >So, essentially, I
have to solve Ax=b. >The corresponding matrix A gets quite big,
say, 10kx70k , with about >1k*70k of non zero entries. Since
the matrix is most probably rank deficient, I think a SVD is
the >right way to go. Am I correct?The matrix is quite large,
so a out-of-core Algorithm seems in order. Does anyone have
any pointers to a library capable of decomposing such >a
matrix? Or some references/suggestions? >---- > Jan C.
Bernauerfor sure you have a 70kx10k problem , i.e. a linear
least squares problemwith 10k unknowns and a sparse matrix.
or, so to say f(x)=norm(Ax-b,2)^2 a quadratic function of 10k
variables to be minimized, with cheap operations to compute
the gradient and Hessian times a vector. no problem. see
http://plato.la.asu.edu/topics/problems/nlolsq.html for ready
to use codes. you may also have a look at the approximation
section of our guide which has pointers to scattered data
fits.
===Subject: Re: SVD of a large matrix or different
solution.X-ID:
Ek1jf-ZTQe8ktlQhUjyAHK6qNOtyuUj9wZ3VGGCDwS-gIA++
OHPQYOnospamspellucci@fb04373.mathematik.tu-darmstadt.de (
Spellucci)>for sure you have a 70kx10k problem , i.e. a linear
least squares problem>with 10k unknowns and a sparse matrix.
or, so to say> f(x)=norm(Ax-b,2)^2 Exactly!> a quadratic
function of 10k variables to be minimized, with cheap >
operations to compute the gradient and Hessian times a vector.
> no problem.> see>
http://plato.la.asu.edu/topics/problems/nlolsq.html> for ready
to use codes. you may also have a look at the> approximation
section of our guide which has pointers to scattered data >
fits.Greetings,----Jan C. Bernauer
===Subject: Re: SVD of a
large matrix or different solution.> Im trying to fit a 4-6
dimensional spline to scattered datapoints.> So, essentially,
I have to solve Ax=b.> The corresponding matrix A gets quite
big, say, 10kx70k , with about> 1k*70k of non zero entries.>
Since the matrix is most probably rank deficient, I think a
SVD is the> right way to go. Am I correct?Probably not. Since
you want to fit 70000 coefficients to 10000 data points, you
may well get a completely spurious solution.(Try fitting 7
coefficients to 1 data point!)Therefore you need to check
whatever approach you use by fitting the model to 95% of your
data and determine the quality of the fit using the remaining
5%. (This is called cross validation.)A better way to proceed
is to add equations kJx=0 where J is an approximation to
second derivative operators and k a suitable scaling constant.
Then solve (A^TA+k^2 J^TJ)x=A^Tb by conjugate gradients
(CG).For background see, e.g.
http://www.mat.univie.ac.at/~neum/papers.html#regtutorial[~ is
a tilde]For CG see any book on iterative methods.> The matrix
is quite large, so a out-of-core Algorithm seems in order.With
CG, you can generate the multiplication maps x to Ax, A^Tx, Jx,
J^Tx on thefly and need not much storage.> Or some
references/suggestions?http://www.mat.univie.ac.at/~neum/
stat.html#fithas pointers to scattered data fit
routines.
===Subject: Re: SVD of a large matrix or different
solution.X-ID:
SAjPYmZeoe7tS5Iq5io8AfwKu1xPOlnrA-i-v3VketgwIuSQ2alkcK Im
trying to fit a 4-6 dimensional spline to scattered
datapoints.> So, essentially, I have to solve Ax=b.> The
corresponding matrix A gets quite big, say, 10kx70k , with
about> 1k*70k of non zero entries.Since the matrix is most
probably rank deficient, I think a SVD is the> right way to
go. Am I correct?>Probably not. Since you want to fit 70000
coefficients to >10000 data points, you may well get a
completely spurious solution.>(Try fitting 7 coefficients to 1
data point!)Sorry, got that one in the wrong order. 70000 data
points to 10000coefficients. Somehow Im used to the screen
dimension ordering ;)>Therefore you need to check whatever
approach you use by fitting >the model to 95% of your data and
determine the quality of the fit >using the remaining 5%. (This
is called cross validation.)>A better way to proceed is to add
equations kJx=0 where J is an >approximation to second
derivative operators and k a suitable >scaling constant. Then
solve (A^TA+k^2 J^TJ)x=A^Tb by conjugate gradients (CG).>For
background see, e.g.
>http://www.mat.univie.ac.at/~neum/papers.html#regtutorial>[~
is a tilde]>For CG see any book on iterative methods.I will
look into that, thanks.> The matrix is quite large, so a
out-of-core Algorithm seems in order.>With CG, you can
generate the multiplication maps x to Ax, A^Tx, Jx, J^Tx on
the>fly and need not much storage.That would make things a lot
easier.> Or some
references/suggestions?>http://www.mat.univie.ac.at/~neum/
stat.html#fit>has pointers to scattered data fit routines.Ill
check that, too!----Jan C. Bernauer
===Subject: the definition
of frequency in fftwHi all, I am using fftw as my fourier
transform tools. I am going to takenumerical derivative with
FFT. i.e.D[u,t] = ifft(i*w*fft(u))where w is the frequency and
i is for the sqrt(-1).According to numerical recipe in C/C++, I
got the angular frequencyis of the form w = 2*pi*n/(N * dt ), n
= -N/2, ... , N/2where N is the length of the signal and dt is
the sampling interval.However, from the document of fftw, the
k-th output corresponds to thefrequency k/There T is the total
sampling period.According to the definition from fftw, I take w
= n/NIt is quite different from the definition in the previous
one (dT and2*pi has gone !!!). A testing program show that the
result is correctonly when I takeD[u, t] = ifftw( i*w*fftw(u) )
and w=n/N;
===Subject: Re: the definition of frequency in fftw
>According to numerical recipe in C/C++, I got the angular
frequency >is of the form > w = 2*pi*n/(N * dt ), n = -N/2,
... , N/2It says nothing of the kind. See the paragraph after
Eq. (12.1.8).> However, from the document of fftw, the k-th
output corresponds to the> frequency k/T> here T is the total
sampling period.This is the *frequency*, not the angular
frequency (which is what you need for differentiation); for
the latter you multiply by 2*pi.> According to the definition
from fftw, I take > w = n/NBesides the 2*pi factor, you are
forgetting about aliasing. You have a choice of whether to
call something a frequency of n/N or (n-N)/N or ... The best
approach for numerical differentiation (e.g. for solving
differential equations) is usually to take the angular
frequency with the smallest possible magnitude, i.e. for FFTW:
w = 2*pi*n/(N*dT), n = 0, ..., [N/2]-1, -[N/2], ..., -1where
[..] denotes the greatest integer (floor) function.(This is
the same as what the Numerical Recipes authors give, but to
actually use it for differentiation in their case you would
multiply by -i*w instead of i*w, because they use the opposite
sign convention from FFTW.)All of this follows simply by
inspection of the DFT (or rather the IDFT) definition, by the
way.
===Subject: Re: the definition of frequency in fftw>
However, from the document of fftw, the k-th output
corresponds to the> frequency k/T> here T is the total
sampling period.> This is the *frequency*, not the angular
frequency (which is what you > need for differentiation); for
the latter you multiply by 2*pi.Yes. In matlab (it also apply
fftw, right?), I take the an angularfrequency as the
differentitator. > According to the definition from fftw, I
take > w = n/N> Besides the 2*pi factor, you are forgetting
about aliasing. You have a > choice of whether to call
something a frequency of n/N or (n-N)/N or ... > The best
approach for numerical differentiation (e.g. for solving >
differential equations) is usually to take the angular
frequency with > the smallest possible magnitude, i.e. for
FFTW:> w = 2*pi*n/(N*dT), n = 0, ..., [N/2]-1, -[N/2], ...,
-1> where [..] denotes the greatest integer (floor)
function.Actually, I do define w like your way. I take N =
2^k, k is anpositive integer.and n = 0, 1, ... N/2, -N/2+1,
..., -1 w = 2*pi*n/(N*dT), (*)> (This is the same as what the
Numerical Recipes authors give, but to > actually use it for
differentiation in their case you would multiply by > -i*w
instead of i*w, because they use the opposite sign convention
from > FFTW.)It is not the same. In fftw, w = n/N. In which ,
2*pi and dT has gone!!! If I take -i*w as the differentitator
and define w as (*), theresult is quite different from what we
expect unless we forget 2*piand dT. ??????
===Subject: Re: the
definition of frequency in fftw> Yes. In matlab (it also apply
fftw, right?), I take the an angular> frequency as the
differentitator.This follows trivially from the
differentiation rule for exponentials; it's not specific to
matlab or fftw or numerical recipes. The fact that you can't
derive this for yourself indicates something important.> It is
not the same. In fftw, w = n/N. In which , 2*pi and dT has
gone> !!!Where in FFTW does it say that? Nowhere. The FFTW
manual says that the *frequency*, not the angular frequency (=
2*pi * frequency), is n/N. As for the factor of dT, that is
just a question of units, and the FFTW manual explicitly gives
the alternative formula n/T, where T is the total sampling time
(i.e. T = N*dT). The FFTW manual also explicitly mentions the
aliasing of positive and negative frequencies.You need to make
some effort (on your own) to have some understanding of what
you are doing; I get the sense that you just want to blindly
copy formulas.
===Subject: Re: Solve Benedict Webb Rubin
Equation of state> Can any body assist me in solving the BWR
equation for all its roots> analytically.> JohnWhat is the
equation? If it is a transcendental equation or a polynomialof
degree larger than 4, there is no analytic (i.e. closed-form)
solution.You then have to do it numerically. Getting ALL the
roots may be hard.
===
^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/
===
Subject: Searching for introductory statistical modeling textI
am looking for an introductory (graduate level physicist) text
thatdescribes modern statistical modeling. In particular, I
want tounderstand blurbs such as:The first technique uses
bootstrapping and is thought to be as accurateusing less data
than the second method which employs only brute force.For the
first method, a data set of size uniformly distributed between
25and 50 was generated and analyzed. The second method has a
data set sizebetween 100 and 200, which was subsequently
generated and analyze it. This process was repeated 1000
times.For variance reduction, we want the random numbers used
in the two methodsto be the same for each of the 1000
comparisons...Statistics and I have never been the best of
friends... words like bootstrapping, variance reduction, and
interference are all lost on me at the moment. I know what a
mean, covariances, kurtosis, skew, pdf, etc
is...Jamie
===Subject: correlation metricI'm looking for some
good metrics which may be derived througha correlation matrix.
I.e. I'm only interested in one singleoutput number as a
function of the correlation matrix.If the correlation matrix
is of size MxM and there is full correlation between all
elements then I would like to get say, 1.With correlation
between only two elements M-1 is a good number whileunder no
correlation M would be a nice number.Any tips ? References
would be nice. Thank you.
===Subject: Re: correlation metric>
I'm looking for some good metrics which may be derived
through> a correlation matrix. I.e. I'm only interested in one
single> output number as a function of the correlation matrix.>
If the correlation matrix is of size MxM and there is full >
correlation between all elements then I would like to get say,
1.> With correlation between only two elements M-1 is a good
number while> under no correlation M would be a nice number.>
Any tips ? References would be nice. Thank you.Try
(trace[R])^2 / trace[R^2], where R is the correlation
matrix.
===Subject: Re: correlation metric> Try (trace[R])^2 /
trace[R^2], where R is the correlation matrix.Is there a name
for this metric ? Is it used in some applications?
===Subject:
Re: correlation metric> Is there a name for this metric ?Not
that I know of. In this case I would call the result the
effective dimensionality of your set of M variables.> ... Is
it used in some applications?(trace[C])^2/trace[C^2], where C
is a double-centered covariance matrix, is Box's adjusted
degrees of freedom (to correct for non-sphericity) in
repeated-measures analysis of variance.For a frequency
distribution {f1,f2,...} over unordered categories,
(Sum[f])^2/Sum[f^2] gives the effective number of categories
over which the observations are spread.For data sets composed
of i.i.d. observations in which each case has associated with
it an arbitrary weight w that is used in calculating summary
statistics (e.g., mean[x] = Sum[w*x]/Sum[w]), the effective
sample size is (Sum[w])^2/Sum[w^2].There are almost certainly
others, but none come to mind just now.
===Subject: Re:
Percentile Estimation in a Monte Carlo Simulation> I am trying
to estimate a fairly extreme percentile (99.7) for a>
distribution obtained from Monte Carlo simulation. I am facing
a> problem because of the instability of data at the extremes
of the> distribution, as in, if the percentile values are
calculated for> different simulation runs, they seem to jump
around a lot. I am> currently using the MATLAB prctile
function to estimate the percentile> and am limited to
estimate it with 1000 iterations.> I would really appreciate
if anybody can suggest me a better algorithm> than what is
being used in the prctile function, which can provide a>
closer estimate to the true percentile.> UttamDavid Jones has
offered sound advice. If you have Excel (it has apercentile
function ... subject to comments users here may like tomake!),
you may wish to try it as an alternative as you will not
belimited to 1000 iterations (rather strange that Matlab
restricts youto this rather useless limitation).
===Subject:
Linux + MS Matlab + Wine = ?X-ZLExpiry: -1X-ZLReceiptConfirm:
NMail-To-News-Contact: postmaster@nym.alias.netIs it possible
to run the windows version of matlab under Linuxwith WINdows
Emulator? I know there is a Linux Matlab, but I only have the
Windows version.BTW doesn't WINE implement the complete
Windows API? How come some software does run on it? Is it
because of bugsin such software or because it relies on
undocumentedfeatures?
===Subject: Re: Linux + MS Matlab + Wine
= ?> Is it possible to run the windows version of matlab under
Linux> with WINdows Emulator? I know there is a Linux Matlab,
but I only> have the Windows version.> BTW doesn't WINE
implement the complete Windows API? How> come some software
does run on it? Is it because of bugs> in such software or
because it relies on undocumented> features?Have you
considered using GNU Octave, a Matlab replacement written by
JohnEaton and others?http://www.octave.org/It is free
software, i.e. available under the GNU General Public
License.Written for Linux and ported to run under Windows
using cygwin, it mightserve your purpose nicely.regards, chip
===Subject: Re: Linux + MS
Matlab + Wine = ?X-Enigmail-Version:
0.76.1.0X-Enigmail-Supports: pgp-inline, pgp-mime> Is it
possible to run the windows version of matlab under Linux>
with WINdows Emulator? I know there is a Linux Matlab, but I
only > have the Windows version.> BTW doesn't WINE implement
the complete Windows API? How > come some software does run on
it? Is it because of bugs> in such software or because it
relies on undocumented> features?I have tried to install it
under various Wine versions in a Mandrake installation and it
is the only program which has been able to freeze the computer
completely - so no I don't think it is worth trying :-)But with
windows 98 running under win4lin there is no problem using the
windows version of matlab.-Torbjrn.
===Subject: Re: Linux + MS
Matlab + Wine = ?> Is it possible to run the windows version
of matlab under Linux> with WINdows Emulator? I know there is
a Linux Matlab, but I only> have the Windows version.If you
encounter problems using Wine, try VMware.Evaulation versions
should be available from their web site.It's pretty sure that
you wouldn't have trouble with Matlab using that emulator. But
it doesn't come for free...Bernhard
===Subject: Re: Linux + MS
Matlab + Wine = ?: Is it possible to run the windows version
of matlab under Linux: with WINdows Emulator? I know there is
a Linux Matlab, but I only : have the Windows version.: BTW
doesn't WINE implement the complete Windows API? How : come
some software does run on it? Is it because of bugs: in such
software or because it relies on undocumented: features?Not
every software uses everything of the Windows API. In fact,
Software notwritten by Microsoft normally only uses only a
small part of the Api...Bye
===Uwe Bonnes
bon@elektron.ikp.physik.tu-darmstadt.deInstitut fuer
Kernphysik Schlossgartenstrasse 9 64289 Darmstadt---------
Tel. 06151 162516 -------- Fax. 06151 164321
----------
===Subject: Re: Linux + MS Matlab + Wine = ?> BTW
doesn't WINE implement the complete Windows API? How > come
some software does run on it? Is it because of bugs> in such
software or because it relies on undocumented> features?Wine
is alpha software, it's under development, so it simply isn't
finished.
===Subject: 2 axis matrix question!first off im only
in grade 9 so i wont know much, if u can provide mewith a
straight forward formula that would be excellent.Ok here's my
drill,I have an object on a 2d matrix and I only know the
position andfacing angle of that object. I want to know which
point is exactly orapproximately 5 units away from the object
going along the angle theobject is facing. the object can be
anywhere on the matrix and canface any angle, but still it's
destination must be 5 units along itsfacing angle.Please
someone help me. I need this for reprogramming a
game.
===Subject: Supreme by support1.mathforum.org
(8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id
h99Hjt410223;I'm in my second week of my first year of
university and hoping for alittle help with my analysis
homework as I am finding it quite astruggle. The question
is:Let A and B be bounded, nonempty sets of real numbers. Is
it always true that sup(AB) = supAsupB, where AB = {xy: x
belongsto A, y belongs to B}?I know that it is not true that
sup(AB) = supAsupB but I am havingdifficulty in finding a
counterexample. Not quite sure where to startlooking. I'd
appreciate any help or hints on this. Also if anyone can
recommend some good books for beginners trying tolearn
Analysis that would be great. Katherine
===Subject: Re:
Supreme> I'm in my second week of my first year of university
and hoping for a> little help with my analysis homework as I
am finding it quite a> struggle. The question is:> Let A and B
be bounded, nonempty sets of real numbers. > Is it always true
that sup(AB) = supAsupB, where AB = {xy: x belongs> to A, y
belongs to B}?> I know that it is not true that sup(AB) =
supAsupB but I am having> difficulty in finding a
counterexample. Not quite sure where to start> looking. > I'd
appreciate any help or hints on this. I'll give one hint.
sup() is all about the >= relation.Is it always true that b>=c
==> ab>=ac?> Also if anyone can recommend some good books for
beginners trying to> learn Analysis that would be great.
Clapham's Introduction to real analysis is very shortand
pretty clearly written, but it seems to be out ofprint.A book
I remember fondly is Ralph Boas's A primer of realfunctions. I
read this as a slightly precocious schoolpupil a year or two
before going to university, and lovedit, which at least
suggests that it's clearly written andaccessible to beginners.
There's been a new edition sincethen, and in any case I haven't
a copy to hand, so I can'tmake any guarantees :-).
===Gareth
McCaughan.sig under construc
===Subject: Help on 3D peak fit by
support1.mathforum.org (8.11.6/8.11.6/The Math Forum,
$Revision: 1.9 primary) id h9A1krV13760;I have ploted data
from a matrix as a 3D graph or a contour map (usingOrigin).
These plots are showing multiple peaks and I would like tofit
these peaks with 3D Gaussians to get numerical parameters of
eachpeak.Does anyone know a way and a software to fit these
peak and get theparameters?Thank You Jean-Pierre
===Subject:
Re: Help on 3D peak fit>I have ploted data from a matrix as a
3D graph or a contour map (using>Origin). These plots are
showing multiple peaks and I would like to>fit these peaks
with 3D Gaussians to get numerical parameters of
each>peak.>Does anyone know a way and a software to fit these
peak and get the>parameters?>Thank You >Jean-Pierrethe way to
go depends....if your model is f(x,y) = sum a_i
exp(-((x-x_i)^2+(y_y_i)^2)/s_i^2) then you have a sum of
radial functions and can use codes for fitting by radialbasis
functions. there are some in http://www.netlib.org/tomsif your
model is more general, namely f(x,y)= sum a_i
exp(-[x-x_i,y-y_i]A_i[x_x_i;y-y_i]) (in matlab notation)with
positive definite 2 by 2 matrices A_i, then you have a
nonlinear least squares problem in a quite general form.. you
can avoid the side condition on A_iby replacing it beforehand
by A_i =L_i L_i' with triangular matrices L_i as the unknowns,
but nevertheless there remainsa condition that the diagonal
elements of the L_i should by >= eps>0 for someeps. solving
this as a least squares problem is in principle routineand for
codes see
http://plato.la.asu.edu/topics/problems/nlolsq.htmlbut I
should warn you: with more than 3 peaks say (i.i. i=1,2,3)you
will get a very hard problem and even good codes may fail or
use a tremendous effort to get some reasonable result. you
also must beable to provide reasonable initial guesses for all
parameters.the a_i appear linear hence you have a so called
separable problem andcould in principle express the a_i as
functions of the other parameters.but this makes the problem
even stronger nonlinear hence i discourage this
===Subject: Re:
Help on 3D peak fit>I have plotted data from a matrix as a 3D
graph or a contour map (using>Origin). These plots are showing
multiple peaks and I would like to>fit these peaks with 3D
Gaussians to get numerical parameters of each>peak.>Does
anyone know a way and a software to fit these peak and get
the>parameters?> but I should warn you: with more than 3 peaks
say (i.i. i=1,2,3)> you will get a very hard problem and even
good codes may fail or use a> tremendous effort to get some
reasonable result. you also must be> able to provide
reasonable initial guesses for all parameters.This is done by
discarding temporarily all data note 'belonging' to each
peak(in turn), and fitting a Gaussian to them, Then solving a
linear least squares problem to weight the resulting
Gaussians.(Requires that you have enough data.)
===Subject: Re:
Help on 3D peak fit>I have plotted data from a matrix as a 3D
graph or a contour map (using>Origin). These plots are showing
multiple peaks and I would like to>fit these peaks with 3D
Gaussians to get numerical parameters of each>peak.You should
look for clustering or mixture modeling. Something like the
(*)EMalgorithm could solve your problem (if you have enough
data).> This is done by discarding temporarily all data note
'belonging' to each peak> (in turn), and fitting a Gaussian to
them, Then solving > a linear least squares problem to weight
the resulting Gaussians.> (Requires that you have enough
data.)Yes. It is faster (in my experience) to start by
assigning uniformprobabilities that each data belong to the
peak i, and by iterating somethingcalled the SEM algorithm. It
will give you a good start for more classicalalgorithm (like
EM).
===
===Subject: Newton-Raphson for complex argumentGiven
the classic Newton-Raphson solution:w_(n+1) = w_n -
f(w_n)/f'(w-n)where the process is repeated until the
difference betweenw_(n+1) and w_n is sufficiently small the
operation isobvious enough when w is a real number.But... what
is the test when w is complex (u+iv)?Is there any in print
reference material for numericalapplications for complex
numbers?
===Subject: Re: Newton-Raphson for complex argument>
Given the classic Newton-Raphson solution:> w_(n+1) = w_n -
f(w_n)/f'(w-n)> where the process is repeated until the
difference between> w_(n+1) and w_n is sufficiently small the
operation is> obvious enough when w is a real number.> But...
what is the test when w is complex (u+iv)?> Is there any in
print reference material for numerical> applications for
complex numbers?The test is still the same you just use
complex arithmetic andcontinue until the differences between
w_(n+1) and w_n; lets call itu+iv; norm (u^2+v^2) or abs sqrt(
norm )is sufficient small.I don't have any print reference
material. But one application isfinding zeros of polynomial
with either real or complex coefficients.If you are into that
you can download a windows baszed polynomialsolver using the
most common method like Newton, Laguerre, Halley andGraeffe to
solve these problems fromwww.hvks.com/numercial.htm
===Subject:
Re: Newton-Raphson for complex argument> Given the classic
Newton-Raphson solution:> w_(n+1) = w_n - f(w_n)/f'(w-n)>
where the process is repeated until the difference between>
w_(n+1) and w_n is sufficiently small the operation is>
obvious enough when w is a real number.> But... what is the
test when w is complex (u+iv)?|w_(n+1)-w_n| is sufficiently
small> Is there any in print reference material for numerical>
applications for complex numbers?You may wish to look at
Section 5.6 Complex zerosof my book Introduction to Numerical
Analysis Cambridge Univ. Press, Cambridge 2001.For a review of
the book, see
===Subject: Re: Newton-Raphson for complex
argument> Given the classic Newton-Raphson solution:> w_(n+1)
= w_n - f(w_n)/f'(w-n) ^ w_(n+1) = w_n - f(w_n)/f'(w_n) >
where the process is repeated until the difference between>
w_(n+1) and w_n is sufficiently small the operation is>
obvious enough when w is a real number.> But... what is the
test when w is complex (u+iv)?> Is there any in print
reference material for numerical> applications for complex
numbers?You do it the same way. However, you can run into the
same instabilitiesas for the real case, esp. near a saddle
point where f(z) -> 0.The basis of the Laguerre method for
finding the roots of polynomialsis complex Newton-Raphson.
Have a look at Acton, Numerical Methodsthat Work (American
Methematical Society).
===
^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/
===
Subject: Looking for a software for an estimation with
confidence intervalsA function a(n) is presumed to be
asymptotically equivalent to(n!)^s*A^n*n^beta*M when
n->+infinityI have data (ni,a(ni)) i=1..N and I would like to
identify s, A, beta andM.I need error estimates on these
values (estimation of these values withconfidence intervals)I
presume that several softwares can be used for that, but I
don't knowwhich one to choose (Matlab ? SAS ? SPSS ? others
?). I have no experiencewith statistical softwares.I've used
Maple LeastSquares function, applied on ln(a(n)), it works
finefor having values of s, A, beta and M; but it does not
give me confidenceintervals.All suggestions are welcome.
Franck
===Subject: Re: Looking for a software for an estimation
with confidence intervalsI've forgotten a precision: the values
of a(n) in the data are very large,and known with hundreds of
digits. Probably high precision will be needed toperform the
estimation.
===Subject: Questions on Function Approximation
(using FPGAs)Dear all,I am working on (compound) function
approximation with one inputvariable using piecewise
polynomial approximation with non-linear joints. These
approximations are implemented in hardware using Xilinx
FPGAs.Example of such functions include: f(x)=sqrt(-ln(x))
orf(x)=x*ln(x) where x = [0,1), which are used for Gaussian
noisegeneration (Box-Muller method) and Entropy
calculationrespectively.Does anyone know any other real-life
applications where compoundfunctions need to be
approximated?My second question is on the function
f(x)=sqrt(-ln(x)) over x =[0,1). This function is highly
non-linear and approaches infinityas x gets close to zero.
This requires floating pointimplementation (due to the large
polynomial coefficients, which Iwant to avoid). Are there any
transformations I am apply to thefunction to decompose it 2 or
more functions that are more linear?(Note that ln(x) is also
highly non-linear over x = [0,1))Dong-U Lee
===Subject: Re:
Questions on Function Approximation (using FPGAs)> Dear all,>
I am working on (compound) function approximation with one
input> variable using piecewise polynomial approximation with
non-linear > joints. These approximations are implemented in
hardware using Xilinx FPGAs.> Example of such functions
include: f(x)=sqrt(-ln(x)) or> f(x)=x*ln(x) where x = [0,1),
which are used for Gaussian noise> generation (Box-Muller
method) and Entropy calculation> respectively.> Does anyone
know any other real-life applications where compound>
functions need to be approximated?sqrt(1/x) is used so often
that more than one modern cpu contains built-in lookup tables
to generate a good starting point for a NR iteration.> My
second question is on the function f(x)=sqrt(-ln(x)) over x =>
[0,1). This function is highly non-linear and approaches
infinity> as x gets close to zero. This requires floating
point> implementation (due to the large polynomial
coefficients, which I> want to avoid). Are there any
transformations I am apply to the> function to decompose it 2
or more functions that are more linear?> (Note that ln(x) is
also highly non-linear over x = [0,1))My first guess would be
to look for some kind of rational approximation, even if this
does require a final division.Terje
===-
almost all programming can be
viewed as an exercise in caching
===Subject: Questions on
Function Approximation (using FPGAs)Dear all,I am working on
(compound) function approximation with one inputvariable using
piecewise polynomial approximation with non-linear joints.
These approximations are implemented in hardware using Xilinx
FPGAs.Example of such functions include: f(x)=sqrt(-ln(x))
orf(x)=x*ln(x) where x = [0,1), which are used for Gaussian
noisegeneration (Box-Muller method) and Entropy
calculationrespectively.Does anyone know any other real-life
applications where compoundfunctions need to be
approximated?My second question is on the function
f(x)=sqrt(-ln(x)) over x =[0,1). This function is highly
non-linear and approaches infinityas x gets close to zero.
This requires floating pointimplementation (due to the large
polynomial coefficients, which Iwant to avoid). Are there any
transformations I am apply to thefunction to decompose it 2 or
more functions that are more linear?(Note that ln(x) is also
highly non-linear over x = [0,1))Dong-U Lee
===Subject: Re
skyline matrix solver by support1.mathforum.org
(8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id
h9ADf1530198;You can use Spooles, among
others:http://www.netlib.org/linalg/spooles/spooles.2.2.html==
=Subject: Maybe it is a new number factorization algorithm. by
support1.mathforum.org (8.11.6/8.11.6/The Math Forum,
$Revision: 1.9 primary) id h9ADf3S30202;Hi everyone,First I
want to sorry for my language mistakes, because it is not
mynation language.It is possible that I had just found a new
number factorizationalgorithm. It is grounded by equation
system solving and nothing more.You create an equation system,
solve it and get results. For example:Your number: N:= 120;You
create equation system: Eqs:=
create_my_equation_system(N);Solve it: solve(Eqs);after my
algorithm solves the equation system it returns this
results(similar):1 * 120, 2 * 60, 3 * 40, 4 * 30, 5 * 24, 6 *
20, 8 * 15, 10 * 12. It works fine with any length of number.
But solving time is notsmall, jet :-)Before I publish this my
algorithm I want to check out if there is noany other similar
algorithm to my. In this case I want to register it.Well my
question would be to ask if you know any existing
similaralgorithm to my (or you think it is similar) please
send it, or writewhere I can find it.Mindaugas
Rukas
===Subject: roots of quadratic matrix equation? by
support1.mathforum.org (8.11.6/8.11.6/The Math Forum,
$Revision: 1.9 primary) id h9ABpO622561;Is there a formula for
the roots of the following quadratic matrixequation?QAQ - BQ +
C = 0where A, B, C and Q are all symmetric nxn matrices,
n>1?thanks
===Subject: Re: roots of quadratic matrix equation?>
Is there a formula for the roots of the following quadratic
matrix> equation?> QAQ - BQ + C = 0> where A, B, C and Q are
all symmetric nxn matrices, n>1?> thanksDo you really need a
closed formula ? you could solve the equationabove by
numerically means: if F(A) = QAQ - BQ + C thendF/dA = 2*AQ -
B. Now you can use Newton iteration: delta_A = - inv (dF/dA) *
F(A) A = A + delta_AGreetings, Uwe.
===Dr. rer. nat. Uwe Schmitt
http://www.procoders.net schmitt@procoders.net A service to
open source is a service to mankind.
===Subject: golf shaft
flex equation by support1.mathforum.org (8.11.6/8.11.6/The
Math Forum, $Revision: 1.9 primary) id h9B3OEc22958;The
reference numbers you are referring to are part of the
BrunswickFrequency Coefficient Measurement system for
measuring/anotating shaftstiffness in a golf club.They are
based on the resonant frequency of a 43.5 length of
shaft,clamped at one end (5) with a weight of 200 grams at the
tip. Thefrequency 'slope' of these shafts is 8.6 cpm per inch.
ie as the shaftgets shorter the frequency increases at that
rate. (or vice versa) TheFCM numbers, such as 2.0, 4.5, 6.0,
7.0 etc are merely related to themeasured cycles per minute
and represent (cpm-200)/10, eg. FCM 4.5 isthe 'standard' of a
43.5 shaft/club system which resonates at 245cpm. You can
arrange a simple chart with 'frequency' on the verticalaxis,
and shaft length on the horizontal axis. Spot the
relevantfrequency standards at the 43.5 length (ie the 2.0,
4.5, 6.0 's etc),and draw straight lines from these points
back down to the minimumlengths, sloping down (from right to
left) at the rate of 8.6 cpm perinch of reducing length. Thus
you can then read off the appropriatefrequencies at the
different shaft lengths, given the 'stiffness'nomenclature of
the actual club, or if you have a particular club ofknown
frequency you can read off from the chart the relevant
stiffness'standard' of that club.
===Subject: Bode Plots from
Data. - Matlab by support1.mathforum.org (8.11.6/8.11.6/The
Math Forum, $Revision: 1.9 primary) id h9BC2QM22328;I have
tons of data of an oscillation of some control laws. I
amhaving a difficult time finding a way to create a bode
plot(frequencyresponse) of the data. I need to do this in
MatlabIf any of you could lend me a hand, you can e-mail me
atfkatzenb@usermail.com or
katzenbergfl@navair.navy.milFrank
===Subject: Discrete
Chebyshev: best algorithm?I'm interested in discrete chebyshev
approximations. I haveto use linear functions, in the meaning
that I have to usey=mx+qand I am interested in using the best
available algorithm.I know that there is an algorithm from
Borrodale andlips (for the general problem) that is
implemented inthe nag library (by the way, how can I get it in
C?.. howalgorithms by Megiddo and Seidel that work inO(n) for
linear programming in small dimension. I amcurrently using M.
Hohmeyer implementation of Seidelrandomized algorithm.But what
is the best one? Is Borrodale' and lips' stillpreferable?Thank
you very muchm
===Subject: Re: Discrete Chebyshev: best
algorithm?> I'm interested in discrete chebyshev
approximations. I have> to use linear functions, in the
meaning that I have to use> y=mx+q> and I am interested in
using the best available algorithm.> I know that there is an
algorithm from Borrodale and> lips (for the general problem)
that is implemented in> the nag library (by the way, how can I
get it in C?.. how> much money?)nag_linf_fit is part of the NAG
C Library - you can buy it.(see
http://www.nag.co.uk/numeric/CL/manual/pdf/E02/e02gcc.pdf )It
seems to use a modified version of ACM algo 495, for which
theFortran source is available
athttp://www.netlib.org/toms/495If neither efficiency nor
readability is an issue for your applicationyou could try to
throw the Fortran source into f2c and see whatcomes out ...>
algorithms by Megiddo and Seidel that work in> O(n) for linear
programming in small dimension. I am> currently using M.
Hohmeyer implementation of Seidel> randomized algorithm.> But
what is the best one? Is Borrodale' and lips' still>
preferable?> Thank you very much> m
===Subject: 4th order
BsplineI have a stupid question about the 4th order Bspline.
In the paper(see link), they use summation to denote a
B-spline. But why they use a vector (4 j)~T.Can anyone help me
to understand
it?thanks.http://www.ri.cmu.edu/pub_files/pub2/wu_yu_te_1998_1
/wu_yu_te_1998_1.pdf
===Subject: Re: 4th order Bsplinethe 4th
order Bspline is denoted in equation 1 in the paper
below.thanks.> I have a stupid question about the 4th order
Bspline. In the paper(> see link), they use summation to
denote a B-spline. But why they use a> vector (4 j)~T.> Can
anyone help me to understand it?>
http://www.ri.cmu.edu/pub_files/pub2/wu_yu_te_1998_1/wu_yu_te_
1998_1.pdf
===Subject: Re: mathematical extrapolation of
measurement values> I need to find an algorithm to extrapolate
measurement values gathered> by a server I am writing. However,
I don't know any rules these values> adhere to, since I might
use that algorithm on several sets of> measurement values with
different maths lying beneath.> Best idea I have come up with
so far is to put a curve of nth potence> through the
measurement values of mine and hope that the resulting> curve
allows some extrapolation into the future. So far I have>
realized this for a 2nd degree curve like f(x)=ax*x+b*x+c
using three> points of my measurements. Doesn't work too well,
though.> Requirements and other stuff:> I know that my
measurements won't show any abrupt changes. The values> I'm
watching tend to change rather slow. Furthermore,
implementation> time and calculation time are both restricted.
I can't spend much time> for either - meaning I can spend about
1-2 days for the implementation> of the algorithm and I'll need
to be able to calculate (roughly> guessed) 1000+ extrapolations
per second on a publically available PC> (making me refrain
from iterations). Additionally, the calculation> should need
only as few measurement values as reasonably possible,> since
I don't have a history of my values and need to store my>
calculation-values separately. Extrapolation should start with
only a> hand full of values gathered. Here's an approach that
worked well in retrofitting the (classical BDF)derivative of
the PID controller fcn. Create a sliding interval withthree
fcn values, fit it with a cubic spline and then use the
endpointderivative for extrapolation.You can adapt spline from
netlib/sfmm dir for maximum efficiencyby return(ing) right
after the spline derivative (b(n), for n=3)
iscalculated.--Dr.B.Voh---------------------------------------
---------------Applied Algorithms
http://sdynamix.com
===Subject: Runge Kutta ABC's inner working
?Do you kwow where I can found a smooth introductionto Runge
Kutta methods (internet site). I know there is google but I
want your opinionabout a site for the non mathematician
layman.Thank you for sharing
===Subject: Orthogonal polynomials
, weight w(x)=exp(-k*x) , (a,b)=(0,1).Content-Length:
941Originator: rusin@vesuviusSuppose that (P_n(x))_{n>=0} ,
(P_n(1)=1 ,n=0,1,...) , is the polynomial system orthogonal on
[0,1] with respect to weight function w(x)=e^{-x} . In other
words, for all non-negative integers n,j
(P_n,P_j):=Integral_{x=0 to x=1} P_n(x)P_j(x)w(x) dx = 0 when
n=/=jand (P_n,P_n) =/=0 .I am interested to find coefficients
(A_n),(B_n) and (C_n) fromthe three-term recurrence P_{n+1}(x)
= (A_n*x+B_n)*P_n(x) - C_n*P_{n-1}(x) , n=1,2,..., with
P_0(x)=1 , P_1(x)=(e-1)*x+(2-e).The same question in a more
general case when w(x)=e^{-k*x},k in R, k=/=0 . Integrals of
form Integral_{x=0 to x=1}e^{-k*x}f(x) dx arise in quantum
mechanical computations (for instance see R. Mach ,
J.Math.Phys 25 (1984)). Of interest is also the symmetric
weight w(x)=e^{-x(1-x)} .I need the coefficients in order to
find some Gaussian quadraturesfor above weights.If possible,
please some references. Thanking in advance,
Alex.
===
===
===
===Subject: Re: Orthogonal polynomials , weight
w(x)=exp(-k*x) , (a,b)=(0,1).Content-Length: 1471Originator:
rusin@vesuvius >Suppose that (P_n(x))_{n>=0} , (P_n(1)=1
,n=0,1,...) , >is the polynomial system orthogonal on [0,1]
with >respect to weight function w(x)=e^{-x} . >In other
words, for all non-negative integers n,j >
(P_n,P_j):=Integral_{x=0 to x=1} P_n(x)P_j(x)w(x) dx = 0 when
n=/=j >and (P_n,P_n) =/=0 .I am interested to find
coefficients (A_n),(B_n) and (C_n) from >the three-term
recurrence >P_{n+1}(x) = (A_n*x+B_n)*P_n(x) - C_n*P_{n-1}(x) ,
n=1,2,..., > with P_0(x)=1 , P_1(x)=(e-1)*x+(2-e). >The same
question in a more general case when w(x)=e^{-k*x},k in R,
k=/=0 . > Integrals of form Integral_{x=0 to x=1}e^{-k*x}f(x)
dx >arise in quantum mechanical computations (for instance see
R. Mach , > J.Math.Phys 25 (1984)). >Of interest is also the
symmetric weight w(x)=e^{-x(1-x)} . >I need the coefficients
in order to find some Gaussian quadratures >for above weights.
>If possible, please some references. Thanking in advance,
Alex. >
===
===
=== >you could compute the coefficients using the
scalar product itself, see e.g. theformulae in stoer-bulirsch.
but it is not a good idea to use them for computing weights
and nodes of quadrature formulae because of problems with
roundoffwalter gautschi did a lot of work on this and his code
orthpol inhttp://www.netlib.org/toms will serve your needs
perfectly.
===Subject: Simplifying logical propositionI would
like to know if anyone could help me understanding how Maple
did thefollowing simplification (using bsimp):Proposition:
(not p or q) and ((q and not r) or (r and not q))Simplified
proposition: r and not p and not qThank you
===Subject: Re:
Simplifying logical proposition> I would like to know if
anyone could help me understanding how Maple didthe> following
simplification (using bsimp):If it does, good question since
what you post below are not equivalent. Sothose trying to
prove it can stop.> Proposition: (not p or q) and ((q and not
r) or (r and not q))> Simplified proposition: r and not p and
not qp q r ~p or q q and ~r r and ~q (q and ~r) (~p or q) or
(r and ~q)and ((q and ~r) or (r and ~q))F F F T F F FFF F T T
F T TTF T F T T F TTF T T T F F FFT F F F F F FFT F T F F T
TFT T F T T F TTT T T T F F FF> Simplified proposition: r and
not p and not qp q r r and ~p (r and ~p) and ~q (parenthesis
not reallyneeded)F F F F FF F T T TF T F F FF T T T FT F F F
FT F T F FT T F F FT T T F F
===Subject: FFT Spectral
Analysis - SpectraLAB, SpectraPLUS,FFT Spectral Analysis -
SpectraLAB, SpectraPLUS,SpectraPRO, SpectraRTA, acoustic
calculation, simulationand sound system designSpectraLAB
v4.32.16cSpectraLAB v4.32.17SpectraPLUS v2.32.03cSpectraPLUS
v2.32.04SpectraPRO v3.32.16cSpectraPRO332 v3.32.17SpectraRTA
v1.32.12cSpectraRTA132 v1.32.13for more info, please send
e-mail
===Subject: calculating log of the ratio of
determinantsI have to calculate the log ratio of two symmetric
p.d 500x500determinants S1 and S2, the determinants are too
large and will notfit into a double variable.Since I cannot
calculate the determinants exactly I can think of thefollowing
ways to do it1]calculate S2^{-1}S1, find out its eigenvalues
and add up the sum oftheir logs2]calculate sum of the logs of
the ev's of S1, find out sum of the ev'sof S2 and subtractIs
there any reason to prefer one over the another?
===Subject:
Re: calculating log of the ratio of determinants> I have to
calculate the log ratio of two symmetric p.d 500x500>
determinants S1 and S2, the determinants are too large and
will not> fit into a double variable.> Since I cannot
calculate the determinants exactly I can think of the>
following ways to do it> 1]> calculate S2^{-1}S1, find out its
eigenvalues and add up the sum of> their logs> 2]> calculate
sum of the logs of the ev's of S1, find out sum of the ev's>
of S2 and subtract> Is there any reason to prefer one over the
another?The second is usually much cheaper and is often used in
practice.you can also add the binary exponents and multiply the
mantissae(but if your matrix would be much larger you would
occasionally need to rescale the mantissa product by some
power of 2.
===Subject: Re: sign functionPlease do not post the
same question _separately_ to different newsgroups.If you feel
that you _must_ post the same question to more than one
group,crossposting is much better.I already responded to your
question in k12.ed.math . Interested readersmay see the
response there, once its moderator posts my message.David
Cantrell
===Subject: How can I get the frequency of a
periodical functionI have n datapoints
{(t1,x1),(t2,x2),...,(tn,xn)} of a periodicalSignal t->x(t).
Now I search the Fourier-Approximation for thisSignal. When I
search the approximation with fft, so I got themagnitudes ak
and bk to the harmonical functions sin(k*w0*t) andcos(k*w0*t)
(where w0=2*pi/(tn-t1) and T=tn-t1).My problem is, that I
don't know the frequency of the signal. Thedatapoints are
dristributed over a range tn-t1=k*T (k>1,real). How canI find
T?
===Subject: FFTW and polynomail MultiplicationI am trying to
use FFTW to square a polynomial of large degree N. Ihave tried
the
following-----------------------------------------------------
-----fftw_complex *in, *out;in = fftw_malloc(N* sizeof
(fftw_complex));out = fftw_malloc(N* sizeof
(fftw_complex));for(int i=0;i < N/2; i++) /* in is my
row-vector of Array coefficients */ in[i][0]= MY ARRAY
COEFFICIENTS; /* all imaginary parts are 0 */ in[i][1]= 0;}p =
fftw_plan_dft_1d(N, in, out, FFTW_FORWARD,
FFTW_MEASURE);fftw_execute(p);for(int i=0;i < N; i++)
out[i][0]= out[i][0]*out[i][0];}p = fftw_plan_dft_1d(N,out,in,
FFTW_BACKWARD,
FFTW_MEASURE);fftw_execute(p);print(in);----------------------
------------------------------------now when I print... I get
all zeroes in the real part of 'in'.what am I doing wrong? Is
there any hope for me?-Manish
===Subject: Re: FFTW and
polynomail Multiplication> I am trying to use FFTW to square a
polynomial of large degree N. I> have tried the following
[....................................]> now when I print... I
get all zeroes in the real part of 'in'.See question 3.16 in
the FFTW FAQ.
===Subject: Numerical PuzzleI have been working
on a puzzle sponsored by my university but areunable to solve
it. I am asking if someone here can solve it. If yousolve it
in time you are eligible for a PocketPC or Microsoftsoftware.
The hint is Fortune Teller Origami and the title is2914849
steps to go? Better save all my powers! The puzzle is
here:http://www.acm.uiuc.edu/puzzlehack/ph2/puzzles2/1/puzzle1
.jpg