SERVING THE QUANTITATIVE FINANCE COMMUNITY

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

quadruple precision complex erfc code (qerfc) in fortran

January 17th, 2018, 8:50 am

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
 
User avatar
outrun
Posts: 4573
Joined: April 29th, 2016, 1:40 pm

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

January 17th, 2018, 9:36 am

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?
 
User avatar
ppauper
Topic Author
Posts: 70239
Joined: November 15th, 2001, 1:29 pm

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

January 17th, 2018, 10:11 am

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
 
User avatar
Cuchulainn
Posts: 59713
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

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

January 17th, 2018, 10:19 am

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.
 
User avatar
Cuchulainn
Posts: 59713
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

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

January 17th, 2018, 10:24 am

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.
 
User avatar
Cuchulainn
Posts: 59713
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

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

January 17th, 2018, 10:25 am

Maybe simple integration with Gaussian Quadrature is easier to implement?
I think this is a non-starter for at least 2 reasons.
 
User avatar
outrun
Posts: 4573
Joined: April 29th, 2016, 1:40 pm

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

January 17th, 2018, 11:16 am

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.
 
User avatar
outrun
Posts: 4573
Joined: April 29th, 2016, 1:40 pm

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

January 17th, 2018, 11:30 am

Ah, .. *complex*, missed that, that makes it even more tricky.
 
User avatar
Cuchulainn
Posts: 59713
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

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

January 17th, 2018, 1:15 pm

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? 
 
User avatar
outrun
Posts: 4573
Joined: April 29th, 2016, 1:40 pm

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

January 17th, 2018, 1:16 pm

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?
 
User avatar
ppauper
Topic Author
Posts: 70239
Joined: November 15th, 2001, 1:29 pm

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

January 17th, 2018, 1:29 pm

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
 
User avatar
outrun
Posts: 4573
Joined: April 29th, 2016, 1:40 pm

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

January 17th, 2018, 1:53 pm

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.
 
User avatar
Cuchulainn
Posts: 59713
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

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

January 18th, 2018, 6:35 am

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/

add additional terms to the series to get the accuracy needed
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.
 
User avatar
ppauper
Topic Author
Posts: 70239
Joined: November 15th, 2001, 1:29 pm

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

January 18th, 2018, 8:39 am

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
 
User avatar
Cuchulainn
Posts: 59713
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

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

January 18th, 2018, 4:51 pm

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

BTW how did they come up with those numbers?
ABOUT WILMOTT

PW by JB

Wilmott.com has been "Serving the Quantitative Finance Community" since 2001. Continued...


Twitter LinkedIn Instagram

JOBS BOARD

JOBS BOARD

Looking for a quant job, risk, algo trading,...? Browse jobs here...


GZIP: On