(*
 AWGQ
Abscissas and Weight of Gaussian Quadrature.
H Fukuda, M Katuya, E O Alt and A V Matveenko 
\int_a^b{f(z(x))W(x)dx}.
W(x)>0 z(x):monotonic 

2003.9.7-9 for Mathematica 2.2.1 
2004.3.22-5.20 for Mathematica 4.0.1 

This program may work 
for any a,b,z(x) and W(x),
   if moment\int_a^b{z(x)^i W(x)dx} 
  is given analytically by Mathematica.
 
 Results are stored in the array,
      xx[i],i=1,..,n:abscissas 
    ww[i],i=1,..,n:weight
 and given in FORTRAN subroutine.
*)

(* define integral *)

(* limits a and b of integration region *)
a=0;
b=1;

(* weight function *)
c=1;
W[x_]:=1/(x + 1/2);

(* argument *)
z=x;

(* inverse relation *)
zi=u;

(* degree of formula *)
n=4;

(* Internal precision *)
(* in number of digits *)
Nm=50; (* moments *)
Nz=25;   (* Jacobi matrix *)

(* final Precision *)
Na=17;

(*
   Table 1. Example of internal precision.
n= 4:Nm=   50 Nz=25     1[sec] (P4 1.7GHz);
   n= 8:Nm=  100 Nz=25     2[sec]; 
   n=16:Nm=  400 Nz=25     5[sec];
   n=32:Nm= 1300 Nz=25    20[sec];
   n=64:Nm= 6000 Nz=34  3246[sec];
   n=96:Nm=12000 Nz=25 32323[sec];
*)

(* FORTRAN output *)
outputww=0; 
    (* 
  =0: FORTRAN subroutine provides weights itself
=1: FORTRAN subroutine provides weights divided by W(x)
 *)

(* FORTRAN filename *)
fortranfname="gqxw.f";

(**********************************************)

Tstart=AbsoluteTime[];

(* Print parameters *)
Print["weight W(x)=",W[x]];
Print["argument z(x)=",z];
Print["integral region {a,b}=",{a,b}];
Print["Degree of formula"];
Print[" n=",n];
Print["Precision"];
Print[" moment Nm=",Nm];
Print[" Jacobi Nz=",Nz];
Print[" final  Na=",Na];
Print[" MachinePrecision=",$MachinePrecision];
Print["FORTRAN file"];
Print[" ",fortranfname];

Print["Evaluating moment z^i W[x]..."];

Clear[s];
Do[
s[i]=0;
,{i,0,2*n-1}
]

Do[
 i=2*n-1-j;
zW=Integrate[z^i W[x],{x,a,b}];
(* check i'th moment *)
If[
NumericQ[zW]==False,
Print[i,"'th moment is not given analytically."];
Abort[];
];
s[i]=Re[N[zW,Nm]];
Clear[zW];
If[Mod[j,Floor[(2*n+1)/10]+1]==1, 
Print["...",Floor[100*j/(2*n+1)],"%"]
];
,{j,0,2*n-1}
]

Print["Orthogonal polynomial and Jacobi matrix ..."];
Clear[aa,bb];
i=0;
p0=1;
i=i+1;
pp0=s[0];
xpp=s[1];
aa[0]=xpp/pp0;
bb[0]=0;
p1=(u-xpp/pp0)p0;
Do[
pp1=Sum[Coefficient[p1,u,j] Coefficient[p1,u,k] s[j+k],{j,0,i-1},{k,0,i-1}];
xpp=Sum[Coefficient[p1,u,j] Coefficient[p1,u,k] s[j+k+1],{j,0,i-1},{k,0,i-1}];
aa[i-1]=xpp/pp1;
bb[i-1]=pp1/pp0;
p2=Expand[(u-xpp/pp1)p1-pp1/pp0 p0,u];
p0=p1;
pp0=pp1;
p1=p2;
Clear[p2,xpp,pp1];
If[Mod[i,Floor[n/10]+1]==0, 
Print["...",Floor[100*i/n],"%"]
];
,{i,2,n}
];

jm=Table[
If[i==j,N[aa[i-1],Nz],
If[j==i+1,N[Sqrt[bb[i]],Nz],
If[i==j+1,N[Sqrt[bb[j]],Nz]
,0]]]
,{i,n},{j,n}];
jeig=Eigensystem[jm];
jz=Reverse[ jeig[[1]]];
Print["abscissas"];
Print[jz];

Print["weights"];
jv=jeig[[2]];
jw=Reverse[Table[s[0]*jv[[i,1]]^2,{i,n}]];
jws=Table[{w->jw[[j]]},{j,1,n}];
Print[jw];

jxs=Table[{x->zi /.{u->jz[[j]]}},{j,1,n}];
Do[
jx[j]=zi/.{u->jz[[j]]},{j,1,n}
];


(* FORTRAN output *)
epart[x_]:=Floor[Log[10,Abs[x]]];
dForm[x_]:=ToString[N[Re[x]/10^epart[x],Na]]<>"D"<>ToString[epart[x]];
fform[x_]:=FortranForm[N[x,Na]];
stopmsg="'"<>fortranfname<>", n=', n";

o=OpenWrite[fortranfname];
WriteString[o,"c     Abscissas and weights of Gaussian Quadrature\n"] ;
WriteString[o, "c     produced by Mathematica code: AWGQ\n"];
WriteString[o, "cccccc\n"];
WriteString[o, "c     for int_a^b{f(z(x))W(x)dx}\n"];
WriteString[o,"c     argument    z(x) = ",FortranForm[z],"\n"];
WriteString[o,"c     weight      W(x) = ",FortranForm[W[x]], "\n"];
WriteString[o,"c     lower limit a = ",FortranForm[a],"\n"];
WriteString[o,"c     upper limit b = ",FortranForm[b],"\n"];
WriteString[o,"c     GQ order n<=200\n"];
WriteString[o, "cccccc\n"];
WriteString[o,"      subroutine gqxw(x,w,n)\n"];
WriteString[o,"c     x: abscissas\n"];
If[outputww==0,
WriteString[o,"c     w: weights\n"];,
WriteString[o,"c     w: weights divided by W(x)\n"];
]
WriteString[o,"      implicit real*8(a-h,o-z)\n"];
WriteString[o,"      dimension x(200),w(200)\n"];
WriteString[o,"c\n"];
WriteString[o,"      if(n.eq.",ToString[n],") then\n"];
Do[
WriteString[o,"        x(",ToString[i],")=",dForm[jx[i]],"\n"];
If[outputww==0,
WriteString[o,"        w(",ToString[i],")=",dForm[jw[[i]]],"\n"];,
WriteString[o,"        w(",ToString[i],")=",dForm[jw[[i]]/W[jx[i]]],"\n"];
]
,{i,1,n}
];
WriteString[o,"        return\n"];
WriteString[o,"      end if\n"];
WriteString[o,"c\n"];
WriteString[o,"      stop ", stopmsg,"\n"];
WriteString[o,"      end\n"];
Close[o];

(* Check result *)
Print["Check results"]
Print["exact moment, moment by GQ, (relative error)"];
 Do[f=z^i;
gaussianm=N[Sum[f w/.jxs[[j]]/.jws[[j]],{j,1,n}], Na];
analyticm=N[s[i],Na];
rerror=Abs[ (gaussianm-analyticm)/analyticm];
Print["  m",i,": ",
analyticm,", ",
gaussianm,", (",
FortranForm[rerror],")"
];
,{i,0,2 n-1}
];
Print[
"*** computation time=",
AbsoluteTime[]-Tstart,"[sec]"
];