stuff
This commit is contained in:
William Ball
parent
a66d06c648
commit
55dbfd8335
6 changed files with 791 additions and 88 deletions
140
fractal_dimension/bai_finch_rust/Cargo.lock
generated
140
fractal_dimension/bai_finch_rust/Cargo.lock
generated
|
|
@ -11,6 +11,15 @@ dependencies = [
|
|||
"num-traits",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "approx"
|
||||
version = "0.5.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "072df7202e63b127ab55acfe16ce97013d5b97bf160489336d3f1840fd78e99e"
|
||||
dependencies = [
|
||||
"num-traits",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "autocfg"
|
||||
version = "1.0.1"
|
||||
|
|
@ -23,8 +32,8 @@ version = "0.11.0"
|
|||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "e5d568d86ec0ca7d24caa67afff9605bc5bc3cae2642d6670ea5924b839ccabc"
|
||||
dependencies = [
|
||||
"nalgebra",
|
||||
"num-complex",
|
||||
"nalgebra 0.24.1",
|
||||
"num-complex 0.3.1",
|
||||
"num-traits",
|
||||
"phf",
|
||||
"phf_codegen",
|
||||
|
|
@ -35,10 +44,10 @@ name = "bai_finch"
|
|||
version = "0.1.0"
|
||||
dependencies = [
|
||||
"bacon-sci",
|
||||
"nalgebra",
|
||||
"nalgebra 0.27.1",
|
||||
"num",
|
||||
"num-complex",
|
||||
"simba",
|
||||
"num-complex 0.3.1",
|
||||
"simba 0.3.1",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
|
|
@ -89,24 +98,60 @@ dependencies = [
|
|||
"rawpointer",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "matrixmultiply"
|
||||
version = "0.3.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "5a8a15b776d9dfaecd44b03c5828c2199cddff5247215858aac14624f8d6b741"
|
||||
dependencies = [
|
||||
"rawpointer",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "nalgebra"
|
||||
version = "0.24.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "1a9002895a0de45e3cde58b36d0cf2f83249e7dba4a43ee64dafbf01bfd464ff"
|
||||
dependencies = [
|
||||
"approx",
|
||||
"approx 0.4.0",
|
||||
"generic-array",
|
||||
"matrixmultiply",
|
||||
"num-complex",
|
||||
"num-rational",
|
||||
"matrixmultiply 0.2.4",
|
||||
"num-complex 0.3.1",
|
||||
"num-rational 0.3.2",
|
||||
"num-traits",
|
||||
"rand",
|
||||
"rand_distr",
|
||||
"simba",
|
||||
"simba 0.3.1",
|
||||
"typenum",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "nalgebra"
|
||||
version = "0.27.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "462fffe4002f4f2e1f6a9dcf12cc1a6fc0e15989014efc02a941d3e0f5dc2120"
|
||||
dependencies = [
|
||||
"approx 0.5.0",
|
||||
"matrixmultiply 0.3.1",
|
||||
"nalgebra-macros",
|
||||
"num-complex 0.4.0",
|
||||
"num-rational 0.4.0",
|
||||
"num-traits",
|
||||
"simba 0.5.1",
|
||||
"typenum",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "nalgebra-macros"
|
||||
version = "0.1.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "01fcc0b8149b4632adc89ac3b7b31a12fb6099a0317a4eb2ebff574ef7de7218"
|
||||
dependencies = [
|
||||
"proc-macro2",
|
||||
"quote",
|
||||
"syn",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "num"
|
||||
version = "0.3.1"
|
||||
|
|
@ -114,10 +159,10 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
|
|||
checksum = "8b7a8e9be5e039e2ff869df49155f1c06bd01ade2117ec783e56ab0932b67a8f"
|
||||
dependencies = [
|
||||
"num-bigint",
|
||||
"num-complex",
|
||||
"num-complex 0.3.1",
|
||||
"num-integer",
|
||||
"num-iter",
|
||||
"num-rational",
|
||||
"num-rational 0.3.2",
|
||||
"num-traits",
|
||||
]
|
||||
|
||||
|
|
@ -141,6 +186,15 @@ dependencies = [
|
|||
"num-traits",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "num-complex"
|
||||
version = "0.4.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "26873667bbbb7c5182d4a37c1add32cdf09f841af72da53318fdb81543c15085"
|
||||
dependencies = [
|
||||
"num-traits",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "num-integer"
|
||||
version = "0.1.44"
|
||||
|
|
@ -174,6 +228,17 @@ dependencies = [
|
|||
"num-traits",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "num-rational"
|
||||
version = "0.4.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "d41702bd167c2df5520b384281bc111a4b5efcf7fbc4c9c222c815b07e0a6a6a"
|
||||
dependencies = [
|
||||
"autocfg",
|
||||
"num-integer",
|
||||
"num-traits",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "num-traits"
|
||||
version = "0.2.14"
|
||||
|
|
@ -234,6 +299,24 @@ version = "0.2.10"
|
|||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "ac74c624d6b2d21f425f752262f42188365d7b8ff1aff74c82e45136510a4857"
|
||||
|
||||
[[package]]
|
||||
name = "proc-macro2"
|
||||
version = "1.0.27"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "f0d8caf72986c1a598726adc988bb5984792ef84f5ee5aa50209145ee8077038"
|
||||
dependencies = [
|
||||
"unicode-xid",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "quote"
|
||||
version = "1.0.9"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "c3d0b9745dc2debf507c8422de05d7226cc1f0644216dfdfead988f9b1ab32a7"
|
||||
dependencies = [
|
||||
"proc-macro2",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rand"
|
||||
version = "0.7.3"
|
||||
|
|
@ -307,8 +390,20 @@ version = "0.3.1"
|
|||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "17bfe642b1728a6e89137ad428ef5d4738eca4efaba9590f9e110b8944028621"
|
||||
dependencies = [
|
||||
"approx",
|
||||
"num-complex",
|
||||
"approx 0.4.0",
|
||||
"num-complex 0.3.1",
|
||||
"num-traits",
|
||||
"paste",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "simba"
|
||||
version = "0.5.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "8e82063457853d00243beda9952e910b82593e4b07ae9f721b9278a99a0d3d5c"
|
||||
dependencies = [
|
||||
"approx 0.5.0",
|
||||
"num-complex 0.4.0",
|
||||
"num-traits",
|
||||
"paste",
|
||||
]
|
||||
|
|
@ -319,12 +414,29 @@ version = "0.3.5"
|
|||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "cbce6d4507c7e4a3962091436e56e95290cb71fa302d0d270e32130b75fbff27"
|
||||
|
||||
[[package]]
|
||||
name = "syn"
|
||||
version = "1.0.73"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "f71489ff30030d2ae598524f61326b902466f72a0fb1a8564c001cc63425bcc7"
|
||||
dependencies = [
|
||||
"proc-macro2",
|
||||
"quote",
|
||||
"unicode-xid",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "typenum"
|
||||
version = "1.13.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "879f6906492a7cd215bfa4cf595b600146ccfac0c79bcbd1f3000162af5e8b06"
|
||||
|
||||
[[package]]
|
||||
name = "unicode-xid"
|
||||
version = "0.2.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "8ccb82d61f80a663efe1f787a51b16b5a51e3314d6ac365b08639f52387b33f3"
|
||||
|
||||
[[package]]
|
||||
name = "version_check"
|
||||
version = "0.9.3"
|
||||
|
|
|
|||
|
|
@ -8,7 +8,8 @@ edition = "2018"
|
|||
|
||||
[dependencies]
|
||||
num = "0.3.1"
|
||||
nalgebra = "0.24.1"
|
||||
# nalgebra = "0.24.1"
|
||||
nalgebra = "0.27.1"
|
||||
bacon-sci = "*"
|
||||
simba = "0.3.1"
|
||||
num-complex = "0.3.1"
|
||||
|
|
|
|||
|
|
@ -1,5 +1,4 @@
|
|||
use num_complex::Complex64;
|
||||
use simba::scalar::ComplexField;
|
||||
use std::f64::consts::PI;
|
||||
|
||||
pub type CFunction<'a> = &'a dyn Fn(Complex64) -> Complex64;
|
||||
|
|
@ -7,7 +6,7 @@ pub type Path<'a> = &'a dyn Fn(f64) -> Complex64;
|
|||
pub type Rule<'a> = &'a dyn Fn(CFunction, Complex64, Complex64) -> Complex64;
|
||||
|
||||
pub fn unit_circle(t: f64) -> Complex64 {
|
||||
num::Complex::new((2.0 * PI * t).cos(), (2.0 * PI * t).sin())
|
||||
Complex64::new((2.0 * PI * t).cos(), (2.0 * PI * t).sin())
|
||||
}
|
||||
|
||||
pub fn trapezoid(f: CFunction, a: Complex64, b: Complex64) -> Complex64 {
|
||||
|
|
@ -32,15 +31,13 @@ fn test(t: f64) -> Complex64 {
|
|||
// }
|
||||
|
||||
pub fn diff(fun: CFunction, a: Complex64, n: i32) -> Result<Complex64, String> {
|
||||
let f = |z: Complex64| -> Complex64 {
|
||||
fun(z) / (z - a).powi(n + 1)
|
||||
};
|
||||
let f = |z: Complex64| -> Complex64 { fun(z) / (z - a).powi(n + 1) };
|
||||
let integral = integrate(&f, &unit_circle, &trapezoid);
|
||||
Ok(integral * Complex64::new(factorial(n) as f64, 0.0) / Complex64::new(0.0, 2.0 * PI))
|
||||
}
|
||||
|
||||
pub fn integrate(f: CFunction, path: Path, rule: Rule) -> Complex64 {
|
||||
let increment = 0.00001;
|
||||
let increment = 0.05;
|
||||
let mut integral = Complex64::new(0.0, 0.0);
|
||||
|
||||
let mut i = increment;
|
||||
|
|
|
|||
171
fractal_dimension/bai_finch_rust/src/integation.rs
Normal file
171
fractal_dimension/bai_finch_rust/src/integation.rs
Normal file
|
|
@ -0,0 +1,171 @@
|
|||
//! The double exponential algorithm is naturally adaptive, it stops calling the integrand when the error is reduced to below the desired threshold.
|
||||
//! It also does not allocate. No box, no vec, etc.
|
||||
//! It has a hard coded maximum of approximately 350 function evaluations. This guarantees that the algorithm will return.
|
||||
//! The error in the algorithm decreases exponentially in the number of function evaluations, specifically O(exp(-cN/log(N))). So if 350 function evaluations is not giving the desired accuracy than the programmer probably needs to give some guidance by splitting up the range at singularities or [other preparation techniques](http://www.johndcook.com/blog/2012/02/21/care-and-treatment-of-singularities/).
|
||||
//!
|
||||
//! This is a port of the [Fast Numerical Integration](https://www.codeproject.com/kb/recipes/fastnumericalintegration.aspx) from c++ to rust. The original code is by John D. Cook, and is licensed under the [BSD](https://opensource.org/licenses/bsd-license.php).
|
||||
|
||||
use num_complex::Complex64;
|
||||
|
||||
/// Integrate an analytic function over a finite interval.
|
||||
/// f is the function to be integrated.
|
||||
/// a is left limit of integration.
|
||||
/// b is right limit of integration
|
||||
/// target_absolute_error is the desired bound on error
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// use quadrature::integrate;
|
||||
/// fn integrand(x: f64) -> f64 {
|
||||
/// (-x / 5.0).exp() * x.powf(-1.0 / 3.0)
|
||||
/// }
|
||||
/// let o = integrate(integrand , 0.0, 10.0, 1e-6);
|
||||
/// assert!((o.integral - 3.6798142583691758).abs() <= 1e-6);
|
||||
/// ```
|
||||
pub fn integrate<F>(f: F, a: f64, b: f64, target_absolute_error: f64) -> Complex64
|
||||
where F: Fn(f64) -> Complex64
|
||||
{
|
||||
// Apply the linear change of variables x = ct + d
|
||||
// $$\int_a^b f(x) dx = c \int_{-1}^1 f( ct + d ) dt$$
|
||||
// c = (b-a)/2, d = (a+b)/2
|
||||
let c = 0.5 * (b - a);
|
||||
let d = 0.5 * (a + b);
|
||||
integrate_core(|x| {
|
||||
let out = f(c * x + d);
|
||||
if out.is_finite() { out } else { 0.0 }
|
||||
},
|
||||
0.25 * target_absolute_error / c)
|
||||
.scale(c)
|
||||
}
|
||||
|
||||
/// Integrate f(x) from [-1.0, 1.0]
|
||||
fn integrate_core<F>(f: F, target_absolute_error: f64) -> Complex64
|
||||
where F: Fn(f64) -> Complex64
|
||||
{
|
||||
let mut error_estimate = ::std::f64::MAX;
|
||||
let mut num_function_evaluations = 1;
|
||||
let mut current_delta = ::std::f64::MAX;
|
||||
|
||||
let mut integral = 2.0 * ::std::f64::consts::FRAC_PI_2 * f(0.0);
|
||||
|
||||
for &weight in &WEIGHTS {
|
||||
let new_contribution = weight.iter()
|
||||
.map(|&(w, x)| w * (f(x) + f(-x)))
|
||||
.fold(0.0, |sum, x| sum + x);
|
||||
num_function_evaluations += 2 * weight.len();
|
||||
|
||||
// difference in consecutive integral estimates
|
||||
let previous_delta_ln = current_delta.ln();
|
||||
current_delta = (0.5 * integral - new_contribution).abs();
|
||||
integral = 0.5 * integral + new_contribution;
|
||||
|
||||
// Once convergence kicks in, error is approximately squared at each step.
|
||||
// Determine whether we're in the convergent region by looking at the trend in the error.
|
||||
if num_function_evaluations <= 13 {
|
||||
// level <= 1
|
||||
continue; // previousDelta meaningless, so cannot check convergence.
|
||||
}
|
||||
|
||||
// Exact comparison with zero is harmless here. Could possibly be replaced with
|
||||
// a small positive upper limit on the size of currentDelta, but determining
|
||||
// that upper limit would be difficult. At worse, the loop is executed more
|
||||
// times than necessary. But no infinite loop can result since there is
|
||||
// an upper bound on the loop variable.
|
||||
if current_delta == 0.0 {
|
||||
error_estimate = 0.0;
|
||||
break;
|
||||
}
|
||||
// previousDelta != 0 or would have been kicked out previously
|
||||
let r = current_delta.ln() / previous_delta_ln;
|
||||
|
||||
if r > 1.9 && r < 2.1 {
|
||||
// If convergence theory applied perfectly, r would be 2 in the convergence region.
|
||||
// r close to 2 is good enough. We expect the difference between this integral estimate
|
||||
// and the next one to be roughly delta^2.
|
||||
error_estimate = current_delta * current_delta;
|
||||
} else {
|
||||
// Not in the convergence region. Assume only that error is decreasing.
|
||||
error_estimate = current_delta;
|
||||
}
|
||||
|
||||
if error_estimate < target_absolute_error {
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
integral
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn trivial_function_works() {
|
||||
let o = integrate(|_| 0.5, -1.0, 1.0, 1e-14);
|
||||
assert!(o.error_estimate <= 1e-14,
|
||||
"error_estimate larger then asked. estimate: {:#?}, asked: {:#?}",
|
||||
o.error_estimate,
|
||||
1e-14);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn demo_function1_works() {
|
||||
let o = integrate(|x| (-x / 5.0).exp() * x.powf(-1.0 / 3.0), 0.0, 10.0, 1e-6);
|
||||
assert!(o.error_estimate <= 1e-6,
|
||||
"error_estimate larger then asked. estimate: {:#?}, asked: {:#?}",
|
||||
o.error_estimate,
|
||||
1e-6);
|
||||
assert!((o.integral - 3.6798142583691758).abs() <= 1e-6,
|
||||
"error larger then error_estimate");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn demo_function2_works() {
|
||||
let o = integrate(|x| (1.0 - x).powf(5.0) * x.powf(-1.0 / 3.0), 0.0, 1.0, 1e-6);
|
||||
assert!(o.error_estimate <= 1e-6,
|
||||
"error_estimate larger then asked. estimate: {:#?}, asked: {:#?}",
|
||||
o.error_estimate,
|
||||
1e-6);
|
||||
assert!((o.integral - 0.41768525592055004).abs() <= o.error_estimate,
|
||||
"error larger then error_estimate");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn demo_function3_works() {
|
||||
let o = integrate(|x| (-x / 5000.0).exp() * (x / 1000.0).powf(-1.0 / 3.0),
|
||||
0.0,
|
||||
10000.0,
|
||||
1e-6);
|
||||
assert!(o.error_estimate <= 1e-6,
|
||||
"error_estimate larger then asked. estimate: {:#?}, asked: {:#?}",
|
||||
o.error_estimate,
|
||||
1e-6);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn demo_bad_function1_works() {
|
||||
let o = integrate(|x| (1.0 - x).powf(0.99), 0.0, 1.0, 1e-6);
|
||||
assert!(o.error_estimate <= 1e-6,
|
||||
"error_estimate larger then asked. estimate: {:#?}, asked: {:#?}",
|
||||
o.error_estimate,
|
||||
1e-6);
|
||||
assert!((o.integral - 0.50251256281407035).abs() <= o.error_estimate,
|
||||
"error larger then error_estimate");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn demo_bad_function2_works() {
|
||||
let o = integrate(|x| x.abs(), -1.0, 1.0, 1e-6);
|
||||
assert!((o.integral - 1.0).abs() <= o.error_estimate,
|
||||
"error larger then error_estimate");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn demo_bad_function3_works() {
|
||||
let o = integrate(|x| (0.5 - x.abs()).abs(), -1.0, 1.0, 1e-6);
|
||||
assert!((o.integral - 0.5).abs() <= o.error_estimate,
|
||||
"error larger then error_estimate");
|
||||
}
|
||||
}
|
||||
358
fractal_dimension/bai_finch_rust/src/integration.rs
Normal file
358
fractal_dimension/bai_finch_rust/src/integration.rs
Normal file
|
|
@ -0,0 +1,358 @@
|
|||
//! The double exponential algorithm is naturally adaptive, it stops calling the integrand when the error is reduced to below the desired threshold.
|
||||
//! It also does not allocate. No box, no vec, etc.
|
||||
//! It has a hard coded maximum of approximately 350 function evaluations. This guarantees that the algorithm will return.
|
||||
//! The error in the algorithm decreases exponentially in the number of function evaluations, specifically O(exp(-cN/log(N))). So if 350 function evaluations is not giving the desired accuracy than the programmer probably needs to give some guidance by splitting up the range at singularities or [other preparation techniques](http://www.johndcook.com/blog/2012/02/21/care-and-treatment-of-singularities/).
|
||||
//!
|
||||
//! This is a port of the [Fast Numerical Integration](https://www.codeproject.com/kb/recipes/fastnumericalintegration.aspx) from c++ to rust. The original code is by John D. Cook, and is licensed under the [BSD](https://opensource.org/licenses/bsd-license.php).
|
||||
|
||||
use num_complex::Complex64;
|
||||
|
||||
/// Integrate an analytic function over a finite interval.
|
||||
/// f is the function to be integrated.
|
||||
/// a is left limit of integration.
|
||||
/// b is right limit of integration
|
||||
/// target_absolute_error is the desired bound on error
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// use quadrature::integrate;
|
||||
/// fn integrand(x: f64) -> f64 {
|
||||
/// (-x / 5.0).exp() * x.powf(-1.0 / 3.0)
|
||||
/// }
|
||||
/// let o = integrate(integrand , 0.0, 10.0, 1e-6);
|
||||
/// assert!((o.integral - 3.6798142583691758).abs() <= 1e-6);
|
||||
/// ```
|
||||
pub fn integrate<F>(f: F, a: f64, b: f64, target_absolute_error: f64) -> Complex64
|
||||
where
|
||||
F: Fn(f64) -> Complex64,
|
||||
{
|
||||
// Apply the linear change of variables x = ct + d
|
||||
// $$\int_a^b f(x) dx = c \int_{-1}^1 f( ct + d ) dt$$
|
||||
// c = (b-a)/2, d = (a+b)/2
|
||||
let c = 0.5 * (b - a);
|
||||
let d = 0.5 * (a + b);
|
||||
integrate_core(
|
||||
|x| {
|
||||
let out = f(c * x + d);
|
||||
if out.is_finite() {
|
||||
out
|
||||
} else {
|
||||
Complex64::new(0.0, 0.0)
|
||||
}
|
||||
},
|
||||
0.25 * target_absolute_error / c,
|
||||
)
|
||||
.scale(c)
|
||||
}
|
||||
|
||||
/// Integrate f(x) from [-1.0, 1.0]
|
||||
fn integrate_core<F>(f: F, target_absolute_error: f64) -> Complex64
|
||||
where
|
||||
F: Fn(f64) -> Complex64,
|
||||
{
|
||||
let mut error_estimate = ::std::f64::MAX;
|
||||
let mut num_function_evaluations = 1;
|
||||
let mut current_delta = ::std::f64::MAX;
|
||||
|
||||
let mut integral = 2.0 * ::std::f64::consts::FRAC_PI_2 * f(0.0);
|
||||
|
||||
for &weight in &WEIGHTS {
|
||||
let new_contribution = weight
|
||||
.iter()
|
||||
.map(|&(w, x)| w * (f(x) + f(-x)))
|
||||
.fold(Complex64::new(0.0, 0.0), |sum, x| sum + x);
|
||||
num_function_evaluations += 2 * weight.len();
|
||||
|
||||
// difference in consecutive integral estimates
|
||||
let previous_delta_ln = current_delta.ln();
|
||||
current_delta = (0.5 * integral - new_contribution).norm();
|
||||
integral = 0.5 * integral + new_contribution;
|
||||
|
||||
// Once convergence kicks in, error is approximately squared at each step.
|
||||
// Determine whether we're in the convergent region by looking at the trend in the error.
|
||||
if num_function_evaluations <= 13 {
|
||||
// level <= 1
|
||||
continue; // previousDelta meaningless, so cannot check convergence.
|
||||
}
|
||||
|
||||
// Exact comparison with zero is harmless here. Could possibly be replaced with
|
||||
// a small positive upper limit on the size of currentDelta, but determining
|
||||
// that upper limit would be difficult. At worse, the loop is executed more
|
||||
// times than necessary. But no infinite loop can result since there is
|
||||
// an upper bound on the loop variable.
|
||||
if current_delta == 0.0 {
|
||||
error_estimate = 0.0;
|
||||
break;
|
||||
}
|
||||
// previousDelta != 0 or would have been kicked out previously
|
||||
let r = current_delta.ln() / previous_delta_ln;
|
||||
|
||||
if r > 1.9 && r < 2.1 {
|
||||
// If convergence theory applied perfectly, r would be 2 in the convergence region.
|
||||
// r close to 2 is good enough. We expect the difference between this integral estimate
|
||||
// and the next one to be roughly delta^2.
|
||||
error_estimate = current_delta * current_delta;
|
||||
} else {
|
||||
// Not in the convergence region. Assume only that error is decreasing.
|
||||
error_estimate = current_delta;
|
||||
}
|
||||
|
||||
if error_estimate < target_absolute_error {
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
integral
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use std::f64::consts::PI;
|
||||
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn trivial_function_works() {
|
||||
let unit_circle =
|
||||
|t: f64| -> Complex64 { Complex64::new((PI * 2.0 * t).cos(), (PI * 2.0 * t).sin()) };
|
||||
let dzdt = |t: f64| {
|
||||
Complex64::new(
|
||||
-2.0 * PI * (2.0 * PI * t).sin(),
|
||||
2.0 * PI * (2.0 * PI * t).cos(),
|
||||
)
|
||||
};
|
||||
let f = |t: f64| dzdt(t) / unit_circle(t);
|
||||
let o = integrate(f, 0.0, 1.0, 1e-14);
|
||||
println!("{}", o);
|
||||
assert!((o - Complex64::new(0.0, 2.0 * PI)).norm() < 1e-14);
|
||||
}
|
||||
}
|
||||
|
||||
pub const WEIGHTS: [&'static [(f64, f64)]; 7] = [
|
||||
&[
|
||||
// First layer weights
|
||||
(0.230022394514788685, 0.95136796407274694573),
|
||||
(0.00026620051375271690866, 0.99997747719246159286),
|
||||
(1.3581784274539090834e-12, 0.99999999999995705839),
|
||||
],
|
||||
&[
|
||||
// 2nd layer weights and abcissas: transformed 1/2, 3/2, 5/2
|
||||
(0.5 * 0.96597657941230114801, 0.67427149224843582608),
|
||||
(0.5 * 0.018343166989927842087, 0.99751485645722438683),
|
||||
(0.5 * 2.1431204556943039358e-7, 0.99999998887566488198),
|
||||
],
|
||||
&[
|
||||
// 3rd layer weights and abcissas: transformed 1/4, 3/4, ...
|
||||
(0.25 * 1.3896147592472563229, 0.37720973816403417379),
|
||||
(0.25 * 0.53107827542805397476, 0.85956905868989663517),
|
||||
(0.25 * 0.076385743570832304188, 0.98704056050737689169),
|
||||
(0.25 * 0.0029025177479013135936, 0.99968826402835320905),
|
||||
(0.25 * 0.000011983701363170720047, 0.99999920473711471266),
|
||||
(0.25 * 1.1631165814255782766e-9, 0.99999999995285644818),
|
||||
],
|
||||
&[
|
||||
// 4th layer weights and abcissas: transformed 1/8, 3/8, ...
|
||||
(0.125 * 1.5232837186347052132, 0.19435700332493543161),
|
||||
(0.125 * 1.1934630258491569639, 0.53914670538796776905),
|
||||
(0.125 * 0.73743784836154784136, 0.78060743898320029925),
|
||||
(0.125 * 0.36046141846934367417, 0.91487926326457461091),
|
||||
(0.125 * 0.13742210773316772341, 0.97396686819567744856),
|
||||
(0.125 * 0.039175005493600779072, 0.99405550663140214329),
|
||||
(0.125 * 0.0077426010260642407123, 0.99906519645578584642),
|
||||
(0.125 * 0.00094994680428346871691, 0.99990938469514399984),
|
||||
(0.125 * 0.000062482559240744082891, 0.99999531604122052843),
|
||||
(0.125 * 1.8263320593710659699e-6, 0.99999989278161241838),
|
||||
(0.125 * 1.8687282268736410132e-8, 0.99999999914270509218),
|
||||
(0.125 * 4.9378538776631926964e-11, 0.99999999999823216531),
|
||||
],
|
||||
&[
|
||||
// 5th layer weights and abcissa: transformed 1/16, 3/16, ...
|
||||
(0.0625 * 1.5587733555333301451, 0.097923885287832333262),
|
||||
(0.0625 * 1.466014426716965781, 0.28787993274271591456),
|
||||
(0.0625 * 1.297475750424977998, 0.46125354393958570440),
|
||||
(0.0625 * 1.0816349854900704074, 0.61027365750063894488),
|
||||
(0.0625 * 0.85017285645662006895, 0.73101803479256151149),
|
||||
(0.0625 * 0.63040513516474369106, 0.82331700550640237006),
|
||||
(0.0625 * 0.44083323627385823707, 0.88989140278426019808),
|
||||
(0.0625 * 0.290240679312454185, 0.93516085752198468323),
|
||||
(0.0625 * 0.17932441211072829296, 0.96411216422354729193),
|
||||
(0.0625 * 0.10343215422333290062, 0.98145482667733517003),
|
||||
(0.0625 * 0.055289683742240583845, 0.99112699244169880223),
|
||||
(0.0625 * 0.027133510013712003219, 0.99610866543750854254),
|
||||
(0.0625 * 0.012083543599157953493, 0.99845420876769773751),
|
||||
(0.0625 * 0.0048162981439284630173, 0.99945143443527460584),
|
||||
(0.0625 * 0.0016908739981426396472, 0.99982882207287494166),
|
||||
(0.0625 * 0.00051339382406790336017, 0.99995387100562796075),
|
||||
(0.0625 * 0.00013205234125609974879, 0.99998948201481850361),
|
||||
(0.0625 * 0.000028110164327940134749, 0.99999801714059543208),
|
||||
(0.0625 * 4.8237182032615502124e-6, 0.99999969889415261122),
|
||||
(0.0625 * 6.4777566035929719908e-7, 0.99999996423908091534),
|
||||
(0.0625 * 6.5835185127183396672e-8, 0.99999999678719909830),
|
||||
(0.0625 * 4.8760060974240625869e-9, 0.99999999978973286224),
|
||||
(0.0625 * 2.5216347918530148572e-10, 0.99999999999039393352),
|
||||
(0.0625 * 8.6759314149796046502e-12, 0.99999999999970809734),
|
||||
],
|
||||
&[
|
||||
// 6th layer weights and abcissas: transformed 1/32, 3/32, ...
|
||||
(0.03125 * 1.5677814313072218572, 0.049055967305077886315),
|
||||
(0.03125 * 1.5438811161769592204, 0.14641798429058794053),
|
||||
(0.03125 * 1.4972262225410362896, 0.24156631953888365838),
|
||||
(0.03125 * 1.4300083548722996676, 0.33314226457763809244),
|
||||
(0.03125 * 1.3452788847662516615, 0.41995211127844715849),
|
||||
(0.03125 * 1.2467012074518577048, 0.50101338937930910152),
|
||||
(0.03125 * 1.1382722433763053734, 0.57558449063515165995),
|
||||
(0.03125 * 1.0240449331118114483, 0.64317675898520470128),
|
||||
(0.03125 * 0.90787937915489531693, 0.70355000514714201566),
|
||||
(0.03125 * 0.79324270082051671787, 0.75669390863372994941),
|
||||
(0.03125 * 0.68306851634426375464, 0.80279874134324126576),
|
||||
(0.03125 * 0.57967810308778764708, 0.84221924635075686382),
|
||||
(0.03125 * 0.48475809121475539287, 0.87543539763040867837),
|
||||
(0.03125 * 0.39938474152571713515, 0.90301328151357387064),
|
||||
(0.03125 * 0.32408253961152890402, 0.92556863406861266645),
|
||||
(0.03125 * 0.258904639514053516, 0.94373478605275715685),
|
||||
(0.03125 * 0.20352399885860174519, 0.95813602271021369012),
|
||||
(0.03125 * 0.15732620348436615027, 0.96936673289691733517),
|
||||
(0.03125 * 0.11949741128869592428, 0.97797623518666497298),
|
||||
(0.03125 * 0.089103139240941462841, 0.98445883116743083087),
|
||||
(0.03125 * 0.065155533432536205042, 0.98924843109013389601),
|
||||
(0.03125 * 0.046668208054846613644, 0.99271699719682728538),
|
||||
(0.03125 * 0.032698732726609031113, 0.99517602615532735426),
|
||||
(0.03125 * 0.022379471063648476483, 0.99688031812819187372),
|
||||
(0.03125 * 0.014937835096050129696, 0.99803333631543375402),
|
||||
(0.03125 * 0.0097072237393916892692, 0.99879353429880589929),
|
||||
(0.03125 * 0.0061300376320830301252, 0.99928111192179195541),
|
||||
(0.03125 * 0.0037542509774318343023, 0.99958475035151758732),
|
||||
(0.03125 * 0.0022250827064786427022, 0.99976797159956083506),
|
||||
(0.03125 * 0.0012733279447082382027, 0.99987486504878034648),
|
||||
(0.03125 * 0.0007018595156842422708, 0.99993501992508242369),
|
||||
(0.03125 * 0.00037166693621677760301, 0.99996759306794345976),
|
||||
(0.03125 * 0.00018856442976700318572, 0.99998451990227082442),
|
||||
(0.03125 * 0.000091390817490710122732, 0.99999293787666288565),
|
||||
(0.03125 * 0.000042183183841757600604, 0.99999693244919035751),
|
||||
(0.03125 * 0.000018481813599879217116, 0.99999873547186590954),
|
||||
(0.03125 * 7.6595758525203162562e-6, 0.99999950700571943689),
|
||||
(0.03125 * 2.9916615878138787094e-6, 0.99999981889371276701),
|
||||
(0.03125 * 1.0968835125901264732e-6, 0.99999993755407837378),
|
||||
(0.03125 * 3.7595411862360630091e-7, 0.99999997987450320175),
|
||||
(0.03125 * 1.1992442782902770219e-7, 0.99999999396413420165),
|
||||
(0.03125 * 3.5434777171421953043e-8, 0.99999999832336194826),
|
||||
(0.03125 * 9.6498888961089633609e-9, 0.99999999957078777261),
|
||||
(0.03125 * 2.4091773256475940779e-9, 0.99999999989927772326),
|
||||
(0.03125 * 5.482835779709497755e-10, 0.99999999997845533741),
|
||||
(0.03125 * 1.1306055347494680536e-10, 0.99999999999582460688),
|
||||
(0.03125 * 2.0989335404511469109e-11, 0.99999999999927152627),
|
||||
(0.03125 * 3.4841937670261059685e-12, 0.99999999999988636130),
|
||||
],
|
||||
&[
|
||||
// 7th layer weights and abcissas: transformed 1/64, 3/64, ...
|
||||
(0.015625 * 1.5700420292795931467, 0.024539763574649160379),
|
||||
(0.015625 * 1.5640214037732320999, 0.073525122985671294475),
|
||||
(0.015625 * 1.5520531698454121192, 0.12222912220155764235),
|
||||
(0.015625 * 1.5342817381543034316, 0.17046797238201051811),
|
||||
(0.015625 * 1.5109197230741697127, 0.21806347346971200463),
|
||||
(0.015625 * 1.48224329788553807, 0.26484507658344795046),
|
||||
(0.015625 * 1.4485862549613225916, 0.31065178055284596083),
|
||||
(0.015625 * 1.4103329714462590129, 0.35533382516507453330),
|
||||
(0.015625 * 1.3679105116808964881, 0.39875415046723775644),
|
||||
(0.015625 * 1.3217801174437728579, 0.44078959903390086627),
|
||||
(0.015625 * 1.2724283455378627082, 0.48133184611690504422),
|
||||
(0.015625 * 1.2203581095793582207, 0.52028805069123015958),
|
||||
(0.015625 * 1.1660798699324345766, 0.55758122826077823080),
|
||||
(0.015625 * 1.1101031939653403796, 0.59315035359195315880),
|
||||
(0.015625 * 1.0529288799552666556, 0.62695020805104287950),
|
||||
(0.015625 * 0.99504180404613271514, 0.65895099174335012438),
|
||||
(0.015625 * 0.93690461274566793366, 0.68913772506166767176),
|
||||
(0.015625 * 0.87895234555278212039, 0.71750946748732412721),
|
||||
(0.015625 * 0.82158803526696470334, 0.74407838354734739913),
|
||||
(0.015625 * 0.7651792989089561367, 0.76886868676824658459),
|
||||
(0.015625 * 0.71005590120546898385, 0.79191549237614211447),
|
||||
(0.015625 * 0.65650824613162753076, 0.81326360850297385168),
|
||||
(0.015625 * 0.60478673057840362158, 0.83296629391941087564),
|
||||
(0.015625 * 0.55510187800363350959, 0.85108400798784873261),
|
||||
(0.015625 * 0.5076251588319080997, 0.86768317577564598669),
|
||||
(0.015625 * 0.4624903980553677613, 0.88283498824466895513),
|
||||
(0.015625 * 0.41979566844501548066, 0.89661425428007602579),
|
||||
(0.015625 * 0.37960556938665160999, 0.90909831816302043511),
|
||||
(0.015625 * 0.3419537959230168323, 0.92036605303195280235),
|
||||
(0.015625 * 0.30684590941791694932, 0.93049693799715340631),
|
||||
(0.015625 * 0.27426222968906810637, 0.93957022393327475539),
|
||||
(0.015625 * 0.24416077786983990868, 0.94766419061515309734),
|
||||
(0.015625 * 0.21648020911729617038, 0.95485549580502268541),
|
||||
(0.015625 * 0.19114268413342749532, 0.96121861515111640753),
|
||||
(0.015625 * 0.16805663794826916233, 0.96682537031235585284),
|
||||
(0.015625 * 0.14711941325785693248, 0.97174454156548730892),
|
||||
(0.015625 * 0.12821973363120098675, 0.97604156025657673933),
|
||||
(0.015625 * 0.11123999898874453035, 0.97977827580061576265),
|
||||
(0.015625 * 0.096058391865189467849, 0.98301279148110110558),
|
||||
(0.015625 * 0.082550788110701737654, 0.98579936302528343597),
|
||||
(0.015625 * 0.070592469906866999352, 0.98818835380074264243),
|
||||
(0.015625 * 0.060059642358636300319, 0.99022624046752774694),
|
||||
(0.015625 * 0.05083075757257047107, 0.99195566300267761562),
|
||||
(0.015625 * 0.042787652157725676034, 0.99341551316926403900),
|
||||
(0.015625 * 0.035816505604196436523, 0.99464105571251119672),
|
||||
(0.015625 * 0.029808628117310126969, 0.99566407681695316965),
|
||||
(0.015625 * 0.024661087314753282511, 0.99651305464025377317),
|
||||
(0.015625 * 0.020277183817500123926, 0.99721334704346870224),
|
||||
(0.015625 * 0.016566786254247575375, 0.99778739195890653083),
|
||||
(0.015625 * 0.013446536605285730674, 0.99825491617199629344),
|
||||
(0.015625 * 0.010839937168255907211, 0.99863314864067747762),
|
||||
(0.015625 * 0.0086773307495391815854, 0.99893703483351217373),
|
||||
(0.015625 * 0.0068957859690660035329, 0.99917944893488591716),
|
||||
(0.015625 * 0.0054388997976239984331, 0.99937140114093768690),
|
||||
(0.015625 * 0.0042565295990178580165, 0.99952223765121720422),
|
||||
(0.015625 * 0.0033044669940348302363, 0.99963983134560036519),
|
||||
(0.015625 * 0.0025440657675291729678, 0.99973076151980848263),
|
||||
(0.015625 * 0.0019418357759843675814, 0.99980048143113838630),
|
||||
(0.015625 * 0.0014690143599429791058, 0.99985347277311141171),
|
||||
(0.015625 * 0.0011011261134519383862, 0.99989338654759256426),
|
||||
(0.015625 * 0.00081754101332469493115, 0.99992317012928932869),
|
||||
(0.015625 * 0.00060103987991147422573, 0.99994518061445869309),
|
||||
(0.015625 * 0.00043739495615911687786, 0.99996128480785666613),
|
||||
(0.015625 * 0.00031497209186021200274, 0.99997294642523223656),
|
||||
(0.015625 * 0.00022435965205008549104, 0.99998130127012072679),
|
||||
(0.015625 * 0.00015802788400701191949, 0.99998722128200062811),
|
||||
(0.015625 * 0.00011002112846666697224, 0.99999136844834487344),
|
||||
(
|
||||
0.015625 * 0.000075683996586201477788,
|
||||
0.99999423962761663478,
|
||||
),
|
||||
(
|
||||
0.015625 * 0.000051421497447658802092,
|
||||
0.99999620334716617675,
|
||||
),
|
||||
(0.015625 * 0.0000344921247593431977, 0.99999752962380516793),
|
||||
(
|
||||
0.015625 * 0.000022832118109036146591,
|
||||
0.99999841381096473542,
|
||||
),
|
||||
(
|
||||
0.015625 * 0.000014908514031870608449,
|
||||
0.99999899541068996962,
|
||||
),
|
||||
(0.015625 * 9.5981941283784710776e-6, 0.99999937270733536947),
|
||||
(0.015625 * 6.0899100320949039256e-6, 0.99999961398855024275),
|
||||
(0.015625 * 3.8061983264644899045e-6, 0.99999976602333243312),
|
||||
(0.015625 * 2.3421667208528096843e-6, 0.99999986037121459941),
|
||||
(0.015625 * 1.4183067155493917523e-6, 0.99999991800479471056),
|
||||
(0.015625 * 8.4473756384859863469e-7, 0.99999995264266446185),
|
||||
(0.015625 * 4.9458288702754198508e-7, 0.99999997311323594362),
|
||||
(0.015625 * 2.8449923659159806339e-7, 0.99999998500307631173),
|
||||
(0.015625 * 1.6069394579076224911e-7, 0.99999999178645609907),
|
||||
(0.015625 * 8.9071395140242387124e-8, 0.99999999558563361584),
|
||||
(0.015625 * 4.8420950198072369669e-8, 0.99999999767323673790),
|
||||
(0.015625 * 2.579956822953589238e-8, 0.99999999879798350040),
|
||||
(0.015625 * 1.3464645522302038796e-8, 0.99999999939177687583),
|
||||
(0.015625 * 6.8784610955899001111e-9, 0.99999999969875436925),
|
||||
(0.015625 * 3.4371856744650090511e-9, 0.99999999985405611550),
|
||||
(0.015625 * 1.6788897682161906807e-9, 0.99999999993088839501),
|
||||
(0.015625 * 8.0099784479729665356e-10, 0.99999999996803321674),
|
||||
(0.015625 * 3.7299501843052790038e-10, 0.99999999998556879008),
|
||||
(0.015625 * 1.6939457789411646876e-10, 0.99999999999364632387),
|
||||
(0.015625 * 7.4967397573818224522e-11, 0.99999999999727404948),
|
||||
(0.015625 * 3.230446433325236576e-11, 0.99999999999886126543),
|
||||
(0.015625 * 1.3542512912336274432e-11, 0.99999999999953722654),
|
||||
(0.015625 * 5.5182369468174885821e-12, 0.99999999999981720098),
|
||||
(0.015625 * 2.1835922099233609052e-12, 0.99999999999992987953),
|
||||
],
|
||||
]; // end weights
|
||||
|
|
@ -1,6 +1,7 @@
|
|||
#![allow(dead_code)]
|
||||
|
||||
mod diff;
|
||||
mod integration;
|
||||
|
||||
use diff::*;
|
||||
use num_complex::Complex64;
|
||||
|
|
@ -8,35 +9,110 @@ use std::f64::consts::PI;
|
|||
// use nalgebra::{SMatrix, SVector};
|
||||
|
||||
type F = f64;
|
||||
// type NComplex = num::Complex<F>;
|
||||
// type Matrix2x2 = SMatrix<num::Complex<F>, 2, 2>;
|
||||
// type Complex = Complex64;
|
||||
type Matrix2x2 = nalgebra::SMatrix<Complex64, 2, 2>;
|
||||
type MatrixBig = nalgebra::SMatrix<Complex64, NC2, NC2>;
|
||||
type MatrixBigr = nalgebra::SMatrix<f64, NC2, NC2>;
|
||||
|
||||
// fn power_method<const N: usize>(
|
||||
// vec: SVector<NComplex, N>,
|
||||
// mat: SMatrix<NComplex, N, N>,
|
||||
// iterations: usize,
|
||||
// ) -> NComplex {
|
||||
// let mut current = vec;
|
||||
// let mut previous = vec;
|
||||
// for _ in 0..iterations {
|
||||
// previous = current;
|
||||
// current = mat * current;
|
||||
// println!("current guess: {}", current[0] / previous[0]);
|
||||
// }
|
||||
const G: Matrix2x2 = Matrix2x2::new(
|
||||
Complex64::new(-2.0 / 12.0, -2.0 / 12.0),
|
||||
Complex64::new(-1.0 / 12.0, -5.0 / 12.0),
|
||||
Complex64::new(4.0 / 12.0, -4.0 / 12.0),
|
||||
Complex64::new(-2.0 / 12.0, -10.0 / 12.0),
|
||||
);
|
||||
|
||||
// current[0] / previous[0]
|
||||
// }
|
||||
const R: Matrix2x2 = Matrix2x2::new(
|
||||
Complex64::new(-6.0 / 12.0, 8.0 / 12.0),
|
||||
Complex64::new(0.0, 11.0 / 12.0),
|
||||
Complex64::new(0.0, 4.0 / 12.0),
|
||||
Complex64::new(-6.0 / 12.0, -8.0 / 12.0),
|
||||
);
|
||||
|
||||
fn secant_method(f: fn(F) -> F, x0: F, x1: F, accuracy: F, iterations: usize) -> F {
|
||||
const NC: usize = 10;
|
||||
const K0: i32 = 100;
|
||||
const LC: usize = 15;
|
||||
const NC2: usize = NC * NC;
|
||||
|
||||
fn fancy_l(q: f64) -> MatrixBigr {
|
||||
MatrixBigr::from_diagonal_element(2.0) * fancy_m(q).map(|z| z.re)
|
||||
}
|
||||
|
||||
fn integral_choose(n: i32, k: i32) -> f64 {
|
||||
match k.cmp(&0) {
|
||||
std::cmp::Ordering::Less => 0.0,
|
||||
std::cmp::Ordering::Equal => 1.0,
|
||||
std::cmp::Ordering::Greater => {
|
||||
(0..k).map(|i|
|
||||
(n - i) as f64 / (k - i) as f64
|
||||
).product()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn non_integral_choose(q: f64, k: i32) -> f64 {
|
||||
match k.cmp(&0) {
|
||||
std::cmp::Ordering::Less => 0.0,
|
||||
std::cmp::Ordering::Equal => 1.0,
|
||||
std::cmp::Ordering::Greater => {
|
||||
(0..k).map(|i|
|
||||
(q - i as f64) / (k - i) as f64
|
||||
).product()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn normal_m(k: i32, q: f64, n: i32, s: i32) -> Complex64 {
|
||||
let ak = R + Matrix2x2::from_diagonal_element(Complex64::new(k as f64, 0.0)) * G;
|
||||
let a11 = ak[(0, 0)];
|
||||
let a12 = ak[(0, 1)];
|
||||
let a21 = ak[(1, 0)];
|
||||
let a22 = ak[(1, 1)];
|
||||
|
||||
(0..=s)
|
||||
.map(|j| {
|
||||
non_integral_choose(-n as f64 - q, j)
|
||||
* integral_choose(n, s - j)
|
||||
* a21.powi(k)
|
||||
* a22.powf(-n as f64 - q - j as f64)
|
||||
* a11.powi(s - j)
|
||||
* a12.powi(n - s + k)
|
||||
})
|
||||
.sum()
|
||||
}
|
||||
|
||||
fn fancy_m(q: f64) -> MatrixBig {
|
||||
let mut fancy_m = MatrixBig::zeros();
|
||||
for m in 0..NC {
|
||||
for n in 0..NC {
|
||||
for r in 0..NC {
|
||||
for s in 0..NC {
|
||||
if m <= n {
|
||||
let mut sum = Complex64::new(0.0, 0.0);
|
||||
for k in 1..K0 {
|
||||
sum += normal_m(k, q, n as i32, s as i32)
|
||||
* normal_m(k as i32, q, m as i32, r as i32).conj();
|
||||
}
|
||||
fancy_m[(m * NC + n, r * NC + s)] = sum;
|
||||
} else {
|
||||
fancy_m[(m * NC + n, r * NC + s)] =
|
||||
fancy_m[(n * NC + m, s * NC + r)].conj();
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
fancy_m
|
||||
}
|
||||
|
||||
fn secant_method(f: fn(F) -> F, target: F, x0: F, x1: F, accuracy: F, iterations: usize) -> F {
|
||||
let mut x0 = x0;
|
||||
let mut x1 = x1;
|
||||
let mut y0 = f(x0);
|
||||
let mut y1 = f(x1);
|
||||
let mut count = 0;
|
||||
|
||||
while y1.abs() >= accuracy && count < iterations {
|
||||
let new_x = x0 - y0 * (x1 - x0) / (y1 - y0);
|
||||
println!("f({}) =\t{}", x1, y1);
|
||||
while (y1 - target).abs() >= accuracy && count < iterations {
|
||||
let new_x = x0 - (y0 - target) * (x1 - x0) / (y1 - y0);
|
||||
x0 = x1;
|
||||
x1 = new_x;
|
||||
y0 = y1;
|
||||
|
|
@ -47,50 +123,38 @@ fn secant_method(f: fn(F) -> F, x0: F, x1: F, accuracy: F, iterations: usize) ->
|
|||
x0
|
||||
}
|
||||
|
||||
fn power_method<const N: usize>(
|
||||
vec: nalgebra::SVector<f64, N>,
|
||||
mat: nalgebra::SMatrix<f64, N, N>,
|
||||
iterations: usize,
|
||||
) -> f64 {
|
||||
let mut current = vec;
|
||||
let mut previous = vec;
|
||||
for _ in 0..iterations {
|
||||
previous = current;
|
||||
current = mat * current;
|
||||
println!("current eigenvalue: {}", current[0] / previous[0]);
|
||||
}
|
||||
|
||||
current[0] / previous[0]
|
||||
}
|
||||
|
||||
fn dzdt(t: f64) -> Complex64 {
|
||||
Complex64::new(
|
||||
-2.0 * PI * (-2.0 * PI * t).sin(),
|
||||
-2.0 * PI * (2.0 * PI * t).sin(),
|
||||
2.0 * PI * (2.0 * PI * t).cos(),
|
||||
)
|
||||
}
|
||||
|
||||
fn main() {
|
||||
// let g = Matrix2x2::new(
|
||||
// Complex::new(2.0 / 12.0, -2.0 / 12.0),
|
||||
// Complex::new(-1.0 / 12.0, -5.0 / 12.0),
|
||||
// Complex::new(4.0 / 12.0, -4.0 / 12.0),
|
||||
// Complex::new(-2.0 / 12.0, -10.0 / 12.0),
|
||||
// );
|
||||
|
||||
// let r = Matrix2x2::new(
|
||||
// Complex::new(-6.0 / 12.0, 8.0 / 12.0),
|
||||
// Complex::new(0.0, 11.0 / 12.0),
|
||||
// Complex::new(0.0, 4.0 / 12.0),
|
||||
// Complex::new(-6.0 / 12.0, -8.0 / 12.0),
|
||||
// );
|
||||
|
||||
println!(
|
||||
"{}",
|
||||
diff(&|z: Complex64| z * z, Complex64::new(0.3, 0.3), 1).unwrap()
|
||||
);
|
||||
|
||||
// println!(
|
||||
// "{}",
|
||||
// bacon_sci::integrate::integrate(
|
||||
// 0.0,
|
||||
// 1.0,
|
||||
// &|t: f64| dzdt(t)
|
||||
// / Complex64::new(
|
||||
// (std::f64::consts::PI * 2.0 * t).cos(),
|
||||
// (std::f64::consts::PI * 2.0 * t).sin()
|
||||
// ),
|
||||
// 0.001
|
||||
// )
|
||||
// .unwrap()
|
||||
// );
|
||||
|
||||
println!(
|
||||
"{}",
|
||||
integrate(&|z: Complex64| 1.0 / z, &unit_circle, &trapezoid)
|
||||
);
|
||||
fn lambda(q: f64) -> f64 {
|
||||
let lq = fancy_l(q);
|
||||
let mut phi0 = nalgebra::SVector::zeros();
|
||||
phi0[0] = 1.0;
|
||||
power_method(phi0, lq, 10)
|
||||
}
|
||||
|
||||
fn main() {
|
||||
println!("{}", non_integral_choose(15.3, 3));
|
||||
println!("{}", normal_m(3, 1.3, 2, 1));
|
||||
// secant_method(lambda, 1.0, 1.3, 1.31, f64::EPSILON, 10);
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in a new issue