// 24_math_science.ruchy - Mathematical and scientific computing
import std::math
import std::stats
import std::linalg
fn main() {
println("=== Mathematical & Scientific Computing ===\n")
// Basic math functions
println("=== Math Functions ===")
let x = 2.5
println(f"x = {x}")
println(f"sqrt({x}) = {sqrt(x)}")
println(f"pow({x}, 3) = {pow(x, 3)}")
println(f"exp({x}) = {exp(x)}")
println(f"ln({x}) = {ln(x)}")
println(f"log10({x}) = {log10(x)}")
println(f"sin({x}) = {sin(x)}")
println(f"cos({x}) = {cos(x)}")
println(f"tan({x}) = {tan(x)}")
println(f"floor({x}) = {floor(x)}")
println(f"ceil({x}) = {ceil(x)}")
println(f"round({x}) = {round(x)}")
println(f"abs(-{x}) = {abs(-x)}")
// Constants
println(f"\nMath constants:")
println(f"PI = {math::PI}")
println(f"E = {math::E}")
println(f"TAU = {math::TAU}")
println(f"PHI = {math::PHI}")
// Complex numbers
println("\n=== Complex Numbers ===")
let z1 = Complex(3, 4) // 3 + 4i
let z2 = Complex(2, -1) // 2 - i
println(f"z1 = {z1}")
println(f"z2 = {z2}")
println(f"z1 + z2 = {z1 + z2}")
println(f"z1 * z2 = {z1 * z2}")
println(f"z1 / z2 = {z1 / z2}")
println(f"|z1| = {z1.abs()}")
println(f"arg(z1) = {z1.arg()}")
println(f"conj(z1) = {z1.conjugate()}")
// Vectors and matrices
println("\n=== Linear Algebra ===")
let vec1 = Vector([1, 2, 3])
let vec2 = Vector([4, 5, 6])
println(f"vec1 = {vec1}")
println(f"vec2 = {vec2}")
println(f"vec1 + vec2 = {vec1 + vec2}")
println(f"vec1 · vec2 = {vec1.dot(vec2)}")
println(f"vec1 × vec2 = {vec1.cross(vec2)}")
println(f"||vec1|| = {vec1.norm()}")
println(f"normalized vec1 = {vec1.normalize()}")
// Matrix operations
let mat_a = Matrix([
[1, 2, 3],
[4, 5, 6],
[7, 8, 9]
])
let mat_b = Matrix([
[9, 8, 7],
[6, 5, 4],
[3, 2, 1]
])
println(f"\nMatrix A:")
println(mat_a.pretty())
println(f"Matrix B:")
println(mat_b.pretty())
println(f"A + B = {mat_a + mat_b}")
println(f"A * B = {mat_a * mat_b}")
println(f"A^T = {mat_a.transpose()}")
println(f"det(A) = {mat_a.determinant()}")
println(f"trace(A) = {mat_a.trace()}")
// Matrix decomposition
let (l, u) = mat_a.lu_decomposition()
println(f"LU decomposition: L = {l}, U = {u}")
let (q, r) = mat_a.qr_decomposition()
println(f"QR decomposition: Q = {q}, R = {r}")
let (eigenvalues, eigenvectors) = mat_a.eigen()
println(f"Eigenvalues: {eigenvalues}")
println(f"Eigenvectors: {eigenvectors}")
// Statistics
println("\n=== Statistics ===")
let data = [2.3, 4.5, 1.8, 3.9, 5.1, 2.7, 4.2, 3.6, 2.9, 4.8]
println(f"Data: {data}")
println(f"Mean: {stats::mean(data)}")
println(f"Median: {stats::median(data)}")
println(f"Mode: {stats::mode(data)}")
println(f"Variance: {stats::variance(data)}")
println(f"Std Dev: {stats::std_dev(data)}")
println(f"Min: {stats::min(data)}")
println(f"Max: {stats::max(data)}")
println(f"Range: {stats::range(data)}")
// Percentiles
println(f"25th percentile: {stats::percentile(data, 25)}")
println(f"75th percentile: {stats::percentile(data, 75)}")
println(f"IQR: {stats::iqr(data)}")
// Correlation and covariance
let x_data = [1, 2, 3, 4, 5]
let y_data = [2, 4, 5, 4, 5]
println(f"\nCorrelation: {stats::correlation(x_data, y_data)}")
println(f"Covariance: {stats::covariance(x_data, y_data)}")
// Probability distributions
println("\n=== Probability Distributions ===")
// Normal distribution
let normal = stats::NormalDist(mean=0, std_dev=1)
println(f"Normal(0,1) PDF at 0: {normal.pdf(0)}")
println(f"Normal(0,1) CDF at 0: {normal.cdf(0)}")
println(f"Normal(0,1) quantile 0.95: {normal.quantile(0.95)}")
// Generate random samples
let samples = normal.sample(100)
println(f"Generated {samples.len()} normal samples")
println(f"Sample mean: {stats::mean(samples)}")
println(f"Sample std: {stats::std_dev(samples)}")
// Other distributions
let poisson = stats::PoissonDist(lambda=3.5)
let binomial = stats::BinomialDist(n=10, p=0.3)
let exponential = stats::ExponentialDist(rate=0.5)
// Hypothesis testing
println("\n=== Hypothesis Testing ===")
let sample1 = [23, 25, 28, 22, 26]
let sample2 = [30, 29, 33, 31, 32]
let t_test = stats::t_test(sample1, sample2)
println(f"T-test statistic: {t_test.statistic}")
println(f"P-value: {t_test.p_value}")
println(f"Significant at 0.05? {t_test.p_value < 0.05}")
// Chi-square test
let observed = [20, 30, 25, 25]
let expected = [25, 25, 25, 25]
let chi2 = stats::chi_square_test(observed, expected)
println(f"Chi-square statistic: {chi2.statistic}")
println(f"P-value: {chi2.p_value}")
// Regression analysis
println("\n=== Regression Analysis ===")
let x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
let y = [2.1, 3.9, 6.1, 7.8, 10.2, 11.9, 14.1, 15.8, 18.2, 19.9]
let regression = stats::linear_regression(x, y)
println(f"Slope: {regression.slope}")
println(f"Intercept: {regression.intercept}")
println(f"R-squared: {regression.r_squared}")
println(f"Equation: y = {regression.slope}x + {regression.intercept}")
// Predictions
let x_new = 11
let y_pred = regression.predict(x_new)
println(f"Prediction for x={x_new}: {y_pred}")
// Numerical integration
println("\n=== Numerical Integration ===")
fn integrate(f, a, b, n = 1000) {
let h = (b - a) / n
let mut sum = 0.0
// Simpson's rule
for i in 0..n {
let x0 = a + i * h
let x1 = x0 + h / 2
let x2 = x0 + h
sum += (f(x0) + 4 * f(x1) + f(x2)) * h / 6
}
sum
}
let result = integrate(x => x * x, 0, 1)
println(f"∫₀¹ x² dx = {result} (expected: 1/3 = {1/3})")
// Differential equations
println("\n=== Differential Equations ===")
fn euler_method(f, y0, t0, t_end, dt) {
let mut t = t0
let mut y = y0
let mut solution = [(t, y)]
while t < t_end {
y = y + dt * f(t, y)
t = t + dt
solution.append((t, y))
}
solution
}
// Solve dy/dt = -2y
let solution = euler_method(
(t, y) => -2 * y,
1.0, // y0
0.0, // t0
2.0, // t_end
0.1 // dt
)
println("Solution to dy/dt = -2y:")
for (t, y) in solution.take(5) {
println(f" t={t:.1}, y={y:.4}")
}
// Fast Fourier Transform
println("\n=== Signal Processing ===")
// Generate signal
let sample_rate = 1000
let duration = 1
let t = [i / sample_rate for i in 0..(sample_rate * duration)]
let signal = t.map(t => sin(2 * math::PI * 50 * t) + 0.5 * sin(2 * math::PI * 120 * t))
// Apply FFT
let frequencies = math::fft(signal)
let power_spectrum = frequencies.map(f => f.abs() * f.abs())
println(f"Dominant frequencies detected")
// Find peaks in power spectrum
let peaks = find_peaks(power_spectrum)
for peak in peaks.take(3) {
println(f" Frequency: {peak.freq} Hz, Power: {peak.power}")
}
}