// 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}")
    }
}