bridge/class-master/depend/pari/src/language/intnum.c
/********************************************************************/
/**             NUMERICAL INTEGRATION (Gauss-Legendre)             **/
/********************************************************************/
NUMERICAL INTEGRATION (Gauss-Legendre)
NUMERICAL INTEGRATION ( Gauss - Legendre )
GEN
intnumgaussinit(long n, long prec)
{
  pari_sp ltop = avma;
  GEN L, dp1, p1, p2, R, W;
  long prec0 = prec + EXTRAPRECWORD;
  long bitprec = prec2nbits(prec), i, d1;
  if (n <= 0) n = (long)(bitprec*0.2258);
  if (odd(n)) n++;
  if (n == 2) n = 4;
  prec = nbits2prec(3*bitprec/2 + 32);
  L = pollegendre(n, 0); 
  p1 = Q_remove_denom(RgX_deflate(L, 2), &dp1);
  d1 = vali(dp1);
  p2 = ZX_deriv(p1); 
  R = ZX_Uspensky(p1, gen_0, 1, 3*bitprec/2 + 32); 
  n >>= 1;
  W = cgetg(n+1, t_VEC);
  for (i = 1; i <= n; ++i)
  {
    GEN t, r2 = gel(R,i);
    if (typ(r2) != t_REAL) r2 = gtofp(r2, prec);
    gel(R,i) = sqrtr_abs(r2); 
    t = mulrr(subrr(r2, sqrr(r2)), sqrr(poleval(p2, r2)));
    shiftr_inplace(t,1-2*d1);
    gel(W,i) = invr(t);
  }
  R = gprec_wtrunc(R,prec0);
  W = gprec_wtrunc(W,prec0);
  return gerepilecopy(ltop, mkvec2(R,W));
}
GEN intnumgaussinit ( long n , long prec ) { pari_sp ltop = avma ; GEN L , dp1 , p1 , p2 , R , W ; long prec0 = prec + EXTRAPRECWORD ; long bitprec = prec2nbits ( prec ) , i , d1 ; if ( n <= 0 ) n = ( long ) ( bitprec * 0.2258 ) ; if ( odd ( n ) ) n ++ ; if ( n == 2 ) n = 4 ; prec = nbits2prec ( 3 * bitprec / 2 + 32 ) ; L = pollegendre ( n , 0 ) ; p1 = Q_remove_denom ( RgX_deflate ( L , 2 ) , & dp1 ) ; d1 = vali ( dp1 ) ; p2 = ZX_deriv ( p1 ) ; R = ZX_Uspensky ( p1 , gen_0 , 1 , 3 * bitprec / 2 + 32 ) ; n >>= 1 ; W = cgetg ( n + 1 , t_VEC ) ; for ( i = 1 ; i <= n ; ++ i ) { GEN t , r2 = gel ( R , i ) ; if ( typ ( r2 ) != t_REAL ) r2 = gtofp ( r2 , prec ) ; gel ( R , i ) = sqrtr_abs ( r2 ) ; t = mulrr ( subrr ( r2 , sqrr ( r2 ) ) , sqrr ( poleval ( p2 , r2 ) ) ) ; shiftr_inplace ( t , 1 - 2 * d1 ) ; gel ( W , i ) = invr ( t ) ; } R = gprec_wtrunc ( R , prec0 ) ; W = gprec_wtrunc ( W , prec0 ) ; return gerepilecopy ( ltop , mkvec2 ( R , W ) ) ; }
NUMERICAL INTEGRATION (Gauss-Legendre)
NUMERICAL INTEGRATION ( Gauss - Legendre )
/* n even >= 4, p1 is even */ /* L_n = p1(x^2) */ /* L_n' = 2x p2(x^2) / 2^d1 */ /* positive roots of p1 */ /* positive root of L_n */ /* 2 / (L'(r)^2(1-r^2)) =  2^(2d1 - 1) / (1-r2)r2 (p2(r2))^2 */