Scientific Computing

/**
 * OpenBLAS matrix multiplication (GEMM - General Matrix Multiply)
 * @param Order Matrix storage order (CblasRowMajor, CblasColMajor)
 * @param TransA Transpose operation for matrix A
 * @param TransB Transpose operation for matrix B  
 * @param M Number of rows in A and C
 * @param N Number of columns in B and C
 * @param K Number of columns in A and rows in B
 * @param alpha Scalar multiplier for A*B
 * @param A Matrix A data
 * @param lda Leading dimension of A
 * @param B Matrix B data
 * @param ldb Leading dimension of B
 * @param beta Scalar multiplier for C
 * @param C Matrix C data (input/output)
 * @param ldc Leading dimension of C
 */
public static native void cblas_dgemm(int Order, int TransA, int TransB,
    int M, int N, int K, double alpha, DoublePointer A, int lda,
    DoublePointer B, int ldb, double beta, DoublePointer C, int ldc);

/**
 * Single precision matrix multiplication
 */
public static native void cblas_sgemm(int Order, int TransA, int TransB,
    int M, int N, int K, float alpha, FloatPointer A, int lda,
    FloatPointer B, int ldb, float beta, FloatPointer C, int ldc);

/**
 * Vector dot product (double precision)
 * @param N Vector length
 * @param X First vector
 * @param incX Stride of X
 * @param Y Second vector  
 * @param incY Stride of Y
 * @return Dot product result
 */
public static native double cblas_ddot(int N, DoublePointer X, int incX,
    DoublePointer Y, int incY);

/**
 * Vector scaling (SCAL)
 * @param N Vector length
 * @param alpha Scale factor
 * @param X Vector to scale
 * @param incX Vector stride
 */
public static native void cblas_dscal(int N, double alpha, DoublePointer X, int incX);

/**
 * Vector addition (AXPY): Y = alpha*X + Y
 * @param N Vector length
 * @param alpha Scale factor for X
 * @param X Source vector
 * @param incX Stride of X
 * @param Y Destination vector
 * @param incY Stride of Y
 */
public static native void cblas_daxpy(int N, double alpha, DoublePointer X, int incX,
    DoublePointer Y, int incY);

/**
 * Solve linear system AX = B using LU decomposition (LAPACK)
 * @param order Matrix storage order
 * @param n Matrix size (A is n×n)
 * @param nrhs Number of right-hand sides
 * @param a Matrix A (modified on output)
 * @param lda Leading dimension of A
 * @param ipiv Pivot indices
 * @param b Right-hand side matrix B (solution on output)
 * @param ldb Leading dimension of B
 * @return 0 on success
 */
public static native int LAPACKE_dgesv(int order, int n, int nrhs,
    DoublePointer a, int lda, IntPointer ipiv, DoublePointer b, int ldb);

/**
 * Compute eigenvalues and eigenvectors of symmetric matrix
 * @param order Matrix storage order
 * @param jobz Compute eigenvectors ('V') or just eigenvalues ('N')
 * @param uplo Upper ('U') or lower ('L') triangular part
 * @param n Matrix size
 * @param a Symmetric matrix (eigenvectors on output if jobz='V')
 * @param lda Leading dimension of A
 * @param w Eigenvalues (output)
 * @return 0 on success
 */
public static native int LAPACKE_dsyev(int order, byte jobz, byte uplo, int n,
    DoublePointer a, int lda, DoublePointer w);

/**
 * FFTW plan for 1D discrete Fourier transform
 * @param n Transform size
 * @param in Input array
 * @param out Output array
 * @param sign Transform direction (FFTW_FORWARD or FFTW_BACKWARD)
 * @param flags Planning flags (FFTW_ESTIMATE, FFTW_MEASURE, etc.)
 * @return FFTW plan
 */
public static native fftw_plan fftw_plan_dft_1d(int n, DoublePointer in,
    DoublePointer out, int sign, int flags);

/**
 * FFTW plan for 2D discrete Fourier transform
 * @param n0 Size of first dimension
 * @param n1 Size of second dimension
 * @param in Input array
 * @param out Output array
 * @param sign Transform direction
 * @param flags Planning flags
 * @return FFTW plan
 */
public static native fftw_plan fftw_plan_dft_2d(int n0, int n1, DoublePointer in,
    DoublePointer out, int sign, int flags);

/**
 * Execute FFTW plan
 * @param plan Prepared FFTW plan
 */
public static native void fftw_execute(fftw_plan plan);

/**
 * Destroy FFTW plan and free resources
 * @param plan Plan to destroy
 */
public static native void fftw_destroy_plan(fftw_plan plan);

/**
 * Real-to-complex 1D FFT plan
 * @param n Transform size
 * @param in Real input array
 * @param out Complex output array
 * @param flags Planning flags
 * @return FFTW plan
 */
public static native fftw_plan fftw_plan_dft_r2c_1d(int n, DoublePointer in,
    fftw_complex out, int flags);

/**
 * Complex-to-real 1D FFT plan (inverse)
 * @param n Transform size
 * @param in Complex input array
 * @param out Real output array
 * @param flags Planning flags
 * @return FFTW plan
 */
public static native fftw_plan fftw_plan_dft_c2r_1d(int n, fftw_complex in,
    DoublePointer out, int flags);

/**
 * FFTW complex number structure
 */
public class fftw_complex extends DoublePointer {
    public fftw_complex(int size);
    public native double real(int index);
    public native double imag(int index);
    public native fftw_complex real(int index, double value);
    public native fftw_complex imag(int index, double value);
}

/**
 * GSL random number generator
 */
public class gsl_rng extends Pointer {
    /**
     * Allocate random number generator
     * @param type RNG type (gsl_rng_mt19937, etc.)
     * @return RNG instance
     */
    public static native gsl_rng gsl_rng_alloc(gsl_rng_type type);
    
    /**
     * Free random number generator
     * @param r RNG to free
     */
    public static native void gsl_rng_free(gsl_rng r);
    
    /**
     * Set random seed
     * @param r Random number generator
     * @param s Seed value
     */
    public static native void gsl_rng_set(gsl_rng r, long s);
    
    /**
     * Generate uniform random number [0,1)
     * @param r Random number generator
     * @return Random double
     */
    public static native double gsl_rng_uniform(gsl_rng r);
    
    /**
     * Generate uniform random integer [0,n)
     * @param r Random number generator
     * @param n Upper bound (exclusive)
     * @return Random integer
     */
    public static native long gsl_rng_uniform_int(gsl_rng r, long n);
}

/**
 * GSL special functions
 */
public static class GSLSpecialFunctions {
    /**
     * Gamma function
     * @param x Input value
     * @return Gamma(x)
     */
    public static native double gsl_sf_gamma(double x);
    
    /**
     * Bessel function J_n(x)
     * @param n Order
     * @param x Argument
     * @return J_n(x)
     */
    public static native double gsl_sf_bessel_Jn(int n, double x);
    
    /**
     * Error function
     * @param x Input value
     * @return erf(x)
     */
    public static native double gsl_sf_erf(double x);
    
    /**
     * Complementary error function
     * @param x Input value
     * @return erfc(x)
     */
    public static native double gsl_sf_erfc(double x);
    
    /**
     * Exponential integral
     * @param x Input value
     * @return Ei(x)
     */
    public static native double gsl_sf_expint_Ei(double x);
}

/**
 * GSL integration routines
 */
public static class GSLIntegration {
    /**
     * Adaptive integration workspace
     */
    public static class gsl_integration_workspace extends Pointer {
        public static native gsl_integration_workspace gsl_integration_workspace_alloc(int n);
        public static native void gsl_integration_workspace_free(gsl_integration_workspace w);
    }
    
    /**
     * Adaptive integration with error control
     * @param f Function to integrate
     * @param a Lower limit
     * @param b Upper limit
     * @param epsabs Absolute error tolerance
     * @param epsrel Relative error tolerance
     * @param limit Maximum subdivisions
     * @param workspace Integration workspace
     * @param result Integration result (output)
     * @param abserr Absolute error estimate (output)  
     * @return Status code
     */
    public static native int gsl_integration_qags(gsl_function f, double a, double b,
        double epsabs, double epsrel, int limit, gsl_integration_workspace workspace,
        DoublePointer result, DoublePointer abserr);
}

/**
 * GSL optimization routines
 */
public static class GSLOptimization {
    /**
     * Minimization algorithm types
     */
    public static native gsl_min_fminimizer_type gsl_min_fminimizer_brent();
    public static native gsl_min_fminimizer_type gsl_min_fminimizer_golden_section();
    
    /**
     * Function minimizer
     */
    public static class gsl_min_fminimizer extends Pointer {
        public static native gsl_min_fminimizer gsl_min_fminimizer_alloc(gsl_min_fminimizer_type T);
        public static native void gsl_min_fminimizer_free(gsl_min_fminimizer s);
        
        /**
         * Set up minimizer
         * @param s Minimizer state
         * @param f Function to minimize
         * @param x_minimum Initial guess for minimum
         * @param x_lower Lower bound
         * @param x_upper Upper bound
         * @return Status code
         */
        public static native int gsl_min_fminimizer_set(gsl_min_fminimizer s,
            gsl_function f, double x_minimum, double x_lower, double x_upper);
        
        /**
         * Perform one minimization iteration
         * @param s Minimizer state
         * @return Status code
         */
        public static native int gsl_min_fminimizer_iterate(gsl_min_fminimizer s);
        
        /**
         * Get current minimum estimate
         * @param s Minimizer state
         * @return Current minimum x value
         */
        public static native double gsl_min_fminimizer_x_minimum(gsl_min_fminimizer s);
    }
}

/**
 * NumPy array object
 */
public class PyArrayObject extends Pointer {
    /**
     * Get array data pointer
     * @return Pointer to array data
     */
    public native Pointer PyArray_DATA();
    
    /**
     * Get number of dimensions
     * @return Number of dimensions
     */
    public native int PyArray_NDIM();
    
    /**
     * Get array dimensions
     * @return Dimensions array
     */
    public native SizeTPointer PyArray_DIMS();
    
    /**
     * Get array strides
     * @return Strides array
     */
    public native SizeTPointer PyArray_STRIDES();
    
    /**
     * Get array data type
     * @return Array dtype
     */
    public native int PyArray_TYPE();
    
    /**
     * Get total number of elements
     * @return Element count
     */
    public native long PyArray_SIZE();
}

/**
 * NumPy array creation functions
 */
public static class NumPyArrays {
    /**
     * Create new array
     * @param nd Number of dimensions
     * @param dims Dimension sizes
     * @param type_num Data type
     * @param fortran Fortran order flag
     * @return New array object
     */
    public static native PyArrayObject PyArray_SimpleNew(int nd, SizeTPointer dims, int type_num);
    
    /**
     * Create array from existing data
     * @param nd Number of dimensions
     * @param dims Dimension sizes
     * @param type_num Data type
     * @param data Existing data buffer
     * @return Array object wrapping data
     */
    public static native PyArrayObject PyArray_SimpleNewFromData(int nd, SizeTPointer dims,
        int type_num, Pointer data);
    
    /**
     * Create array filled with zeros
     * @param nd Number of dimensions
     * @param dims Dimension sizes
     * @param type_num Data type
     * @param fortran Fortran order flag
     * @return Zero-filled array
     */
    public static native PyArrayObject PyArray_ZEROS(int nd, SizeTPointer dims,
        int type_num, int fortran);
    
    /**
     * Create array filled with ones
     * @param nd Number of dimensions
     * @param dims Dimension sizes  
     * @param type_num Data type
     * @param fortran Fortran order flag
     * @return One-filled array
     */
    public static native PyArrayObject PyArray_ONES(int nd, SizeTPointer dims,
        int type_num, int fortran);
}

/**
 * Intel MKL BLAS routines (optimized implementations)
 */
public static class MKL_BLAS {
    /**
     * MKL matrix multiplication (same interface as CBLAS but optimized)
     */
    public static native void cblas_dgemm(int Order, int TransA, int TransB,
        int M, int N, int K, double alpha, DoublePointer A, int lda,
        DoublePointer B, int ldb, double beta, DoublePointer C, int ldc);
}

/**
 * Intel MKL Vector Math Library (VML)
 */
public static class MKL_VML {
    /**
     * Vector exponential function
     * @param n Vector length
     * @param a Input vector
     * @param y Output vector
     */
    public static native void vdExp(int n, DoublePointer a, DoublePointer y);
    
    /**
     * Vector logarithm function
     * @param n Vector length
     * @param a Input vector
     * @param y Output vector
     */
    public static native void vdLn(int n, DoublePointer a, DoublePointer y);
    
    /**
     * Vector sine function
     * @param n Vector length
     * @param a Input vector
     * @param y Output vector
     */
    public static native void vdSin(int n, DoublePointer a, DoublePointer y);
    
    /**
     * Vector cosine function
     * @param n Vector length
     * @param a Input vector
     * @param y Output vector
     */
    public static native void vdCos(int n, DoublePointer a, DoublePointer y);
    
    /**
     * Vector square root function
     * @param n Vector length
     * @param a Input vector
     * @param y Output vector
     */
    public static native void vdSqrt(int n, DoublePointer a, DoublePointer y);
}

/**
 * Intel MKL Random Number Generators
 */
public static class MKL_VSL {
    /**
     * Initialize random stream
     * @param stream Stream handle (output)
     * @param brng Basic random number generator
     * @param seed Random seed
     * @return Status code
     */
    public static native int vslNewStream(VSLStreamStatePtr stream, int brng, int seed);
    
    /**
     * Generate uniform random numbers
     * @param stream Random stream
     * @param n Number of values to generate
     * @param r Output array
     * @param a Lower bound
     * @param b Upper bound
     * @return Status code
     */
    public static native int vdRngUniform(int method, VSLStreamStatePtr stream,
        int n, DoublePointer r, double a, double b);
    
    /**
     * Generate Gaussian random numbers
     * @param stream Random stream
     * @param n Number of values to generate
     * @param r Output array
     * @param a Mean
     * @param sigma Standard deviation
     * @return Status code
     */
    public static native int vdRngGaussian(int method, VSLStreamStatePtr stream,
        int n, DoublePointer r, double a, double sigma);
}

import org.bytedeco.openblas.*;
import static org.bytedeco.openblas.global.openblas.*;

public class LinearAlgebra {
    static {
        Loader.load(openblas.class);
    }
    
    public static void matrixMultiplication() {
        try (PointerScope scope = new PointerScope()) {
            int M = 3, N = 3, K = 3;
            
            // Create matrices
            double[] aData = {1, 2, 3, 4, 5, 6, 7, 8, 9};
            double[] bData = {9, 8, 7, 6, 5, 4, 3, 2, 1};
            double[] cData = new double[M * N];
            
            DoublePointer A = new DoublePointer(aData);
            DoublePointer B = new DoublePointer(bData);
            DoublePointer C = new DoublePointer(cData);
            
            // Perform matrix multiplication: C = A * B
            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                M, N, K, 1.0, A, K, B, N, 0.0, C, N);
            
            // Print result
            System.out.println("Matrix multiplication result:");
            for (int i = 0; i < M; i++) {
                for (int j = 0; j < N; j++) {
                    System.out.printf("%.1f ", cData[i * N + j]);
                }
                System.out.println();
            }
        }
    }
    
    public static void solveLinearSystem() {
        try (PointerScope scope = new PointerScope()) {
            int n = 3;
            
            // System: Ax = b
            double[] aData = {3, 2, -1, 2, -2, 4, -1, 0.5, -1};  // A matrix
            double[] bData = {1, -2, 0};  // b vector
            
            DoublePointer A = new DoublePointer(aData);
            DoublePointer b = new DoublePointer(bData);
            IntPointer ipiv = new IntPointer(n);
            
            // Solve using LU decomposition
            int info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, A, n, ipiv, b, 1);
            
            if (info == 0) {
                System.out.println("Solution:");
                for (int i = 0; i < n; i++) {
                    System.out.printf("x[%d] = %.3f\n", i, bData[i]);
                }
            } else {
                System.err.println("Failed to solve system, info = " + info);
            }
        }
    }
}

import org.bytedeco.fftw3.*;
import static org.bytedeco.fftw3.global.fftw3.*;

public class FFTExample {
    static {
        Loader.load(fftw3.class);
    }
    
    public static void performFFT() {
        try (PointerScope scope = new PointerScope()) {
            int N = 8;
            
            // Create input signal (real part only)
            double[] inputReal = {1, 1, 1, 1, 0, 0, 0, 0};  // Square wave
            double[] inputImag = new double[N];
            
            // Allocate FFTW arrays (interleaved complex format)
            DoublePointer in = new DoublePointer(2 * N);
            DoublePointer out = new DoublePointer(2 * N);
            
            // Fill input array (real, imag, real, imag, ...)
            for (int i = 0; i < N; i++) {
                in.put(2 * i, inputReal[i]);      // Real part
                in.put(2 * i + 1, inputImag[i]); // Imaginary part
            }
            
            // Create FFT plan
            fftw_plan plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
            
            // Execute FFT
            fftw_execute(plan);
            
            // Print results
            System.out.println("FFT Results:");
            for (int i = 0; i < N; i++) {
                double real = out.get(2 * i);
                double imag = out.get(2 * i + 1);
                double magnitude = Math.sqrt(real * real + imag * imag);
                System.out.printf("Bin %d: %.3f + %.3fi (mag: %.3f)\n", 
                    i, real, imag, magnitude);
            }
            
            // Cleanup
            fftw_destroy_plan(plan);
        }
    }
    
    public static void realFFT() {
        try (PointerScope scope = new PointerScope()) {
            int N = 8;
            double[] signal = {0, 1, 0, -1, 0, 1, 0, -1};  // Sine-like signal
            
            DoublePointer in = new DoublePointer(signal);
            fftw_complex out = new fftw_complex(N/2 + 1);  // N/2+1 complex outputs for real input
            
            // Create real-to-complex FFT plan
            fftw_plan plan = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
            
            fftw_execute(plan);
            
            System.out.println("Real FFT Results:");
            for (int i = 0; i <= N/2; i++) {
                double real = out.real(i);
                double imag = out.imag(i);
                System.out.printf("Frequency %d: %.3f + %.3fi\n", i, real, imag);
            }
            
            fftw_destroy_plan(plan);
        }
    }
}

import org.bytedeco.gsl.*;
import static org.bytedeco.gsl.global.gsl.*;

public class StatisticalAnalysis {
    static {
        Loader.load(gsl.class);
    }
    
    public static void randomNumberGeneration() {
        try (PointerScope scope = new PointerScope()) {
            // Initialize random number generator
            gsl_rng_type rngType = gsl_rng_mt19937();
            gsl_rng rng = gsl_rng_alloc(rngType);
            gsl_rng_set(rng, 12345);  // Set seed
            
            System.out.println("Random numbers:");
            for (int i = 0; i < 10; i++) {
                double uniform = gsl_rng_uniform(rng);
                long uniformInt = gsl_rng_uniform_int(rng, 100);
                
                // Generate normal distribution using Box-Muller
                double normal = gsl_ran_gaussian(rng, 1.0);  // mean=0, sigma=1
                
                System.out.printf("Uniform: %.4f, Int: %2d, Normal: %.4f\n",
                    uniform, uniformInt, normal);
            }
            
            gsl_rng_free(rng);
        }
    }
    
    public static void specialFunctions() {
        try (PointerScope scope = new PointerScope()) {
            double x = 2.5;
            
            System.out.println("Special function values for x = " + x + ":");
            System.out.printf("Gamma(x) = %.6f\n", gsl_sf_gamma(x));
            System.out.printf("erf(x) = %.6f\n", gsl_sf_erf(x));
            System.out.printf("erfc(x) = %.6f\n", gsl_sf_erfc(x));
            
            // Bessel functions
            for (int n = 0; n <= 3; n++) {
                double bessel = gsl_sf_bessel_Jn(n, x);
                System.out.printf("J_%d(%.1f) = %.6f\n", n, x, bessel);
            }
        }
    }
    
    public static void numericalIntegration() {
        try (PointerScope scope = new PointerScope()) {
            // Define function to integrate: f(x) = x^2
            gsl_function f = new gsl_function() {
                @Override
                public double function(double x, Pointer params) {
                    return x * x;
                }
            };
            
            gsl_integration_workspace workspace = 
                gsl_integration_workspace.gsl_integration_workspace_alloc(1000);
            
            DoublePointer result = new DoublePointer(1);
            DoublePointer abserr = new DoublePointer(1);
            
            // Integrate x^2 from 0 to 2 (analytical result = 8/3 ≈ 2.667)
            int status = gsl_integration_qags(f, 0.0, 2.0, 1e-7, 1e-7, 1000,
                workspace, result, abserr);
            
            if (status == GSL_SUCCESS) {
                System.out.printf("Integration result: %.6f ± %.2e\n",
                    result.get(), abserr.get());
                System.out.printf("Analytical result: %.6f\n", 8.0/3.0);
            }
            
            gsl_integration_workspace.gsl_integration_workspace_free(workspace);
        }
    }
}