diff --git a/fractal_dimension/bai_finch_rust/Cargo.lock b/fractal_dimension/bai_finch_rust/Cargo.lock index 1aaf551..14ff1df 100644 --- a/fractal_dimension/bai_finch_rust/Cargo.lock +++ b/fractal_dimension/bai_finch_rust/Cargo.lock @@ -2,15 +2,6 @@ # It is not intended for manual editing. version = 3 -[[package]] -name = "approx" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3f2a05fd1bd10b2527e20a2cd32d8873d115b8b39fe219ee25f42a8aca6ba278" -dependencies = [ - "num-traits", -] - [[package]] name = "approx" version = "0.5.0" @@ -26,28 +17,14 @@ version = "1.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "cdb031dd78e28731d87d56cc8ffef4a8f36ca26c38fe2de700543e627f8a464a" -[[package]] -name = "bacon-sci" -version = "0.11.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e5d568d86ec0ca7d24caa67afff9605bc5bc3cae2642d6670ea5924b839ccabc" -dependencies = [ - "nalgebra 0.24.1", - "num-complex 0.3.1", - "num-traits", - "phf", - "phf_codegen", -] - [[package]] name = "bai_finch" version = "0.1.0" dependencies = [ - "bacon-sci", - "nalgebra 0.27.1", + "nalgebra", "num", "num-complex 0.3.1", - "simba 0.3.1", + "rayon", ] [[package]] @@ -57,47 +34,76 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" [[package]] -name = "generic-array" -version = "0.14.4" +name = "crossbeam-channel" +version = "0.5.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "501466ecc8a30d1d3b7fc9229b122b2ce8ed6e9d9223f1138d4babb253e51817" +checksum = "06ed27e177f16d65f0f0c22a213e17c696ace5dd64b14258b52f9417ccb52db4" dependencies = [ - "typenum", - "version_check", + "cfg-if", + "crossbeam-utils", ] [[package]] -name = "getrandom" -version = "0.1.16" +name = "crossbeam-deque" +version = "0.8.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8fc3cb4d91f53b50155bdcfd23f6a4c39ae1969c2ae85982b135750cccaf5fce" +checksum = "94af6efb46fef72616855b036a624cf27ba656ffc9be1b9a3c931cfc7749a9a9" dependencies = [ "cfg-if", - "libc", - "wasi", + "crossbeam-epoch", + "crossbeam-utils", ] +[[package]] +name = "crossbeam-epoch" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ec02e091aa634e2c3ada4a392989e7c3116673ef0ac5b72232439094d73b7fd" +dependencies = [ + "cfg-if", + "crossbeam-utils", + "lazy_static", + "memoffset", + "scopeguard", +] + +[[package]] +name = "crossbeam-utils" +version = "0.8.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d82cfc11ce7f2c3faef78d8a684447b40d503d9681acebed6cb728d45940c4db" +dependencies = [ + "cfg-if", + "lazy_static", +] + +[[package]] +name = "either" +version = "1.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e78d4f1cc4ae33bbfc157ed5d5a5ef3bc29227303d595861deb238fcec4e9457" + +[[package]] +name = "hermit-abi" +version = "0.1.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "62b467343b94ba476dcb2500d242dadbb39557df889310ac77c5d99100aaac33" +dependencies = [ + "libc", +] + +[[package]] +name = "lazy_static" +version = "1.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" + [[package]] name = "libc" version = "0.2.97" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "12b8adadd720df158f4d70dfe7ccc6adb0472d7c55ca83445f6a5ab3e36f8fb6" -[[package]] -name = "libm" -version = "0.2.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c7d73b3f436185384286bd8098d17ec07c9a7d2388a6599f824d8502b529702a" - -[[package]] -name = "matrixmultiply" -version = "0.2.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "916806ba0031cd542105d916a97c8572e1fa6dd79c9c51e7eb43a09ec2dd84c1" -dependencies = [ - "rawpointer", -] - [[package]] name = "matrixmultiply" version = "0.3.1" @@ -108,21 +114,12 @@ dependencies = [ ] [[package]] -name = "nalgebra" -version = "0.24.1" +name = "memoffset" +version = "0.6.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1a9002895a0de45e3cde58b36d0cf2f83249e7dba4a43ee64dafbf01bfd464ff" +checksum = "59accc507f1338036a0477ef61afdae33cde60840f4dfe481319ce3ad116ddf9" dependencies = [ - "approx 0.4.0", - "generic-array", - "matrixmultiply 0.2.4", - "num-complex 0.3.1", - "num-rational 0.3.2", - "num-traits", - "rand", - "rand_distr", - "simba 0.3.1", - "typenum", + "autocfg", ] [[package]] @@ -131,13 +128,13 @@ 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", + "approx", + "matrixmultiply", "nalgebra-macros", "num-complex 0.4.0", "num-rational 0.4.0", "num-traits", - "simba 0.5.1", + "simba", "typenum", ] @@ -246,7 +243,16 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9a64b1ec5cda2586e284722486d802acf1f7dbdc623e2bfc57e65ca1cd099290" dependencies = [ "autocfg", - "libm", +] + +[[package]] +name = "num_cpus" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "05499f3756671c15885fee9034446956fff3f243d6077b91e5767df161f766b3" +dependencies = [ + "hermit-abi", + "libc", ] [[package]] @@ -255,50 +261,6 @@ version = "1.0.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "acbf547ad0c65e31259204bd90935776d1c693cec2f4ff7abb7a1bbbd40dfe58" -[[package]] -name = "phf" -version = "0.8.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3dfb61232e34fcb633f43d12c58f83c1df82962dcdfa565a4e866ffc17dafe12" -dependencies = [ - "phf_shared", -] - -[[package]] -name = "phf_codegen" -version = "0.8.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cbffee61585b0411840d3ece935cce9cb6321f01c45477d30066498cd5e1a815" -dependencies = [ - "phf_generator", - "phf_shared", -] - -[[package]] -name = "phf_generator" -version = "0.8.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "17367f0cc86f2d25802b2c26ee58a7b23faeccf78a396094c13dced0d0182526" -dependencies = [ - "phf_shared", - "rand", -] - -[[package]] -name = "phf_shared" -version = "0.8.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c00cf8b9eafe68dde5e9eaa2cef8ee84a9336a47d566ec55ca16589633b65af7" -dependencies = [ - "siphasher", -] - -[[package]] -name = "ppv-lite86" -version = "0.2.10" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ac74c624d6b2d21f425f752262f42188365d7b8ff1aff74c82e45136510a4857" - [[package]] name = "proc-macro2" version = "1.0.27" @@ -317,67 +279,6 @@ dependencies = [ "proc-macro2", ] -[[package]] -name = "rand" -version = "0.7.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6a6b1679d49b24bbfe0c803429aa1874472f50d9b363131f0e89fc356b544d03" -dependencies = [ - "getrandom", - "libc", - "rand_chacha", - "rand_core", - "rand_hc", - "rand_pcg", -] - -[[package]] -name = "rand_chacha" -version = "0.2.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f4c8ed856279c9737206bf725bf36935d8666ead7aa69b52be55af369d193402" -dependencies = [ - "ppv-lite86", - "rand_core", -] - -[[package]] -name = "rand_core" -version = "0.5.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "90bde5296fc891b0cef12a6d03ddccc162ce7b2aff54160af9338f8d40df6d19" -dependencies = [ - "getrandom", -] - -[[package]] -name = "rand_distr" -version = "0.3.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c9e9532ada3929fb8b2e9dbe28d1e06c9b2cc65813f074fcb6bd5fbefeff9d56" -dependencies = [ - "num-traits", - "rand", -] - -[[package]] -name = "rand_hc" -version = "0.2.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ca3129af7b92a17112d59ad498c6f81eaf463253766b90396d39ea7a39d6613c" -dependencies = [ - "rand_core", -] - -[[package]] -name = "rand_pcg" -version = "0.2.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "16abd0c1b639e9eb4d7c50c0b8100b0d0f849be2349829c740fe8e6eb4816429" -dependencies = [ - "rand_core", -] - [[package]] name = "rawpointer" version = "0.2.1" @@ -385,35 +286,48 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" [[package]] -name = "simba" -version = "0.3.1" +name = "rayon" +version = "1.5.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "17bfe642b1728a6e89137ad428ef5d4738eca4efaba9590f9e110b8944028621" +checksum = "c06aca804d41dbc8ba42dfd964f0d01334eceb64314b9ecf7c5fad5188a06d90" dependencies = [ - "approx 0.4.0", - "num-complex 0.3.1", - "num-traits", - "paste", + "autocfg", + "crossbeam-deque", + "either", + "rayon-core", ] +[[package]] +name = "rayon-core" +version = "1.9.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d78120e2c850279833f1dd3582f730c4ab53ed95aeaaaa862a2a5c71b1656d8e" +dependencies = [ + "crossbeam-channel", + "crossbeam-deque", + "crossbeam-utils", + "lazy_static", + "num_cpus", +] + +[[package]] +name = "scopeguard" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd" + [[package]] name = "simba" version = "0.5.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8e82063457853d00243beda9952e910b82593e4b07ae9f721b9278a99a0d3d5c" dependencies = [ - "approx 0.5.0", + "approx", "num-complex 0.4.0", "num-traits", "paste", ] -[[package]] -name = "siphasher" -version = "0.3.5" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cbce6d4507c7e4a3962091436e56e95290cb71fa302d0d270e32130b75fbff27" - [[package]] name = "syn" version = "1.0.73" @@ -436,15 +350,3 @@ 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" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5fecdca9a5291cc2b8dcf7dc02453fee791a280f3743cb0905f8822ae463b3fe" - -[[package]] -name = "wasi" -version = "0.9.0+wasi-snapshot-preview1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cccddf32554fecc6acb585f82a32a72e28b48f8c4c1883ddfeeeaa96f7d8e519" diff --git a/fractal_dimension/bai_finch_rust/Cargo.toml b/fractal_dimension/bai_finch_rust/Cargo.toml index f5d367a..36764bb 100644 --- a/fractal_dimension/bai_finch_rust/Cargo.toml +++ b/fractal_dimension/bai_finch_rust/Cargo.toml @@ -8,8 +8,6 @@ edition = "2018" [dependencies] num = "0.3.1" -# nalgebra = "0.24.1" nalgebra = "0.27.1" -bacon-sci = "*" -simba = "0.3.1" num-complex = "0.3.1" +rayon = "1.5.1" diff --git a/fractal_dimension/bai_finch_rust/src/diff.rs b/fractal_dimension/bai_finch_rust/src/diff.rs deleted file mode 100644 index 92d296f..0000000 --- a/fractal_dimension/bai_finch_rust/src/diff.rs +++ /dev/null @@ -1,54 +0,0 @@ -use num_complex::Complex64; -use std::f64::consts::PI; - -pub type CFunction<'a> = &'a dyn Fn(Complex64) -> Complex64; -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 { - Complex64::new((2.0 * PI * t).cos(), (2.0 * PI * t).sin()) -} - -pub fn trapezoid(f: CFunction, a: Complex64, b: Complex64) -> Complex64 { - 0.5 * (b - a) * (f(a) + f(b)) -} - -fn factorial(n: i32) -> i32 { - (1..=n).product() -} - -fn test(t: f64) -> Complex64 { - 1.0 / unit_circle(t) -} - -// pub fn diff(fun: CFunction, a: Complex64, n: i32) -> Result { -// let f = |t: f64| -> Complex64 { -// let z = unit_circle(t); -// fun(z) / (z - a).powi(n + 1) -// }; -// let integral = bacon_sci::integrate::integrate_fixed(0f64, 1f64, f, 10)?; -// Ok(integral * Complex64::new(factorial(n) as f64, 0.0) / Complex64::new(0.0, 2.0 * PI)) -// } - -pub fn diff(fun: CFunction, a: Complex64, n: i32) -> Result { - 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.05; - let mut integral = Complex64::new(0.0, 0.0); - - let mut i = increment; - let mut previous_point = path(0.0); - while i <= 1.0 { - let point = path(i); - integral += rule(f, previous_point, point); - previous_point = point; - - i += increment; - } - - integral -} diff --git a/fractal_dimension/bai_finch_rust/src/integration.rs b/fractal_dimension/bai_finch_rust/src/integration.rs deleted file mode 100644 index a006382..0000000 --- a/fractal_dimension/bai_finch_rust/src/integration.rs +++ /dev/null @@ -1,358 +0,0 @@ -//! 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, 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, 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 diff --git a/fractal_dimension/bai_finch_rust/src/main.rs b/fractal_dimension/bai_finch_rust/src/main.rs index 45ec4a0..9fedbca 100644 --- a/fractal_dimension/bai_finch_rust/src/main.rs +++ b/fractal_dimension/bai_finch_rust/src/main.rs @@ -1,12 +1,7 @@ #![allow(dead_code)] -mod diff; -mod integration; - -use diff::*; use num_complex::Complex64; -use std::f64::consts::PI; -// use nalgebra::{SMatrix, SVector}; +use rayon::prelude::*; type F = f64; type Matrix2x2 = nalgebra::SMatrix; @@ -14,7 +9,7 @@ type MatrixBig = nalgebra::SMatrix; type MatrixBigr = nalgebra::SMatrix; const G: Matrix2x2 = Matrix2x2::new( - Complex64::new(-2.0 / 12.0, -2.0 / 12.0), + 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), @@ -27,39 +22,120 @@ const R: Matrix2x2 = Matrix2x2::new( Complex64::new(-6.0 / 12.0, -8.0 / 12.0), ); -const NC: usize = 10; +const NC: usize = 15; const K0: i32 = 100; -const LC: usize = 15; +const LC: usize = 7; const NC2: usize = NC * NC; fn fancy_l(q: f64) -> MatrixBigr { - MatrixBigr::from_diagonal_element(2.0) * fancy_m(q).map(|z| z.re) + MatrixBigr::from_diagonal_element(2.0) * (fancy_m(q) + fancy_f(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() - } + 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() - } + 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_f(q: f64, n: i32, s: i32, l: i32) -> Complex64 { + let r11 = R[(0, 0)]; + let r12 = R[(0, 1)]; + let r21 = R[(1, 0)]; + let r22 = R[(1, 1)]; + + let g11 = G[(0, 0)]; + let g12 = G[(0, 1)]; + let g21 = G[(1, 0)]; + let g22 = G[(1, 1)]; + (0..=s) + .into_par_iter() + .map(|j: i32| -> Complex64 { + non_integral_choose(-n as f64 - q, j) + * integral_choose(n, s - j) + * (0..=j) + .map(|l1: i32| -> Complex64 { + (0..=(s - j)) + .map(|l3: i32| -> Complex64 { + (0..=(n - s + j)) + .map(|l4: i32| -> Complex64 { + integral_choose(j, l1) + * non_integral_choose( + -n as f64 - q - j as f64, + l - l1 - l3 - l4, + ) + * integral_choose(s - j, l3) + * integral_choose(n - s + j, l4) + * r21.powi(l1) + * r22.powi(l - l1 - l3 - l4) + * r11.powi(l3) + * r12.powi(l4) + * g21.powi(j - l1) + * g22.powf( + -n as f64 - q - j as f64 - l as f64 + + l1 as f64 + + l3 as f64 + + l4 as f64, + ) + * g11.powi(s - j - l3) + * g12.powi(n - s + j - l4) + }) + .sum() + }) + .sum() + }) + .sum::() + }) + .sum() +} + +fn zeta(s: f64, k0: i32) -> f64 { + const UPPER_BOUND: i32 = 40_000; + (k0..=UPPER_BOUND) + .into_par_iter() + .map(|j| (j as f64).powf(-s)) + .sum::() + - (UPPER_BOUND as f64).powf(1.0 - s) / (1.0 - s) +} + +fn fancy_f(q: f64) -> MatrixBig { + let mut fancy_f = 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 { + fancy_f[(m * NC + n, r * NC + s)] = (0..=LC) + .map(|l: usize| -> Complex64 { + zeta(l as f64 + 2.0 * q, K0) + * (0..=l) + .map(|lp| { + normal_f(q, m as i32, r as i32, (l - lp) as i32) + * normal_f(q, n as i32, s as i32, lp as i32).conj() + }) + .sum::() + }) + .sum::(); + } else { + fancy_f[(m * NC + n, r * NC + s)] = + fancy_f[(n * NC + m, s * NC + r)].conj(); + } + } + } + } + } + fancy_f +} + 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)]; @@ -68,37 +144,38 @@ fn normal_m(k: i32, q: f64, n: i32, s: i32) -> Complex64 { 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() + .into_par_iter() + .map(|j| { + non_integral_choose(-n as f64 - q, j) + * integral_choose(n, s - j) + * a21.powi(j) + * a22.powf(-n as f64 - q - j as f64) + * a11.powi(s - j) + * a12.powi(n - s + j) + }) + .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(); - } - } - } - } + 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 } @@ -112,13 +189,13 @@ fn secant_method(f: fn(F) -> F, target: F, x0: F, x1: F, accuracy: F, iterations 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; - y1 = f(x1); - println!("f({}) =\t{}", x1, y1); - count += 1; + let new_x = x0 - (y0 - target) * (x1 - x0) / (y1 - y0); + x0 = x1; + x1 = new_x; + y0 = y1; + y1 = f(x1); + println!("f({}) =\t{}", x1, y1); + count += 1; } x0 } @@ -127,34 +204,32 @@ fn power_method( vec: nalgebra::SVector, mat: nalgebra::SMatrix, iterations: usize, + tolerance: f64, ) -> 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]); + let mut previous_val = 0f64; + let mut current = mat * vec; + let mut current_val = current[0] / previous[0]; + let mut count = 0; + while count < iterations && (current_val - previous_val).abs() > tolerance { + previous_val = current_val; + previous = current; + current = mat * current; + current_val = current[0] / previous[0]; + count += 1; } - 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).cos(), - ) + println!("found eigenvalue {} in {} iterations", current_val, count); + current_val } 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) + power_method(phi0, lq, 50, std::f64::EPSILON) } 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); + secant_method(lambda, 1.0, 1.3, 1.31, f64::EPSILON, 10); }