it's for a pricing problem, and I had been doing it it maple (with Digits=31) but there's a linear algebra problem involved and my laptop can't cope with the memory needed for accurate solutions (accurate = agrees with solutions found by other methods) so I had to port it to fortran

the double precision erfc code I've been using (off the web, and no, there wasn't a quadruple precision version there) is

which works great in double precision but isn't accurate enough for my purposes.function erfc(x)

pv= 9.03777677D0

ph= 4.76888839D0

p0= 3.92213346D-1

p1= 1.49181256D-1

p2= 2.15823157D-2

p3= 1.18760964D-3

p4= 2.48565745D-5

p5= 1.97879468D-7

q0= 1.20830487D-1

q1= 1.08747438D0

q2= 3.02076217D0

q3= 5.92069385D0

q4= 9.78726942D0

q5= 1.46204889D1

r0= 2.21293361D-1

r1= 2.72957057D-1

r2= 6.40298026D-2

r3= 5.71296931D-3

r4= 1.93880223D-4

r5= 2.50263215D-6

r6= 1.22871857D-8

s1= 4.83321947D-1

s2= 1.93328779D0

s3= 4.34989752D0

s4= 7.73315115D0

s5= 1.20830487D1

s6= 1.73995901D1

if(abs(real(x))+abs(aimag(x)).lt.ph) then

z=exp(pv*x)

if(real(z).ge.0) then

y=exp(-y)*x*(p5/(y+q5)+p4/(y+q4)+p3/(y+q3)+p2/(y+q2)+p1/(y+q1)+p0/(y+q0))+2/(1+z)

else

y=exp(-y)*x*(r6/(y+s6)+r5/(y+s5)+r4/(y+s4)+r3/(y+s3)+r2/(y+s2)+r1/(y+s1)+r0/y)+2/(1-z)

endif

else

y=exp(-y)*x*(p5/(y+q5)+p4/(y+q4)+p3/(y+q3)+p2/(y+q2)+p1/(y+q1)+p0/(y+q0))

if(real(x).lt.0) y=y+2

endif

werfc=y

If I knew where the constants here came from, I could probably calculate more accurate versions and/or add additional terms to the series to get the accuracy needed