Serving the Quantitative Finance Community

• 1
• 2

ppauper
Topic Author
Posts: 70239
Joined: November 15th, 2001, 1:29 pm

### quadruple precision complex erfc code (qerfc) in fortran

does anyone have quadruple precision fortran code for complex erfc (complementary error function) they can let me have

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
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
which works great in double precision but isn't accurate enough for my purposes.
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

outrun
Posts: 4573
Joined: April 29th, 2016, 1:40 pm

### Re: quadruple precision complex erfc code (qerfc) in fortran

Tricky!

Abramowitz Stegun doesn't list it.

In the c++ boost math library they have a version that has a max error of 3E-35 for quadrupal precision, the source code is here http://www.boost.org/doc/libs/1_65_0/bo ... ns/erf.hpp

it's a bit involved to transform in to Fortran I guess?

Maybe simple integration with Gaussian Quadrature is easier to implement?

ppauper
Topic Author
Posts: 70239
Joined: November 15th, 2001, 1:29 pm

### Re: quadruple precision complex erfc code (qerfc) in fortran

porting that code to fortran is beyond my capabilities

And quadrature doesn't seem to be an option as we'd be here to Christmas while it runs.
I've got a (2N)x(2N) matrix, and each of the elements has a term of the form $e^{x^2}\textrm{erfc}(x)$ and it looks like I need N to be at least 200 which would mean 160,000 calls

the algorithm I've already got would be perfect if only it were quadruple precision, because it involves $e^{-x^2}$ so I can throw the $e^{x^2}$ factor into the algorithm and cut the computational time considerably

Cuchulainn
Posts: 64298
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

### Re: quadruple precision complex erfc code (qerfc) in fortran

First impressions.
This is a problem similar to a discussion with Alan regarding Kummere's Hypergeometric fuinction in  the most general case in which all parameters are complex.

viewtopic.php?f=10&t=100245&start=240

The current case is a bit easier.

A possible route is erfc(z) -> incomplete gamma function -> Kummer's function

https://ac.els-cdn.com/0377042785900342 ... 22da0164f2

I used to do Fortran a lot and C code is similar. Boost won't be of much hep IMO because it does not support complex argumetns. We discussed this on Wilmott already.  So, you need to us the above or maybe a modified form of Cornelius Lanczos' approximation.

With a good algorithm Quad precision is not an issue.
Last edited by Cuchulainn on January 17th, 2018, 10:42 am, edited 1 time in total.
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl

Cuchulainn
Posts: 64298
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

### Re: quadruple precision complex erfc code (qerfc) in fortran

the algorithm I've already got would be perfect if only it were quadruple precision,

Can you post the Fortran code/algo; It should be easy to port to C.

BTW would gfortran not solve the problem?(?)

https://people.sc.fsu.edu/~jburkardt/f7 ... dmath.html

You might have to split into real and complex parts which is the approach I used for $M(a,b.z)$.
Last edited by Cuchulainn on January 17th, 2018, 10:35 am, edited 4 times in total.
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl

Cuchulainn
Posts: 64298
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

### Re: quadruple precision complex erfc code (qerfc) in fortran

Maybe simple integration with Gaussian Quadrature is easier to implement?
I think this is a non-starter for at least 2 reasons.
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl

outrun
Posts: 4573
Joined: April 29th, 2016, 1:40 pm

### Re: quadruple precision complex erfc code (qerfc) in fortran

the algorithm I've already got would be perfect if only it were quadruple precision,

Can you post the Fortran code/algo; It should be easy to port to C.

BTW would gfortran not solve the problem?(?)

https://people.sc.fsu.edu/~jburkardt/f7 ... dmath.html

You might have to split into real and complex parts which is the approach I used for $M(a,b.z)$.
Why port the low precision Fortran code to C? He needs high precision Fortran code.

outrun
Posts: 4573
Joined: April 29th, 2016, 1:40 pm

### Re: quadruple precision complex erfc code (qerfc) in fortran

Ah, .. *complex*, missed that, that makes it even more tricky.

Cuchulainn
Posts: 64298
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

### Re: quadruple precision complex erfc code (qerfc) in fortran

the algorithm I've already got would be perfect if only it were quadruple precision,

Can you post the Fortran code/algo; It should be easy to port to C.

BTW would gfortran not solve the problem?(?)

https://people.sc.fsu.edu/~jburkardt/f7 ... dmath.html

You might have to split into real and complex parts which is the approach I used for $M(a,b.z)$.
Why port the low precision Fortran code to C? He needs high precision Fortran code.
Are you sure?
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl

outrun
Posts: 4573
Joined: April 29th, 2016, 1:40 pm

### Re: quadruple precision complex erfc code (qerfc) in fortran

the algorithm I've already got would be perfect if only it were quadruple precision,

Can you post the Fortran code/algo; It should be easy to port to C.

BTW would gfortran not solve the problem?(?)

https://people.sc.fsu.edu/~jburkardt/f7 ... dmath.html

You might have to split into real and complex parts which is the approach I used for $M(a,b.z)$.
Why port the low precision Fortran code to C? He needs high precision Fortran code.
Are you sure?
What does the first sentence say?

ppauper
Topic Author
Posts: 70239
Joined: November 15th, 2001, 1:29 pm

### Re: quadruple precision complex erfc code (qerfc) in fortran

I don't speak C++ and if the suggestion is that I port the code to C++, that's a recipe for disaster as I won't understand the code

as per the title of the thread, I need a quadruple precision complex erfc code (qerfc) in fortran
I have a working double precision code but it needs to be more accurate, hence I need to replace the double precision erfc code with a quadruple precision code

The test case I'm running has an answer in the literature of 1.50448.
With the double precision error function and a 300x300 matrix, I'm getting 1.50377

if anyone recognizes the erfc algorithm I posted earlier and can direct me to the underlying paper so I can find and/or compute the constants for a quadruple precision version, that would work

outrun
Posts: 4573
Joined: April 29th, 2016, 1:40 pm

### Re: quadruple precision complex erfc code (qerfc) in fortran

I think you can find it here (the new constants) there seems to be quadruple.. "complex quad erfc(complex quad x) "

https://github.com/philscher/gkc/blob/m ... alMath.cpp

Found it by googling constants from your code.

Cuchulainn
Posts: 64298
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

### Re: quadruple precision complex erfc code (qerfc) in fortran

Plan D: find the algorithm for this and implement it in FORTRAN yourself, or find code.  Or an MSc student to help.

If I were to do it I would first have a good look at A&S 7..2.12 where erfc(z) in written in terms of the gamma function and Kummer's hypergeometric function.The latter function can be done in FORTRAN.

http://www.ece.mtu.edu/faculty/wfp/arti ... h_soft.pdf

There is also a question of the rationale for quad precision in the first place, unless to resolve issues with REAL*16 etc. A better algo might avoid having to quad? It can be a sledgehammer.

The test case I'm running has an answer in the literature of 1.50448.
With the double precision error function and a 300x300 matrix, I'm getting 1.50377

What if you do 600X600 (1200/1200, whatever) ? Does convergence improve at all?  A hypothesis is the numerical algorithms is losing digits or something, caused by series summation, catastrophic cancellation, that kind of thing.

That's one way to approach it.

And CONHYP source code (in FORTRAN) seems to be here. It claims accuracy to 13 significant digits.

http://computers.interactiva.org/Progra ... Functions/

Ah, series, always a pain in the derry air. How do you sum it?
IMO it might be an idea to analyse the current algo in use with a view to providing a plan B alternative.

// I can't imagine that it hasn't been done at e.g. CERN, NIST etc.
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl

ppauper
Topic Author
Posts: 70239
Joined: November 15th, 2001, 1:29 pm

### Re: quadruple precision complex erfc code (qerfc) in fortran

I think you can find it here (the new constants) there seems to be quadruple.. "complex quad erfc(complex quad x) "

https://github.com/philscher/gkc/blob/m ... alMath.cpp

Found it by googling constants from your code.
that seems to work perfectly, thanks

Cuchulainn
Posts: 64298
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

### Re: quadruple precision complex erfc code (qerfc) in fortran

Was erfc() good enough or was qerfc() necessary?

BTW how did they come up with those numbers?
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl