On 9/8/10 8:39 PM, KvS wrote:
Thanks all, however I am not very successful so far :(.
I tried both options mentioned before:
- only optimize the loops in Cython and keep using symbolic
expressions/infinite precision, but this is unfortunately rather
slow;
- fully optimize in Cython by turning to doubles everywhere, although
it gets very fast here happens what I was already afraid of: the
algorithm goes belly up after 15 steps :( (I would need a lot more
steps). In step 14 values like 2.51885860e-34 appear and in
combination with the numerical noise gathered this turns out to
produce very wrong values at the next step.
@Robert: I'm trying to implement an algorithm that computes a sequence
of functions v_k according to the rule v_k(x) = \max \{ K-e^x, \int
v_{k-1}(y) p(x+y) dy \}, v_0(x)=K-e^x, where p is a prob. density. The
background is an optimal stopping problem in stochastics. So at each
step I essentially first compute the integral and then determine where
the integral and x -> K-e^x intersect, this just by the basic
find_root function in Sage.
It's not difficult (however very annoying...) to write down a formula
for the integral (given the specific density p I'm using) and work out
formulae for all the coefficients that appear in this formula for the
integral in terms of those appearing in v_{k-1}. The number of
coefficients (and hence computations) involved in the formula for v_k
increases very quickly with k, that explains why the algorithm gets
slow in 'symbolic mode'. However 'double mode' is not good enough as
mentioned above.
I will try your suggestion for using mpfr, thanks Robert and Greg!
(however it looks a bit scary to me ;)). I need to store the huge
amounts of coefficients somewhere, I guess numpy is my best guess even
though I have to declare the dtype as object then, no?
One of the things you might do is use RR (which is based on mpfr), and
when you need to do a sequence of calculations really fast, reach in to
the RR number, do stuff using mpfr, and then create a new RR number with
the new multiprecision result. The RR _add method shows how to do this,
for example:
cdef RealNumber x
x = self._new()
mpfr_add(x.value, self.value,
(<RealNumber>other).value,(<RealField_class>self._parent).rnd)
return x
Notice how it creates a new real number x, does something using mpfr to
manipulate the mpfr value, but then returns the RR number.
You can find lots more examples in devel/sage/sage/rings/real_mpfr.pyx
Thanks,
Jason
--
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
