Skip to content

Commit

Permalink
Merge pull request #3 from Cactus-proj/dev
Browse files Browse the repository at this point in the history
test: exp_int & bessel/Airy
inkydragon authored Dec 8, 2024

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature.
2 parents 7f7f2dd + 8d3ab11 commit bd773ce
Showing 10 changed files with 370 additions and 104 deletions.
8 changes: 8 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
@@ -17,6 +17,14 @@ Therefore, contributions regarding precision and error-related contributions are
- [Markdown Docs](doc/markdown/index.md)

## Test

```sh
# On Linux
cmake -DCMAKE_BUILD_TYPE=Coverage -S . -B build && cmake --build build
cd build/ && ctest && make coverage_html
```

## License

Using a `BSD-3-Clause` like [LICENSE](LICENSE.txt).
44 changes: 44 additions & 0 deletions doc/markdown/bessel.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
# cephes/bessel

## Airy Functions

* **airy**, [Airy functions](doubldoc.md#airy)

## Hypergeometric Functions

* **hyp2f1**, [Gauss hypergeometric function](doubldoc.md#hyp2f1)
* **hyperg**, [Confluent hypergeometric function](doubldoc.md#hyperg)

## Bessel functions

* **j0**, [Bessel function of order zero](doubldoc.md#j0)
* **j1**, [Bessel function of order one](doubldoc.md#j1)
* **jn**, [Bessel function of integer order](doubldoc.md#jn)
* **jv**, [Bessel function of noninteger order](doubldoc.md#jv)

## Modified Bessel functions

* **i0**, [Modified Bessel function of order zero](doubldoc.md#i0)
* **i0e**, [Exponentially scaled modified Bessel function of order zero](doubldoc.md#i0e)
* **i1**, [Modified Bessel function of order one](doubldoc.md#i1)
* **i1e**, [Exponentially scaled modified Bessel function of order one](doubldoc.md#i1e)
* **iv**, [Modified Bessel function of noninteger order](doubldoc.md#iv)

## Modified Bessel functions - second kind

* **y0**, [Bessel function of the second kind, order zero](doubldoc.md#y0)
* **y1**, [Bessel function of the second kind, order one](doubldoc.md#y1)
* **yn**, [Bessel function of second kind of integer order](doubldoc.md#yn)

## Modified Bessel functions - third kind

* **k0**, [Modified Bessel function, third kind, order zero](doubldoc.md#k0)
* **k0e**, [Modified Bessel function, third kind, order zero, exponentially scaled](doubldoc.md#k0e)
* **k1**, [Modified Bessel function, third kind, order one](doubldoc.md#k1)
* **k1e**, [Modified Bessel function, third kind, order one, exponentially scaled](doubldoc.md#k1e)
* **kn**, [Modified Bessel function, third kind, integer order](doubldoc.md#kn)

## Others

* **psi**, [Psi (digamma) function](doubldoc.md#psi)
* **struve**, [Struve function](doubldoc.md#struve)
8 changes: 8 additions & 0 deletions doc/markdown/exp_int.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
# cephes/Exponential integral

* **expn**, [Exponential integral En](doubldoc.md#expn)

## Sine and Cosine integral

* **sici**, [Sine and cosine integrals](doubldoc.md#sici)
* **shichi**, [Hyperbolic sine and cosine integrals](doubldoc.md#shichi)
75 changes: 0 additions & 75 deletions doc/markdown/index.md
Original file line number Diff line number Diff line change
@@ -2,63 +2,6 @@

Provides and wraps the mathematical functions from the [Cephes mathematical library](http://www.netlib.org/cephes/), developed by [Stephen L. Moshier](http://www.moshier.net). This C library provides a <b>lot</b> of mathematical functions. It is used, among many other places, [at the heart of SciPy](https://github.com/scipy/scipy/tree/master/scipy/special/cephes).

## Example

### Simple call on a number

The wrapped functions can be called from Lua with the same synopsis as their C coutnerpart, and will then return a number, for example:

```lua
require 'cephes'
x = cephes.igam(2, 3) -- returns a number
```

### Calling on tensors

Our wrappers for cephes functions are vectorized, meaning they can

* take tensors as arguments, evaluating the function for each element of the arguments, and return the result into a vector.
* mix tensors and numbers as arguments, numbers are automatically expanded
* **shape does not matter**, only the number of elements.

Like most torch functions, they also accept an optional Tensor as first argument to store the result into.

```lua
require 'cephes'
-- Call over a whole tensor of parameters
result = cephes.ndtr(torch.randn(10)) -- returns a new tensor
-- of 10 elements

-- Several tensor arguments, pairing them map-like
-- Below returns a vector of 100 elements
x = torch.rand(100)
y = torch.rand(100)
result = cephes.igam(x, y)

-- Mix number and tensors: numbers are automatically expanded
-- Below returns a vector of 100 elements
result = cephes.igam(4, y)

-- Can also use matrices: only the number of elements matters
-- Below with a 3D array and a vector, return a vector of 100 elts
z = torch.rand(2,5,10)
result = cephes.igam(z, y)

-- And of course you can store the result into an
-- existing tensor of the right number of elements
-- Below stores into an existing 3D tensors of 100 elements
result = torch.Tensor(2,5,10)
cephes.igam(result, x, y)
```

## Installation

From a terminal:

```bash
torch-rocks install cephes
```

## Error Handling

By default, Torch-Cephes **does not signal any error** (domain, singularity, overflow, underflow, precision). It is as non-intrusive as possible and tries to return a value which is hopefully usable: it might be NaN, it might be inf.
@@ -320,21 +263,3 @@ Checks if a number is finite.
> * `im` : any number, to initalize the imaginary part
>
>**Returns:** a pointer to a new Cephes FFI complex number with real part `r` and imaginary part `im`.
## Unit Tests

Last but not least, the unit tests are in the folder
[`luasrc/tests`](https://github.com/jucor/torch-cephes/tree/master/luasrc/tests). You can run them from your local clone of the repostiory with:

```bash
git clone https://www.github.com/jucor/torch-cephes
find torch-cephes/luasrc/tests -name "test*lua" -exec torch {} \;
```

Those tests will soone be automatically installed with the package, once I sort out a bit of CMake resistance.

## Direct access to FFI

### cephes.ffi.*

Functions directly accessible at the top of the `cephes` table are Lua wrappers to the actual C functions from Cephes, with extra error checking. If, for any reason, you want to get rid of this error checking and of a possible overhead, the FFI-wrapper functions can be called directly via `cephes.ffi.myfunction()` instead of `cephes.myfunction()`.
37 changes: 20 additions & 17 deletions include/cephes/double.h → include/cephes/bessel.h
Original file line number Diff line number Diff line change
@@ -1,48 +1,51 @@
#ifndef CEPHES_DOUBLE_H
/** Cephes double precision special functions suite */
#define CEPHES_DOUBLE_H
#ifndef CEPHES_BESSEL_H
/** Cephes double precision special functions suite
*
* cephes/bessel
*/
#define CEPHES_BESSEL_H

#if defined(__cplusplus)
extern "C" {
#endif

/** cephes/bessel */

/** Airy Functions */
int airy(double x, double *ai, double *aip, double *bi, double *bip);

/** Hypergeometric Functions */
double hyp2f1(double a, double b, double c, double x);
double hyperg(double a, double b, double x);
double hyp2f0(double a, double b, double x, int type, double *err);

/** Bessel functions */
double j0(double x);
double j1(double x);
double jn(int n, double x);
double jv(double n, double x);

/** Modified Bessel functions */
double i0(double x);
double i0e(double x);
double i1(double x);
double i1e(double x);
double iv(double v, double x);

double j0(double x);
/** Modified Bessel functions - second kind */
double y0(double x);
double j1(double x);
double y1(double x);
double jn(int n, double x);
double jv(double n, double x);
double yn(int n, double x);

/** Modified Bessel functions - third kind */
double k0(double x);
double k0e(double x);
double k1(double x);
double k1e(double x);
double kn(int nn, double x);

/** Others */
double psi(double x);

double onef2(double a, double b, double c, double x, double *err);
double threef0(double a, double b, double c, double x, double *err);
double struve(double v, double x);
double yv(double v, double x);

double yn(int n, double x);

#if defined(__cplusplus)
}
#endif
#endif // CEPHES_DOUBLE_H
#endif // CEPHES_BESSEL_H
25 changes: 25 additions & 0 deletions include/cephes/exp_int.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
#ifndef CEPHES_EXP_INT_H
/** Cephes double precision special functions suite
*
* cephes/Exponential integral
*/
#define CEPHES_EXP_INT_H

#if defined(__cplusplus)
extern "C" {
#endif

/** Exponential integral */
/* misc/expn.c */
double expn(int n, double x);

/** Sine and Cosine integral */
/* misc/sici.c */
int sici(double x, double *si, double *ci);
/** misc/shichi.c */
int shichi(double x, double *si, double *ci);

#if defined(__cplusplus)
}
#endif
#endif // CEPHES_EXP_INT_H
4 changes: 4 additions & 0 deletions tests/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -12,6 +12,9 @@ FetchContent_MakeAvailable(googletest)

include(GoogleTest)

# Add include for tests
include_directories(${CMAKE_SOURCE_DIR}/tests/include)

# Add Tests
function(add_gtest test_name)
add_executable(${test_name} ${test_name}.cpp)
@@ -20,3 +23,4 @@ function(add_gtest test_name)
endfunction()

add_gtest(airy)
add_gtest(exp_int)
84 changes: 72 additions & 12 deletions tests/airy.cpp
Original file line number Diff line number Diff line change
@@ -1,18 +1,78 @@
#include <limits>
#include <complex>
#include <gtest/gtest.h>
#include <cephes/double.h>
#include <cephes/bessel.h>
#include <xtest.hpp>

TEST(Airy, BasicAssertions) {
const double nan64 = std::numeric_limits<double>::quiet_NaN();
int ret;
double x, ai, aip, bi, bip;
const double nan64 = std::numeric_limits<double>::quiet_NaN();
int ret;
double x, ai, aip, bi, bip;

x = 26.0;
ret = airy(x, &ai, &aip, &bi, &bip);
EXPECT_EQ(ret, -1);
EXPECT_EQ(ai, 0.0);
EXPECT_EQ(aip, 0.0);
EXPECT_TRUE(bi > 1e308);
EXPECT_TRUE(bip > 1e308);
x = 26.0;
ret = airy(x, &ai, &aip, &bi, &bip);
EXPECT_EQ(ret, -1);
EXPECT_EQ(ai, 0.0);
EXPECT_EQ(aip, 0.0);
EXPECT_TRUE(bi > 1e308);
EXPECT_TRUE(bip > 1e308);
}

TEST(Airy, Branches) {
int ret;
double x, ai, aip, bi, bip, ref;
double ai_ref, aip_ref, bi_ref, bip_ref;

// x < -2.09
x = -3.0;
ret = airy(x, &ai, &aip, &bi, &bip);
EXPECT_EQ(ret, 0);
/*
x = -3.0;
{NumberForm[AiryAi[x], 16], NumberForm[AiryAiPrime[x], 16],
NumberForm[AiryBi[x], 16], NumberForm[AiryBiPrime[x], 16]}
*/
ai_ref = -0.378814293677658;
aip_ref = 0.3145837692165989;
bi_ref = -0.1982896263749266;
bip_ref = -0.6756112226852585;
XTEST_ISAPPROX_F64(ai);
XTEST_ISAPPROX_F64(aip);
XTEST_ISAPPROX_F64(bi);
XTEST_ISAPPROX_F64(bip);

// x >= 2.09 (cbrt(9))
x = 5.0;
ret = airy(x, &ai, &aip, &bi, &bip);
EXPECT_EQ(ret, 0);
/*
x = 5.0;
{NumberForm[AiryAi[x], 16], NumberForm[AiryAiPrime[x], 16],
NumberForm[AiryBi[x], 16], NumberForm[AiryBiPrime[x], 16]}
*/
ai_ref = 0.0001083444281360744;
aip_ref = -0.0002474138908684624;
bi_ref = 657.7920441711713;
bip_ref = 1435.81908021797;
XTEST_ISAPPROX_F64(ai);
XTEST_ISAPPROX_F64(aip);
XTEST_ISAPPROX_F64(bi);
XTEST_ISAPPROX_F64(bip);

// x > 8.3203353
x = 8.5;
ret = airy(x, &ai, &aip, &bi, &bip);
EXPECT_EQ(ret, 0);
/*
x = 8.5;
{NumberForm[AiryAi[x], 16], NumberForm[AiryAiPrime[x], 16],
NumberForm[AiryBi[x], 16], NumberForm[AiryBiPrime[x], 16]}
*/
ai_ref = 1.099700975519552e-8;
aip_ref = -3.237725440447604e-8;
bi_ref = 4.965319541471301e6;
bip_ref = 1.432630103066198e7;
XTEST_ISAPPROX_F64(ai);
XTEST_ISAPPROX_F64(aip);
XTEST_ISAPPROX_F64(bi);
XTEST_ISAPPROX_F64(bip);
}
Loading

0 comments on commit bd773ce

Please sign in to comment.