[R] integrate and quadratic forms

[rdporto1]

Hi all.

I'm trying to numerically invert the characteristic function
of a quadratic form following Imhof's (1961, Biometrika 48)
procedure.

The parameters are:

lambda=c(.6,.3,.1)
h=c(2,2,2)
sigma=c(0,0,0)
q=3

I've implemented Imhof's procedure two ways that, for me,
should give the same answer:

#more legible
integral1 = function(u) {
  o=(1/2)*sum(h*atan(lambda*u)+sigma^2*lambda*u/(1+lambda^2*u^2)) - q*u/2
  rho=prod((1+lambda^2*u^2)^(h/4))*exp( (1/2)*sum((sigma*lambda*u)^2/(1+lambda^2*u^2)) )
  integrand = sin(o)/(u*rho) 
}

#same as above
integral2= function(u) {
((1/2)*sum(h*atan(lambda*u)+sigma^2*lambda*u/(1+lambda^2*u^2)) - q*u/2)/
(u*(prod((1+lambda^2*u^2)^(h/4))*
exp( (1/2)*sum((sigma*lambda*u)^2/(1+lambda^2*u^2)) )))
}

The following should be near 0.18. However, nor the answers are near this
value neither they agree each other!

> 1/2+(1/pi)*integrate(integral1,0,Inf)$value
[1] 1.022537
> 1/2+(1/pi)*integrate(integral2,0,Inf)$value
[1] 1.442720

What's happening? Is this a bug or OS specific? Shouldn't they give the 
same answer? Why do I get results so different from 0.18? In time:
the procedure works fine for q=.2.

I'm running R 2.4.1 in a PC with Windows XP 32bits. Other ways (in R) to
find the distribution of general quadratic forms are welcome.

Thanks in advance.

Rogerio.