If F is a quadratic function on a n-dimensional vector space,
then F(x) is affine equivalent to one of the following:

  Sum( epsilon[j]*x[j]^2, j=1 ..r  ), 
  Sum( epsilon[j]*x[j]^2, j=1 ..r  ) + alpha, 

This provides a Maple solution to compute the bivariate normal distribution by recursions for numerical inputs. It works even for extreme cases and handles situations, where usual integration with Maple has serious problems (even after reducing to dimension 1), it seems to be reliable and fast and works in 'arbitrary' precision.

To use it call N2_as_sum(1.0, 2.0, 0.8,  200) to compute the BVN for x = 1.0, y = 2.0 and correlation rho = 0.8 with at most 200 recursion steps (it will stop earlier, if no more improvements can be seen).

In a concurrent thread I posted the following simplification procedure

  Tryhard:= proc(expr)
    global E_in_Tryhard;
    subs(pow= `^`,
      codegen[optimize](subs(E_in_Tryhard= expr, ()-> E_in_Tryhard), tryhard))()
  end proc;

and now put into an extra blog post, as it might be helpful for others. Note, that this does not work on all constructs in Maple and certainly the package is a bit dated.

'Int(x^(k+v-1)/(1+x),x = 0 .. infinity)': '%'= value(%);

             infinity
            /          (k + v - 1)
           |          x
           |   ...

Based on a problem / bug with 'fibonacci', which Plouffe states in the symbolic newsgroup, I isolated from his example a situation, where elementary operations on lrage interger sends Maple to death on WIN (but seems to work on a MAC).

www.mapleprimes.com/files/102_fibonacci_bug.mws.zip (190 KB)