C++ implementation of Longman's method, and special functions

[muppet]

Hi all,

I'm a programming newbie teaching myself C++ mainly for interest/ because I might want a real job after my physics PhD, but I have a problem in my research some code might be useful for.

I have some functions defined by integrals of the form
$$A(q)=\int db~ b* J_{0}(b*q) (exp(i*f(b))-1)$$
where I have to consider many different functions f(b), which are typically sums of products of special functions, so an analytic approach is right out. (J is a Bessel function.) My supervisor doesn't completely trust the numerics of Mathematica so he wants an independent check on the results it gives me, one where we understand the methods used. The integration runs from 0 to infinity, and the functions f(b) sometimes encounters logarithmic singularities for small arguments, and it's of course highly oscillatory, so I'd expect that any kind of standard off-the-shelf integration routine would run into problems.

Looking around I came across
http://scicomp.stackexchange.com/qu...rt-in-highly-oscillatory-integral-computation
which suggests a couple of methods tailored to oscillatory integrals over infinite regions. It seems though that if this is a standard problem that people often have, there should be some existing code in a library somewhere for it, but I can't find any. Additionally, I can't find a library of special functions that contains all the ones I need. The best I've found so far is http://www.gnu.org/software/gsl/manual/html_node/, but this doesn't have all the hypergeometric functions I need- in particular, it doesn't have 2F3 or 1F2. Can anyone recommend a good place to look for either a library for these functions or for an implementation of Longman's method for dealing with oscillatory integrals over infinite regions? As a total novice I'd rather not go trying to reinvent the wheel, badly, when the result actually matters...

Thanks in advance.

Mod note: I fixed the integral in the text above.

 

[chiro]

Hey muppet.

There are books like the one called Numerical Recipes:

http://www.nr.com/

You're problem seems rather specialized but the above kind of resource is one that is popular when you need to look for numerical routines.

Also I don't know if you are going to find something that meets the specifications exactly for your type of problem. One question I have that you might want to consider is the payoff between time searching for an algorithm to get something with a desired accuracy and the time used to say evaluate the integral using a good enough algorithm with finely tuned parameters (even if this took say a week to run).

If you could for example take a good algorithm with known error metrics and calculate the time it would take to run this with the adequate parameters, then this could give you an indication of whether it might be useful to run the simulation (even if you did it multiple times in parallel on different machines so that you had a few results to compare) with these parameters or to go to the trouble of searching, coding up, and finally evaluting the problem.

I don't know the nature of your work and if this is a one off thing or not, but it's certainly something to consider.

 

[muppet]

Hi chiro,

Thanks for your reply. As a starting point, I've been trying to get the GNU library working. (documentation: http://linux.math.tifr.res.in/programming-doc/gsl/gsl-ref_toc.html )

I've been getting an error message: undefined reference to gsl_sf_bessel_J0.

Looking at the .h file it tells me to include, gsl_sf_bessel_J0. is declared there, but not defined there, and I can't see it defined in any of the .h files it includes there either. I'm trying to compile it using the command
g++ -Wall -L/usr/lib -lgsl -lm -lgslcblas SourceCode/gsltest.cpp -o gsltest.exe
-which should I think include all the relevant libraries? (I know it's usually /usr/local/lib/; for me that's just Python seemingly. I installed the package via ubuntu software centre if that's relevant.)
My code:
//gsltest.cpp
#include <iostream>
#include <gsl/gsl_sf_bessel.h>

double test(double x);

int main()
{
double x = 5.0;
double y = gsl_sf_bessel_J0(x);

std::cout<< y ;
return 0;
}

This should be compared to http://www.gnu.org/software/gsl/manual/gsl-ref.html#An-Example-Program
- this is C code that I've tried to adapt to the C++ I'm learning. The library is also written in C, but the documentation claims that this is accounted for in the header such that it can be called from C++ code straightforwardly; I think the relevant bit of code is
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS/* Regular Bessel Function J_0(x)
*
* exceptions: none
*/

double gsl_sf_bessel_J0(const double x);
__END_DECLS

The rest is mostly just more declarations. Can anyone help me work out how to include the actual definition of J0?

EDIT: P.S. thanks mod!

 

[chiro]

Usually in libraries and external packages, the library/package comes with a "super header" which includes not only all the definitions, but includes them in the right order (which is important for avoiding compilation errors).

I'm not familiar with your package, but does your documentation tell you about a main header that includes everything?

 

[muppet]

Not as far as I can make out. If there were, I'd expect it to appear in the example code linked above...

I've also come across http://ubuntuforums.org/showthread.php?t=711169 which to be honest leaves me slightly more confused than ever.

 

[chiro]

That thread basically tells you how to add directories for include files and for library files (i.e. the compiled files that represent the actual compiled code for these modules).

It's the same thing as if you used Visual Studio and you added directories for include files and for library files except that this tells you how to do it from the command line of your compiler software (compiler, linker, etc) and not using an IDE.

Basically you need to add a path to your directory of header files using -l and add a path to all the lib files using -L. If the compiler finds the right include files and the right library files and has all the definitions, it will compile everything A-OK and you'll get some compiled module that you pass to your linker to finally generate an executable.

 

[muppet]

You see I've used all of those commands to try and include all the relevant libraries. (See post 1.) I'm wondering if there's just a chain of include statements to get the right definition that assumes that something somewhere is in some default location that's different to where the ubuntu software centre has put it. I might try reinstalling the library by a different method, once Mathematica finishes what it's been chewing over for 19hrs and counting...

 

[chiro]

Can you do a file search? I haven't used linux in a long time but I think the command is grep if I'm not mistaken.

In Visual studio you can do a find in files that will search all subdirectories recursively. If the definition exists, it should be able to be found.

 

[muppet]

Sure; in ubuntu I can just search the file system the same way I could search C: in Windows. The trouble is that I don't know what I'm looking for!

The documentation for the special functions claims that the Bessel functions are defined in the header file `gsl_sf_bessel.h', but this doesn't seem to be true; as far as I can tell they're only declared there. If I take away the #include statement for this file it changes the error message to one about the function not being declared, so I'm pretty sure it's finding it correctly. When I look in the .h file in a text editor, I only see declarations and not a definition.

 

[chiro]

Search for the definitions that are giving you errors and that the compiler can't find.

 

[muppet]

Ah I mistook your meaning before- thanks for that idea.

It turns out that the problem is pathetically easy to solve
http://ubuntuforums.org/archive/index.php/t-1833136.html

It seems that in earlier versions of Ubuntu the order in which you wrote out all the compiler options didn't matter, but in later versions it does. Most of the documentation I'd read about how to compile things including the gsl library was more than 12 months old, so I hadn't picked up on this. The command
g++ SourceCode/gsltest.cpp -o gsltest.exe -L/usr/lib -lgsl -lm -lgslcblas
works fine!

One final question, if I may. I don't really understand well how compiling something with a library works. Including a header, I get. You can't be bothered typing out "float function1(float, float)... int function100(int)" every time you want to write a program using some or all of the functions 1-100, so you write a header file function.h once and the command #include<function.h> tells your preprocessor to stick those lines of code into your source code at that point.

The motivation for using a library seems to be similar, but how does it work? If it's the library that contains the function definitions, why do I have to bother declaring the function at the beginning- why doesn't the compiler just stick the definition there? And what's the advantage of having the definition incorporated into your library somehow, rather than just writing the definition in the header file?

If you know a good resource to try and understand this, I'd be grateful. All I can find via google is either "All you have to do to use a library is..." or "This is how you create a shared library; static libraries are easy."

Thanks again for your help.

 

[jtbell]

Keep in mind that there are two separate stages involved in generating your executable program, that are invoked in sequence by the g++ command.

1. First, the compiler translates your C++ source code, along with any #included files, into object code (machine language code). The #included header (.h) file for a library contains only enough information, in the form of function prototypes, to allow the compiler to verify that you are calling (invoking) the library's functions properly: number and types of arguments, and the function's return type. The compiled libraries that actually contain the library code don't come into play here.

It's important to note here that #include <whatever> does not actually include the source code for the functions defined in library "whatever," only the prototypes for the functions. The functions have already been compiled, and their object code is in the corresponding compiled library files. (Exception: in C++, template functions do have their source code in the header files.)

2. Then, the linker combines your compiled object code with the object code from the compiled libraries to produce your complete program. This is where the "-lxxx" switches come in; they tell the linker to look for a compiled library file named something like libxxx.o or libxxx.dylib in whatever standard location (directory) the compiler expects to find libraries. At this point the header files and indeed the rest of your source code have effectively "disappeared;" they're not relevant to this step.

 

[muppet]

Ah I think I understand now. Thanks!