\\ Approximate integral(0,inf,x,f(x)*sin(c*x+q*Pi)) \\ by the Ooura-Mori transformation using original step size Pi/M {q_ossinfM(f,M,c,q)= local(h,alpha,beta,n,wsign,active,s,tn, t,dt,phi,dphi, het,et,dx,x,hex,ex,dy,y,A,w,w0,s0,term); h=Pi/M; beta=0.25; alpha=beta/sqrt(1+M*log(1+M)/(4*Pi)); n=0; wsign=vector(2); active=vector(2,k,1); s=0; tn=vector(2); while(active[1]+active[2]>0, tn[1]=-n-q; tn[2]=n+1-q; n=n+1; for (j=1, 2, if(active[j], t=h*tn[j]; dt=1; if (t==0, ex=1; phi=1/(2+alpha+beta); dphi=((alpha-beta)*phi^2+1)/2, /*else*/ het=exp(t/2); et=het*het; dx=dt*(-2-alpha/et-beta*et); x=-2*(t+(alpha/et+beta)*het*sinh(t/2)); hex=exp(x/2); ex=hex*hex; dy=-dx*ex; y=2*hex*sinh(-x/2); phi=t/y; dphi=dt/y-phi*(dy/y)); A=M*phi; x=A/c; w=dphi*sin(A+q*Pi)*(M/c); if (abs(ex)<0.1, w0=w; w=dphi*sin(A*ex)*(M/c); if (wsign[j]==0, wsign[j]=if(sign(w0)==sign(w), 1, -1)); w=w*wsign[j]; wsign[j]=-wsign[j]); s0=s; term=eval(Str(f "(x)")); s=s+w*term; if (s0==s, active[j]=0); ))); /*if, for, while */ s*h } \\ Approximate integral(0,inf,x,f(x)*sin(c*x+q*Pi)) \\ by the Ooura-Mori transformation using initial step size Pi/M, \\ with repeated step-halving until convergence to tolerance tol. {q_ossinf(f,M,c,q,tol)= local(olds,s,olderr,err); olds=0; olderr=10000; while(1, s=q_ossinfM(f,M,c,q); err=abs(s-olds); if (err<=tol, break); if (err<1e-9 & err>olderr, s=olds; break); olds=s; olderr=err; M=M*2); s }