0
# Scientific Computing
1

2
Mathematical and scientific computing capabilities through OpenBLAS, Intel MKL, NumPy, SciPy, FFTW, GSL, and specialized numerical libraries.
3

4
## Capabilities
5

6
### Linear Algebra Operations
7

8
High-performance linear algebra through OpenBLAS and Intel MKL implementations.
9

10
```java { .api }
11
/**
12
 * OpenBLAS matrix multiplication (GEMM - General Matrix Multiply)
13
 * @param Order Matrix storage order (CblasRowMajor, CblasColMajor)
14
 * @param TransA Transpose operation for matrix A
15
 * @param TransB Transpose operation for matrix B  
16
 * @param M Number of rows in A and C
17
 * @param N Number of columns in B and C
18
 * @param K Number of columns in A and rows in B
19
 * @param alpha Scalar multiplier for A*B
20
 * @param A Matrix A data
21
 * @param lda Leading dimension of A
22
 * @param B Matrix B data
23
 * @param ldb Leading dimension of B
24
 * @param beta Scalar multiplier for C
25
 * @param C Matrix C data (input/output)
26
 * @param ldc Leading dimension of C
27
 */
28
public static native void cblas_dgemm(int Order, int TransA, int TransB,
29
    int M, int N, int K, double alpha, DoublePointer A, int lda,
30
    DoublePointer B, int ldb, double beta, DoublePointer C, int ldc);
31

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

39
/**
40
 * Vector dot product (double precision)
41
 * @param N Vector length
42
 * @param X First vector
43
 * @param incX Stride of X
44
 * @param Y Second vector  
45
 * @param incY Stride of Y
46
 * @return Dot product result
47
 */
48
public static native double cblas_ddot(int N, DoublePointer X, int incX,
49
    DoublePointer Y, int incY);
50

51
/**
52
 * Vector scaling (SCAL)
53
 * @param N Vector length
54
 * @param alpha Scale factor
55
 * @param X Vector to scale
56
 * @param incX Vector stride
57
 */
58
public static native void cblas_dscal(int N, double alpha, DoublePointer X, int incX);
59

60
/**
61
 * Vector addition (AXPY): Y = alpha*X + Y
62
 * @param N Vector length
63
 * @param alpha Scale factor for X
64
 * @param X Source vector
65
 * @param incX Stride of X
66
 * @param Y Destination vector
67
 * @param incY Stride of Y
68
 */
69
public static native void cblas_daxpy(int N, double alpha, DoublePointer X, int incX,
70
    DoublePointer Y, int incY);
71

72
/**
73
 * Solve linear system AX = B using LU decomposition (LAPACK)
74
 * @param order Matrix storage order
75
 * @param n Matrix size (A is n×n)
76
 * @param nrhs Number of right-hand sides
77
 * @param a Matrix A (modified on output)
78
 * @param lda Leading dimension of A
79
 * @param ipiv Pivot indices
80
 * @param b Right-hand side matrix B (solution on output)
81
 * @param ldb Leading dimension of B
82
 * @return 0 on success
83
 */
84
public static native int LAPACKE_dgesv(int order, int n, int nrhs,
85
    DoublePointer a, int lda, IntPointer ipiv, DoublePointer b, int ldb);
86

87
/**
88
 * Compute eigenvalues and eigenvectors of symmetric matrix
89
 * @param order Matrix storage order
90
 * @param jobz Compute eigenvectors ('V') or just eigenvalues ('N')
91
 * @param uplo Upper ('U') or lower ('L') triangular part
92
 * @param n Matrix size
93
 * @param a Symmetric matrix (eigenvectors on output if jobz='V')
94
 * @param lda Leading dimension of A
95
 * @param w Eigenvalues (output)
96
 * @return 0 on success
97
 */
98
public static native int LAPACKE_dsyev(int order, byte jobz, byte uplo, int n,
99
    DoublePointer a, int lda, DoublePointer w);
100
```
101

102
### Fast Fourier Transform
103

104
High-performance FFT operations through FFTW library.
105

106
```java { .api }
107
/**
108
 * FFTW plan for 1D discrete Fourier transform
109
 * @param n Transform size
110
 * @param in Input array
111
 * @param out Output array
112
 * @param sign Transform direction (FFTW_FORWARD or FFTW_BACKWARD)
113
 * @param flags Planning flags (FFTW_ESTIMATE, FFTW_MEASURE, etc.)
114
 * @return FFTW plan
115
 */
116
public static native fftw_plan fftw_plan_dft_1d(int n, DoublePointer in,
117
    DoublePointer out, int sign, int flags);
118

119
/**
120
 * FFTW plan for 2D discrete Fourier transform
121
 * @param n0 Size of first dimension
122
 * @param n1 Size of second dimension
123
 * @param in Input array
124
 * @param out Output array
125
 * @param sign Transform direction
126
 * @param flags Planning flags
127
 * @return FFTW plan
128
 */
129
public static native fftw_plan fftw_plan_dft_2d(int n0, int n1, DoublePointer in,
130
    DoublePointer out, int sign, int flags);
131

132
/**
133
 * Execute FFTW plan
134
 * @param plan Prepared FFTW plan
135
 */
136
public static native void fftw_execute(fftw_plan plan);
137

138
/**
139
 * Destroy FFTW plan and free resources
140
 * @param plan Plan to destroy
141
 */
142
public static native void fftw_destroy_plan(fftw_plan plan);
143

144
/**
145
 * Real-to-complex 1D FFT plan
146
 * @param n Transform size
147
 * @param in Real input array
148
 * @param out Complex output array
149
 * @param flags Planning flags
150
 * @return FFTW plan
151
 */
152
public static native fftw_plan fftw_plan_dft_r2c_1d(int n, DoublePointer in,
153
    fftw_complex out, int flags);
154

155
/**
156
 * Complex-to-real 1D FFT plan (inverse)
157
 * @param n Transform size
158
 * @param in Complex input array
159
 * @param out Real output array
160
 * @param flags Planning flags
161
 * @return FFTW plan
162
 */
163
public static native fftw_plan fftw_plan_dft_c2r_1d(int n, fftw_complex in,
164
    DoublePointer out, int flags);
165

166
/**
167
 * FFTW complex number structure
168
 */
169
public class fftw_complex extends DoublePointer {
170
    public fftw_complex(int size);
171
    public native double real(int index);
172
    public native double imag(int index);
173
    public native fftw_complex real(int index, double value);
174
    public native fftw_complex imag(int index, double value);
175
}
176
```
177

178
### GNU Scientific Library (GSL)
179

180
Comprehensive collection of mathematical functions and algorithms.
181

182
```java { .api }
183
/**
184
 * GSL random number generator
185
 */
186
public class gsl_rng extends Pointer {
187
    /**
188
     * Allocate random number generator
189
     * @param type RNG type (gsl_rng_mt19937, etc.)
190
     * @return RNG instance
191
     */
192
    public static native gsl_rng gsl_rng_alloc(gsl_rng_type type);
193
    
194
    /**
195
     * Free random number generator
196
     * @param r RNG to free
197
     */
198
    public static native void gsl_rng_free(gsl_rng r);
199
    
200
    /**
201
     * Set random seed
202
     * @param r Random number generator
203
     * @param s Seed value
204
     */
205
    public static native void gsl_rng_set(gsl_rng r, long s);
206
    
207
    /**
208
     * Generate uniform random number [0,1)
209
     * @param r Random number generator
210
     * @return Random double
211
     */
212
    public static native double gsl_rng_uniform(gsl_rng r);
213
    
214
    /**
215
     * Generate uniform random integer [0,n)
216
     * @param r Random number generator
217
     * @param n Upper bound (exclusive)
218
     * @return Random integer
219
     */
220
    public static native long gsl_rng_uniform_int(gsl_rng r, long n);
221
}
222

223
/**
224
 * GSL special functions
225
 */
226
public static class GSLSpecialFunctions {
227
    /**
228
     * Gamma function
229
     * @param x Input value
230
     * @return Gamma(x)
231
     */
232
    public static native double gsl_sf_gamma(double x);
233
    
234
    /**
235
     * Bessel function J_n(x)
236
     * @param n Order
237
     * @param x Argument
238
     * @return J_n(x)
239
     */
240
    public static native double gsl_sf_bessel_Jn(int n, double x);
241
    
242
    /**
243
     * Error function
244
     * @param x Input value
245
     * @return erf(x)
246
     */
247
    public static native double gsl_sf_erf(double x);
248
    
249
    /**
250
     * Complementary error function
251
     * @param x Input value
252
     * @return erfc(x)
253
     */
254
    public static native double gsl_sf_erfc(double x);
255
    
256
    /**
257
     * Exponential integral
258
     * @param x Input value
259
     * @return Ei(x)
260
     */
261
    public static native double gsl_sf_expint_Ei(double x);
262
}
263

264
/**
265
 * GSL integration routines
266
 */
267
public static class GSLIntegration {
268
    /**
269
     * Adaptive integration workspace
270
     */
271
    public static class gsl_integration_workspace extends Pointer {
272
        public static native gsl_integration_workspace gsl_integration_workspace_alloc(int n);
273
        public static native void gsl_integration_workspace_free(gsl_integration_workspace w);
274
    }
275
    
276
    /**
277
     * Adaptive integration with error control
278
     * @param f Function to integrate
279
     * @param a Lower limit
280
     * @param b Upper limit
281
     * @param epsabs Absolute error tolerance
282
     * @param epsrel Relative error tolerance
283
     * @param limit Maximum subdivisions
284
     * @param workspace Integration workspace
285
     * @param result Integration result (output)
286
     * @param abserr Absolute error estimate (output)  
287
     * @return Status code
288
     */
289
    public static native int gsl_integration_qags(gsl_function f, double a, double b,
290
        double epsabs, double epsrel, int limit, gsl_integration_workspace workspace,
291
        DoublePointer result, DoublePointer abserr);
292
}
293

294
/**
295
 * GSL optimization routines
296
 */
297
public static class GSLOptimization {
298
    /**
299
     * Minimization algorithm types
300
     */
301
    public static native gsl_min_fminimizer_type gsl_min_fminimizer_brent();
302
    public static native gsl_min_fminimizer_type gsl_min_fminimizer_golden_section();
303
    
304
    /**
305
     * Function minimizer
306
     */
307
    public static class gsl_min_fminimizer extends Pointer {
308
        public static native gsl_min_fminimizer gsl_min_fminimizer_alloc(gsl_min_fminimizer_type T);
309
        public static native void gsl_min_fminimizer_free(gsl_min_fminimizer s);
310
        
311
        /**
312
         * Set up minimizer
313
         * @param s Minimizer state
314
         * @param f Function to minimize
315
         * @param x_minimum Initial guess for minimum
316
         * @param x_lower Lower bound
317
         * @param x_upper Upper bound
318
         * @return Status code
319
         */
320
        public static native int gsl_min_fminimizer_set(gsl_min_fminimizer s,
321
            gsl_function f, double x_minimum, double x_lower, double x_upper);
322
        
323
        /**
324
         * Perform one minimization iteration
325
         * @param s Minimizer state
326
         * @return Status code
327
         */
328
        public static native int gsl_min_fminimizer_iterate(gsl_min_fminimizer s);
329
        
330
        /**
331
         * Get current minimum estimate
332
         * @param s Minimizer state
333
         * @return Current minimum x value
334
         */
335
        public static native double gsl_min_fminimizer_x_minimum(gsl_min_fminimizer s);
336
    }
337
}
338
```
339

340
### NumPy Integration
341

342
Multi-dimensional array operations and mathematical functions.
343

344
```java { .api }
345
/**
346
 * NumPy array object
347
 */
348
public class PyArrayObject extends Pointer {
349
    /**
350
     * Get array data pointer
351
     * @return Pointer to array data
352
     */
353
    public native Pointer PyArray_DATA();
354
    
355
    /**
356
     * Get number of dimensions
357
     * @return Number of dimensions
358
     */
359
    public native int PyArray_NDIM();
360
    
361
    /**
362
     * Get array dimensions
363
     * @return Dimensions array
364
     */
365
    public native SizeTPointer PyArray_DIMS();
366
    
367
    /**
368
     * Get array strides
369
     * @return Strides array
370
     */
371
    public native SizeTPointer PyArray_STRIDES();
372
    
373
    /**
374
     * Get array data type
375
     * @return Array dtype
376
     */
377
    public native int PyArray_TYPE();
378
    
379
    /**
380
     * Get total number of elements
381
     * @return Element count
382
     */
383
    public native long PyArray_SIZE();
384
}
385

386
/**
387
 * NumPy array creation functions
388
 */
389
public static class NumPyArrays {
390
    /**
391
     * Create new array
392
     * @param nd Number of dimensions
393
     * @param dims Dimension sizes
394
     * @param type_num Data type
395
     * @param fortran Fortran order flag
396
     * @return New array object
397
     */
398
    public static native PyArrayObject PyArray_SimpleNew(int nd, SizeTPointer dims, int type_num);
399
    
400
    /**
401
     * Create array from existing data
402
     * @param nd Number of dimensions
403
     * @param dims Dimension sizes
404
     * @param type_num Data type
405
     * @param data Existing data buffer
406
     * @return Array object wrapping data
407
     */
408
    public static native PyArrayObject PyArray_SimpleNewFromData(int nd, SizeTPointer dims,
409
        int type_num, Pointer data);
410
    
411
    /**
412
     * Create array filled with zeros
413
     * @param nd Number of dimensions
414
     * @param dims Dimension sizes
415
     * @param type_num Data type
416
     * @param fortran Fortran order flag
417
     * @return Zero-filled array
418
     */
419
    public static native PyArrayObject PyArray_ZEROS(int nd, SizeTPointer dims,
420
        int type_num, int fortran);
421
    
422
    /**
423
     * Create array filled with ones
424
     * @param nd Number of dimensions
425
     * @param dims Dimension sizes  
426
     * @param type_num Data type
427
     * @param fortran Fortran order flag
428
     * @return One-filled array
429
     */
430
    public static native PyArrayObject PyArray_ONES(int nd, SizeTPointer dims,
431
        int type_num, int fortran);
432
}
433
```
434

435
### Intel MKL Integration
436

437
Optimized mathematical functions and linear algebra for Intel processors.
438

439
```java { .api }
440
/**
441
 * Intel MKL BLAS routines (optimized implementations)
442
 */
443
public static class MKL_BLAS {
444
    /**
445
     * MKL matrix multiplication (same interface as CBLAS but optimized)
446
     */
447
    public static native void cblas_dgemm(int Order, int TransA, int TransB,
448
        int M, int N, int K, double alpha, DoublePointer A, int lda,
449
        DoublePointer B, int ldb, double beta, DoublePointer C, int ldc);
450
}
451

452
/**
453
 * Intel MKL Vector Math Library (VML)
454
 */
455
public static class MKL_VML {
456
    /**
457
     * Vector exponential function
458
     * @param n Vector length
459
     * @param a Input vector
460
     * @param y Output vector
461
     */
462
    public static native void vdExp(int n, DoublePointer a, DoublePointer y);
463
    
464
    /**
465
     * Vector logarithm function
466
     * @param n Vector length
467
     * @param a Input vector
468
     * @param y Output vector
469
     */
470
    public static native void vdLn(int n, DoublePointer a, DoublePointer y);
471
    
472
    /**
473
     * Vector sine function
474
     * @param n Vector length
475
     * @param a Input vector
476
     * @param y Output vector
477
     */
478
    public static native void vdSin(int n, DoublePointer a, DoublePointer y);
479
    
480
    /**
481
     * Vector cosine function
482
     * @param n Vector length
483
     * @param a Input vector
484
     * @param y Output vector
485
     */
486
    public static native void vdCos(int n, DoublePointer a, DoublePointer y);
487
    
488
    /**
489
     * Vector square root function
490
     * @param n Vector length
491
     * @param a Input vector
492
     * @param y Output vector
493
     */
494
    public static native void vdSqrt(int n, DoublePointer a, DoublePointer y);
495
}
496

497
/**
498
 * Intel MKL Random Number Generators
499
 */
500
public static class MKL_VSL {
501
    /**
502
     * Initialize random stream
503
     * @param stream Stream handle (output)
504
     * @param brng Basic random number generator
505
     * @param seed Random seed
506
     * @return Status code
507
     */
508
    public static native int vslNewStream(VSLStreamStatePtr stream, int brng, int seed);
509
    
510
    /**
511
     * Generate uniform random numbers
512
     * @param stream Random stream
513
     * @param n Number of values to generate
514
     * @param r Output array
515
     * @param a Lower bound
516
     * @param b Upper bound
517
     * @return Status code
518
     */
519
    public static native int vdRngUniform(int method, VSLStreamStatePtr stream,
520
        int n, DoublePointer r, double a, double b);
521
    
522
    /**
523
     * Generate Gaussian random numbers
524
     * @param stream Random stream
525
     * @param n Number of values to generate
526
     * @param r Output array
527
     * @param a Mean
528
     * @param sigma Standard deviation
529
     * @return Status code
530
     */
531
    public static native int vdRngGaussian(int method, VSLStreamStatePtr stream,
532
        int n, DoublePointer r, double a, double sigma);
533
}
534
```
535

536
## Usage Examples
537

538
### Matrix Operations with OpenBLAS
539

540
```java
541
import org.bytedeco.openblas.*;
542
import static org.bytedeco.openblas.global.openblas.*;
543

544
public class LinearAlgebra {
545
    static {
546
        Loader.load(openblas.class);
547
    }
548
    
549
    public static void matrixMultiplication() {
550
        try (PointerScope scope = new PointerScope()) {
551
            int M = 3, N = 3, K = 3;
552
            
553
            // Create matrices
554
            double[] aData = {1, 2, 3, 4, 5, 6, 7, 8, 9};
555
            double[] bData = {9, 8, 7, 6, 5, 4, 3, 2, 1};
556
            double[] cData = new double[M * N];
557
            
558
            DoublePointer A = new DoublePointer(aData);
559
            DoublePointer B = new DoublePointer(bData);
560
            DoublePointer C = new DoublePointer(cData);
561
            
562
            // Perform matrix multiplication: C = A * B
563
            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
564
                M, N, K, 1.0, A, K, B, N, 0.0, C, N);
565
            
566
            // Print result
567
            System.out.println("Matrix multiplication result:");
568
            for (int i = 0; i < M; i++) {
569
                for (int j = 0; j < N; j++) {
570
                    System.out.printf("%.1f ", cData[i * N + j]);
571
                }
572
                System.out.println();
573
            }
574
        }
575
    }
576
    
577
    public static void solveLinearSystem() {
578
        try (PointerScope scope = new PointerScope()) {
579
            int n = 3;
580
            
581
            // System: Ax = b
582
            double[] aData = {3, 2, -1, 2, -2, 4, -1, 0.5, -1};  // A matrix
583
            double[] bData = {1, -2, 0};  // b vector
584
            
585
            DoublePointer A = new DoublePointer(aData);
586
            DoublePointer b = new DoublePointer(bData);
587
            IntPointer ipiv = new IntPointer(n);
588
            
589
            // Solve using LU decomposition
590
            int info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, A, n, ipiv, b, 1);
591
            
592
            if (info == 0) {
593
                System.out.println("Solution:");
594
                for (int i = 0; i < n; i++) {
595
                    System.out.printf("x[%d] = %.3f\n", i, bData[i]);
596
                }
597
            } else {
598
                System.err.println("Failed to solve system, info = " + info);
599
            }
600
        }
601
    }
602
}
603
```
604

605
### FFT with FFTW
606

607
```java
608
import org.bytedeco.fftw3.*;
609
import static org.bytedeco.fftw3.global.fftw3.*;
610

611
public class FFTExample {
612
    static {
613
        Loader.load(fftw3.class);
614
    }
615
    
616
    public static void performFFT() {
617
        try (PointerScope scope = new PointerScope()) {
618
            int N = 8;
619
            
620
            // Create input signal (real part only)
621
            double[] inputReal = {1, 1, 1, 1, 0, 0, 0, 0};  // Square wave
622
            double[] inputImag = new double[N];
623
            
624
            // Allocate FFTW arrays (interleaved complex format)
625
            DoublePointer in = new DoublePointer(2 * N);
626
            DoublePointer out = new DoublePointer(2 * N);
627
            
628
            // Fill input array (real, imag, real, imag, ...)
629
            for (int i = 0; i < N; i++) {
630
                in.put(2 * i, inputReal[i]);      // Real part
631
                in.put(2 * i + 1, inputImag[i]); // Imaginary part
632
            }
633
            
634
            // Create FFT plan
635
            fftw_plan plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
636
            
637
            // Execute FFT
638
            fftw_execute(plan);
639
            
640
            // Print results
641
            System.out.println("FFT Results:");
642
            for (int i = 0; i < N; i++) {
643
                double real = out.get(2 * i);
644
                double imag = out.get(2 * i + 1);
645
                double magnitude = Math.sqrt(real * real + imag * imag);
646
                System.out.printf("Bin %d: %.3f + %.3fi (mag: %.3f)\n", 
647
                    i, real, imag, magnitude);
648
            }
649
            
650
            // Cleanup
651
            fftw_destroy_plan(plan);
652
        }
653
    }
654
    
655
    public static void realFFT() {
656
        try (PointerScope scope = new PointerScope()) {
657
            int N = 8;
658
            double[] signal = {0, 1, 0, -1, 0, 1, 0, -1};  // Sine-like signal
659
            
660
            DoublePointer in = new DoublePointer(signal);
661
            fftw_complex out = new fftw_complex(N/2 + 1);  // N/2+1 complex outputs for real input
662
            
663
            // Create real-to-complex FFT plan
664
            fftw_plan plan = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
665
            
666
            fftw_execute(plan);
667
            
668
            System.out.println("Real FFT Results:");
669
            for (int i = 0; i <= N/2; i++) {
670
                double real = out.real(i);
671
                double imag = out.imag(i);
672
                System.out.printf("Frequency %d: %.3f + %.3fi\n", i, real, imag);
673
            }
674
            
675
            fftw_destroy_plan(plan);
676
        }
677
    }
678
}
679
```
680

681
### GSL Statistical Functions
682

683
```java
684
import org.bytedeco.gsl.*;
685
import static org.bytedeco.gsl.global.gsl.*;
686

687
public class StatisticalAnalysis {
688
    static {
689
        Loader.load(gsl.class);
690
    }
691
    
692
    public static void randomNumberGeneration() {
693
        try (PointerScope scope = new PointerScope()) {
694
            // Initialize random number generator
695
            gsl_rng_type rngType = gsl_rng_mt19937();
696
            gsl_rng rng = gsl_rng_alloc(rngType);
697
            gsl_rng_set(rng, 12345);  // Set seed
698
            
699
            System.out.println("Random numbers:");
700
            for (int i = 0; i < 10; i++) {
701
                double uniform = gsl_rng_uniform(rng);
702
                long uniformInt = gsl_rng_uniform_int(rng, 100);
703
                
704
                // Generate normal distribution using Box-Muller
705
                double normal = gsl_ran_gaussian(rng, 1.0);  // mean=0, sigma=1
706
                
707
                System.out.printf("Uniform: %.4f, Int: %2d, Normal: %.4f\n",
708
                    uniform, uniformInt, normal);
709
            }
710
            
711
            gsl_rng_free(rng);
712
        }
713
    }
714
    
715
    public static void specialFunctions() {
716
        try (PointerScope scope = new PointerScope()) {
717
            double x = 2.5;
718
            
719
            System.out.println("Special function values for x = " + x + ":");
720
            System.out.printf("Gamma(x) = %.6f\n", gsl_sf_gamma(x));
721
            System.out.printf("erf(x) = %.6f\n", gsl_sf_erf(x));
722
            System.out.printf("erfc(x) = %.6f\n", gsl_sf_erfc(x));
723
            
724
            // Bessel functions
725
            for (int n = 0; n <= 3; n++) {
726
                double bessel = gsl_sf_bessel_Jn(n, x);
727
                System.out.printf("J_%d(%.1f) = %.6f\n", n, x, bessel);
728
            }
729
        }
730
    }
731
    
732
    public static void numericalIntegration() {
733
        try (PointerScope scope = new PointerScope()) {
734
            // Define function to integrate: f(x) = x^2
735
            gsl_function f = new gsl_function() {
736
                @Override
737
                public double function(double x, Pointer params) {
738
                    return x * x;
739
                }
740
            };
741
            
742
            gsl_integration_workspace workspace = 
743
                gsl_integration_workspace.gsl_integration_workspace_alloc(1000);
744
            
745
            DoublePointer result = new DoublePointer(1);
746
            DoublePointer abserr = new DoublePointer(1);
747
            
748
            // Integrate x^2 from 0 to 2 (analytical result = 8/3 ≈ 2.667)
749
            int status = gsl_integration_qags(f, 0.0, 2.0, 1e-7, 1e-7, 1000,
750
                workspace, result, abserr);
751
            
752
            if (status == GSL_SUCCESS) {
753
                System.out.printf("Integration result: %.6f ± %.2e\n",
754
                    result.get(), abserr.get());
755
                System.out.printf("Analytical result: %.6f\n", 8.0/3.0);
756
            }
757
            
758
            gsl_integration_workspace.gsl_integration_workspace_free(workspace);
759
        }
760
    }
761
}
762
```