\ GNU Scientific Library interface Mon Sep 12 14:40:15 MDT 2005 \ Copyright (C) 2007, Sergey Plis \ \ This program is free software; you can redistribute it and/or modify \ it under the terms of the GNU General Public License as published by \ the Free Software Foundation; either version 2 of the License, or \ (at your option) any later version. \ \ This program is distributed in the hope that it will be useful, \ but WITHOUT ANY WARRANTY; without even the implied warranty of \ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the \ GNU General Public License for more details. \needs float import float warning off \needs locals| include locals.fs warning on \needs atlas include atlas.fs \needs callback include callback.fs \needs vectors include vectors.fs \needs complex include complex.fb Module GSL \ stole the hash function from hash.fs of bigforth | &14 Value Hashbits | 1 Hashbits << Value Hashlen Label (hash ( SI:string -- AX:key ) :R DX push .b lods $1F # AX and AX CX mov DX DX xor CX 1 # shr b IF SI ) AH mov SI inc THEN CX dec 0>= IF BEGIN .w SI ) DX mov 2 # SI add CX dec DX AX *2 I) AX lea 0< UNTIL THEN & Hashbits A#) CX mov AX DX mov AX shr DX AX add & Hashlen A#) CX mov CX dec CX AX and DX pop ret | Code Hash ( string -- key ) R: SI push AX SI mov (hash rel) call SI pop Next end-code also dos also complex also float also atlas also vectors s" libptcblas.so" getlib 0<> [IF] library libblas libptcblas.so depends libatlas [ELSE] library libblas libcblas.so depends libatlas [THEN] library libgsl libgsl.so.0 depends libblas legacy off \ some functions libgsl gsl_log1p df (fp) gsl_log1p ( df -- df ) libgsl gsl_acosh df (fp) gsl_acosh ( df -- df ) \ error handling Wed Sep 21 23:04:06 MDT 2005 libgsl gsl_set_error_handler ptr (int) gsl_set_error_handler ( function -- function ) libgsl gsl_strerror int (ptr) gsl_strerror callback 4:0 (void) int int int int callback; : cstr-fstr ( addr -- addr len ) 0 begin 2dup + c@ 0 = not while 1+ repeat ; | : .bold-red ." " ; | : .red ." " ; | : .reset ." " ; | : cb-test cr \ .bold-red ." GSL ERROR: " cr \ .reset cr 10 spaces gsl_strerror cstr-fstr type cr drop \ ." at line: " . cr drop \ ." of file: " cstr-fstr type cr 10 spaces cstr-fstr type cr \ .red -1 abort" failed at" ; ' cb-test 4:0 c_plus \ 1 2 c_plus 2:1call . | variable old_handler c_plus gsl_set_error_handler old_handler ! 0 Constant GSL_SUCCESS \ Special functions \ Digamma Function libgsl gsl_sf_psi_int int (fp) gsl_sf_psi_int ( n -- f:\psi[n] ) libgsl gsl_sf_psi_int_e int ptr (int) gsl_sf_psi_int_e ( n *r -- n ) libgsl gsl_sf_psi df (fp) gsl_sf_psi libgsl gsl_sf_psi_e int ptr (int) gsl_sf_psi_e ( n *r -- n ) ' gsl_sf_psi_int alias ipsi ' gsl_sf_psi alias fpsi \ Gamma function libgsl gsl_sf_gamma df (fp) gsl_sf_gamma ' gsl_sf_gamma alias fgamma \ random number generation Mon Sep 12 22:06:01 MDT 2005 libgsl gsl_rng_types_setup (ptr) gsl_rng_types_setup ( -- *gsl_rng_type) libgsl gsl_rng_env_setup (ptr) gsl_rng_env_setup ( -- *gsl_rng) libgsl gsl_rng_alloc int (int) gsl_rng_alloc ( *gsl_rng_type -- *gsl_rng ) libgsl gsl_rng_name int (int) gsl_rng_name ( *gsl_rng -- string ) libgsl gsl_rng_set int int (void) gsl_rng_set ( *gsl_rng int -- ) libgsl gsl_rng_uniform int (fp) gsl_rng_uniform ( *gsl_rng -- df ) libgsl gsl_rng_uniform_pos int (fp) gsl_rng_uniform_pos ( *gsl_rng -- df ) libgsl gsl_rng_uniform_int int int (int) gsl_rng_uniform_int ( *gsl_rng n --n ) libgsl gsl_rng_get int (int) gsl_rng_get ( *gsl_rng -- int ) libgsl gsl_rng_max int (int) gsl_rng_max ( *gsl_rng -- int ) libgsl gsl_rng_min int (int) gsl_rng_min ( *gsl_rng -- int ) libgsl gsl_rng_clone int (int) gsl_rng_clone ( *gsl_rng -- *gsl_rng ) libgsl gsl_rng_free int (void) gsl_rng_free ( *gsl_rng -- ) libgsl _gsl_rng_borosh13 (int) gsl_rng_borosh13 libgsl _gsl_rng_coveyou (int) gsl_rng_coveyou libgsl _gsl_rng_cmrg (int) gsl_rng_cmrg libgsl _gsl_rng_fishman18 (int) gsl_rng_fishman18 libgsl _gsl_rng_fishman20 (int) gsl_rng_fishman20 libgsl _gsl_rng_fishman2x (int) gsl_rng_fishman2x libgsl _gsl_rng_gfsr4 (int) gsl_rng_gfsr4 \ 2nd fastest libgsl _gsl_rng_knuthran (int) gsl_rng_knuthran libgsl _gsl_rng_knuthran2 (int) gsl_rng_knuthran2 libgsl _gsl_rng_knuthran2002 (int) gsl_rng_knuthran2002 libgsl _gsl_rng_lecuyer21 (int) gsl_rng_lecuyer21 libgsl _gsl_rng_minstd (int) gsl_rng_minstd libgsl _gsl_rng_mrg (int) gsl_rng_mrg libgsl _gsl_rng_mt19937 (int) gsl_rng_mt19937 \ 3rd fastest libgsl _gsl_rng_mt19937_1999 (int) gsl_rng_mt19937_1999 libgsl _gsl_rng_mt19937_1998 (int) gsl_rng_mt19937_1998 libgsl _gsl_rng_r250 (int) gsl_rng_r250 libgsl _gsl_rng_ran0 (int) gsl_rng_ran0 libgsl _gsl_rng_ran1 (int) gsl_rng_ran1 libgsl _gsl_rng_ran2 (int) gsl_rng_ran2 libgsl _gsl_rng_ran3 (int) gsl_rng_ran3 libgsl _gsl_rng_rand (int) gsl_rng_rand libgsl _gsl_rng_rand48 (int) gsl_rng_rand48 libgsl _gsl_rng_random128_bsd (int) gsl_rng_random128_bsd libgsl _gsl_rng_random128_glibc2 (int) gsl_rng_random128_glibc2 libgsl _gsl_rng_random128_libc5 (int) gsl_rng_random128_libc5 libgsl _gsl_rng_random256_bsd (int) gsl_rng_random256_bsd libgsl _gsl_rng_random256_glibc2 (int) gsl_rng_random256_glibc2 libgsl _gsl_rng_random256_libc5 (int) gsl_rng_random256_libc5 libgsl _gsl_rng_random32_bsd (int) gsl_rng_random32_bsd libgsl _gsl_rng_random32_glibc2 (int) gsl_rng_random32_glibc2 libgsl _gsl_rng_random32_libc5 (int) gsl_rng_random32_libc5 libgsl _gsl_rng_random64_bsd (int) gsl_rng_random64_bsd libgsl _gsl_rng_random64_glibc2 (int) gsl_rng_random64_glibc2 libgsl _gsl_rng_random64_libc5 (int) gsl_rng_random64_libc5 libgsl _gsl_rng_random8_bsd (int) gsl_rng_random8_bsd libgsl _gsl_rng_random8_glibc2 (int) gsl_rng_random8_glibc2 libgsl _gsl_rng_random8_libc5 (int) gsl_rng_random8_libc5 libgsl _gsl_rng_random_bsd (int) gsl_rng_random_bsd libgsl _gsl_rng_random_glibc2 (int) gsl_rng_random_glibc2 libgsl _gsl_rng_random_libc5 (int) gsl_rng_random_libc5 libgsl _gsl_rng_randu (int) gsl_rng_randu libgsl _gsl_rng_ranf (int) gsl_rng_ranf libgsl _gsl_rng_ranlux (int) gsl_rng_ranlux \ statistically best libgsl _gsl_rng_ranlux389 (int) gsl_rng_ranlux389 \ libgsl _gsl_rng_ranlxd1 (int) gsl_rng_ranlxd1 \ libgsl _gsl_rng_ranlxd2 (int) gsl_rng_ranlxd2 libgsl _gsl_rng_ranlxs0 (int) gsl_rng_ranlxs0 libgsl _gsl_rng_ranlxs1 (int) gsl_rng_ranlxs1 libgsl _gsl_rng_ranlxs2 (int) gsl_rng_ranlxs2 \ up to here libgsl _gsl_rng_ranmar (int) gsl_rng_ranmar libgsl _gsl_rng_slatec (int) gsl_rng_slatec libgsl _gsl_rng_taus (int) gsl_rng_taus \ 1st fastest libgsl _gsl_rng_taus2 (int) gsl_rng_taus2 libgsl _gsl_rng_taus113 (int) gsl_rng_taus113 libgsl _gsl_rng_transputer (int) gsl_rng_transputer libgsl _gsl_rng_tt800 (int) gsl_rng_tt800 libgsl _gsl_rng_uni (int) gsl_rng_uni libgsl _gsl_rng_uni32 (int) gsl_rng_uni32 libgsl _gsl_rng_vax (int) gsl_rng_vax libgsl _gsl_rng_waterman14 (int) gsl_rng_waterman14 libgsl _gsl_rng_zuf (int) gsl_rng_zuf libgsl _gsl_rng_default (int) gsl_rng_default \ random number distributions Tue Sep 13 00:44:35 MDT 2005 \ Gaussian libgsl gsl_ran_gaussian df int (fp) gsl_ran_gaussian ( *gsl_rng df -- df ) libgsl gsl_ran_gaussian_ratio_method df int (fp) gsl_ran_gaussian_ratio_method ( *gsl_rng df -- df ) libgsl gsl_ran_gaussian_pdf df df (fp) gsl_ran_gaussian_pdf ( df df -- df ) \ sigma = 1 libgsl gsl_ran_ugaussian int (fp) gsl_ran_ugaussian ( *gsl_rng -- df ) libgsl gsl_ran_ugaussian_ratio_method int (fp) gsl_ran_ugaussian_ratio_method ( *gsl_rng -- df ) libgsl gsl_ran_ugaussian_pdf df (fp) gsl_ran_ugaussian_pdf ( df df -- df ) libgsl gsl_ran_discrete_preproc int int (int) gsl_ran_discrete_preproc ( int int -- int ) libgsl gsl_ran_discrete int int (int) gsl_ran_discrete libgsl gsl_ran_discrete_free int (void) gsl_ran_discrete_free libgsl gsl_ran_shuffle int int ptr ptr (void) gsl_ran_shuffle libgsl gsl_ran_choose int int ptr int ptr ptr (int) gsl_ran_choose libgsl gsl_ran_sample int int ptr int ptr ptr (void) gsl_ran_sample \ cdf P(x) = \int_{-\infty}^{x} p(x)dx Q(x) = \int_{x}^{\infty} p(x)dx libgsl gsl_cdf_gaussian_P df df (fp) gsl_cdf_gaussian_P ( df df -- df ) libgsl gsl_cdf_gaussian_Q df df (fp) gsl_cdf_gaussian_Q ( df df -- df ) libgsl gsl_cdf_gaussian_Pinv df df (fp) gsl_cdf_gaussian_Pinv ( df df -- df ) libgsl gsl_cdf_gaussian_Qinv df df (fp) gsl_cdf_gaussian_Qinv ( df df -- df ) \ sigma = 1 cdf libgsl gsl_cdf_ugaussian_P df (fp) gsl_cdf_ugaussian_P ( df -- df ) libgsl gsl_cdf_ugaussian_Q df (fp) gsl_cdf_ugaussian_Q ( df -- df ) libgsl gsl_cdf_ugaussian_Pinv df (fp) gsl_cdf_ugaussian_Pinv ( df -- df ) libgsl gsl_cdf_ugaussian_Qinv df (fp) gsl_cdf_ugaussian_Qinv ( df -- df ) \ statistics Tue Sep 13 01:17:35 MDT 2005 libgsl gsl_stats_mean int int int (fp) gsl_stats_mean ( array{ step size -- df ) libgsl gsl_stats_variance int int int (fp) gsl_stats_variance ( array{ step size -- df ) libgsl gsl_stats_variance_m df int int int (fp) gsl_stats_variance_m ( df array{ step size -- df ) libgsl gsl_stats_sd int int int (fp) gsl_stats_sd ( array{ step size -- df ) libgsl gsl_stats_sd_m df int int int (fp) gsl_stats_sd_m ( df array{ step size -- df ) libgsl gsl_stats_skew int int int (fp) gsl_stats_skew ( array{ step size -- df ) libgsl gsl_stats_kurtosis int int int (fp) gsl_stats_kurtosis ( array{ step size -- df ) libgsl gsl_stats_lag1_autocorrelation int int int (fp) gsl_stats_lag1_autocorrelation ( array{ step size -- df ) libgsl gsl_stats_max int int int (fp) gsl_stats_max ( array{ step size -- df ) libgsl gsl_stats_min int int int (fp) gsl_stats_min ( array{ step size -- df ) libgsl gsl_stats_max_index int int int (int) gsl_stats_max_index ( array{ step size -- n ) libgsl gsl_stats_min_index int int int (int) gsl_stats_min_index ( array{ step size -- n ) \ vectors and matrices Wed Sep 14 00:15:36 MDT 2005 \ Vectors libgsl gsl_block_alloc int (int) gsl_block_alloc ( n -- addr ) libgsl gsl_block_calloc int (int) gsl_block_calloc ( n -- addr ) libgsl gsl_block_free int (int) gsl_block_free ( n -- addr ) libgsl gsl_vector_alloc int (int) gsl_vector_alloc ( n -- addr ) libgsl gsl_vector_calloc int (int) gsl_vector_calloc ( n -- addr ) libgsl gsl_vector_alloc_from_vector int int int ptr (int) gsl_vector_alloc_from_vector libgsl gsl_vector_free int (void) gsl_vector_free ( addr -- ) libgsl gsl_vector_get int int (fp) gsl_vector_get ( addr i -- df ) libgsl gsl_vector_set df int int (void/fp) gsl_vector_set ( df addr i -- ) libgsl gsl_vector_set_all df int (void) gsl_vector_set_all ( df addr -- ) libgsl gsl_vector_set_zero int (void) gsl_vector_set_zero ( addr -- ) libgsl gsl_vector_memcpy int int (int) gsl_vector_memcpy ( dest_addr src_addr -- n ) libgsl gsl_vector_add int int (int) gsl_vector_add ( addr addr -- n ) libgsl gsl_vector_sub int int (int) gsl_vector_sub ( addr addr -- n ) libgsl gsl_vector_mul int int (int) gsl_vector_mul ( addr addr -- n ) libgsl gsl_vector_div int int (int) gsl_vector_div ( addr addr -- n ) libgsl gsl_vector_scale df int (int) gsl_vector_scale ( df addr -- n ) libgsl gsl_vector_add_constant df int (int) gsl_vector_add_constant ( df addr -- n ) libgsl gsl_vector_max int (fp) gsl_vector_max ( addr -- df ) libgsl gsl_vector_min int (fp) gsl_vector_min ( addr -- df ) libgsl gsl_vector_minmax ptr ptr ptr (void) gsl_vector_minmax libgsl gsl_vector_max_index int (int) gsl_vector_max_index ( a -- n ) libgsl gsl_vector_min_index int (int) gsl_vector_min_index ( a -- n ) libgsl gsl_vector_minmax_index ptr ptr ptr (void) gsl_vector_minmax_index libgsl gsl_vector_subvector int int int (int) gsl_vector_subvector \ Vector properties libgsl gsl_vector_isnull ptr (int) gsl_vector_isnull libgsl gsl_vector_ispos ptr (int) gsl_vector_ispos libgsl gsl_vector_isneg ptr (int) gsl_vector_isneg \ Sorting vectors libgsl gsl_sort int int ptr (void) gsl_sort libgsl gsl_sort_vector ptr (void) gsl_sort_vector libgsl gsl_sort_index int int ptr ptr (void) gsl_sort_index libgsl gsl_sort_vector_index ptr ptr (int) gsl_sort_vector_index libgsl gsl_sort_vector_smallest_index ptr int ptr (int) gsl_sort_vector_smallest_index \ permutations libgsl gsl_permutation_alloc int (int) gsl_permutation_alloc ( n -- *gsl_prm) libgsl gsl_permutation_calloc int (int) gsl_permutation_calloc ( n -- *gsl_prm) libgsl gsl_permutation_init int (void) gsl_permutation_init ( *gsl_prm -- ) libgsl gsl_permutation_free int (void) gsl_permutation_free ( *gsl_prm -- ) libgsl gsl_permutation_get int int (int) gsl_permutation_get ( *gsl_prm i -- n) libgsl gsl_permutation_reverse ptr (void) gsl_permutation_reverse libgsl gsl_permutation_inverse ptr ptr (int) gsl_permutation_inverse \ Matrices \ Allocation libgsl gsl_matrix_alloc int int (int) gsl_matrix_alloc libgsl gsl_matrix_calloc int int (int) gsl_matrix_calloc libgsl gsl_matrix_alloc_from_block [ 5 ] ints (int) gsl_matrix_alloc_from_block libgsl gsl_matrix_alloc_from_matrix [ 5 ] ints (int) gsl_matrix_alloc_from_matrix libgsl gsl_matrix_free ( *gsl_matrix -- ) int (void) gsl_matrix_free \ Accessing matrix elements libgsl gsl_matrix_get int int int (fp) gsl_matrix_get ( *m i j -- df ) libgsl gsl_matrix_set df int int int (void) gsl_matrix_set ( df *m i j -- ) libgsl gsl_matrix_ptr int int int (int) gsl_matrix_ptr ( *m i j -- *[i,j] ) \ Initializing matrix elements libgsl gsl_matrix_set_all df int (void) gsl_matrix_set_all ( *m df -- n ) libgsl gsl_matrix_set_zero df int (void) gsl_matrix_set_zero ( *m df -- n ) libgsl gsl_matrix_set_identity df int (void) gsl_matrix_set_identity ( *m df -- n ) \ Reading and writing matrices libgsl gsl_matrix_fwrite ptr ptr (int) gsl_matrix_fwrite libgsl gsl_matrix_fread ptr ptr (int) gsl_matrix_fread libgsl gsl_matrix_fprintf ptr ptr ptr (int) gsl_matrix_fprintf libgsl gsl_matrix_fscanf ptr ptr (int) gsl_matrix_fscanf \ Copying matrices libgsl gsl_matrix_memcpy int int (int) gsl_matrix_memcpy ( *m *m -- n ) libgsl gsl_matrix_swap int int (int) gsl_matrix_swap ( *m *m -- n ) \ Copying Rows and columns libgsl gsl_matrix_get_row int int int (int) gsl_matrix_get_row libgsl gsl_matrix_set_row int int int (int) gsl_matrix_set_row libgsl gsl_matrix_get_col int int int (int) gsl_matrix_get_col libgsl gsl_matrix_set_col int int int (int) gsl_matrix_set_col \ Exchanging rows and columns libgsl gsl_matrix_swap_rows int int ptr (int) gsl_matrix_swap_rows libgsl gsl_matrix_swap_columns int int ptr (int) gsl_matrix_swap_columns libgsl gsl_matrix_swap_rowcol int int ptr (int) gsl_matrix_swap_rowcol libgsl gsl_matrix_transpose_memcpy int int (int) gsl_matrix_transpose_memcpy libgsl gsl_matrix_transpose int (int) gsl_matrix_transpose \ Matrix operations libgsl gsl_matrix_add int int (int) gsl_matrix_add libgsl gsl_matrix_sub int int (int) gsl_matrix_sub libgsl gsl_matrix_mul_elements int int (int) gsl_matrix_mul_elements libgsl gsl_matrix_div_elements int int (int) gsl_matrix_div_elements libgsl gsl_matrix_scale df int (int) gsl_matrix_scale libgsl gsl_matrix_add_constant df int (int) gsl_matrix_add_constant \ Finding maximum and minimum elements of matrices libgsl gsl_matrix_max ptr (fp) gsl_matrix_max libgsl gsl_matrix_min ptr (fp) gsl_matrix_min libgsl gsl_matrix_minmax ptr ptr ptr (void) gsl_matrix_minmax libgsl gsl_matrix_min_index ptr ptr ptr (void) gsl_matrix_min_index libgsl gsl_matrix_max_index ptr ptr ptr (void) gsl_matrix_max_index libgsl gsl_matrix_minmax_index ptr ptr ptr ptr ptr (void) gsl_matrix_minmax_index \ Matrix properties libgsl gsl_matrix_isnull ptr (int) gsl_matrix_isnull libgsl gsl_matrix_ispos ptr (int) gsl_matrix_ispos libgsl gsl_matrix_isneg ptr (int) gsl_matrix_isneg \ libgsl gsl_matrix_isnonneg ptr (int) gsl_matrix_isnonneg libgsl gsl_matrix_submatrix int int int int int (int) gsl_matrix_submatrix ( *gsl_matrix k1 k2 n1 n2 -- n ) libgsl gsl_matrix_row int int (int) gsl_matrix_row ( *gsl_matrix idx -- *gsl_vector ) libgsl gsl_matrix_column int int (int) gsl_matrix_column ( *gsl_matrix idx -- *gsl_vector ) libgsl gsl_matrix_diagonal int (int) gsl_matrix_diagonal ( *gsl_matrix -- *gsl_vector ) \ BLAS Wed Sep 14 16:10:34 MDT 2005 \ libblas cblas_dgemm int int df int int int \ int df int int int int int int (void/fp) cblas_dgemm libblas cblas_dgemv int int int int df int int df int int int int (void/fp) cblas_dgemv libgsl gsl_blas_ddot int int int (int) gsl_blas_ddot ( *gsl_vector *gsl_vector df -- n ) libgsl gsl_blas_dgemm int df int int df int int (int/fp) gsl_blas_dgemm libgsl gsl_blas_dger int int int df (int/fp) gsl_blas_dger ( alpha *gsl_vector *gsl_vector *gsl_matrix -- n ) ( A=\alpha x y^T+A ) libgsl gsl_blas_dgemv int df int int df int (int/fp) gsl_blas_dgemv ( n alpha *gsl_matrix *gsl_vector beta *gsl_vector -- n ) \ Linear ALgebra Wed Sep 14 13:39:22 MDT 2005 libgsl gsl_linalg_LU_decomp int int int (int) gsl_linalg_LU_decomp ( *gsl_matrix *gsl_permutation *variable -- n ) libgsl gsl_linalg_LU_invert int int int (int) gsl_linalg_LU_invert ( *gsl_matrix *gsl_permutation *gsl_matrix -- n ) libgsl gsl_linalg_SV_decomp int int int int (int) gsl_linalg_SV_decomp ( *gsl_matrix *gsl_matrix *gsl_vector *gsl_vector -- n ) libgsl gsl_linalg_SV_decomp_mod int int int int int (int) gsl_linalg_SV_decomp_mod ( *gsl_matrix *gsl_matrix *gsl_matrix *gsl_vector *gsl_vector -- n ) \ ----------------------------------------------------------------------------- \ *** Ordinary Differential Equations *** \ --- ODE system struct{ cell func \ (* function) \ (double t, const double y[], double dydt[], void * params); cell jac \ (* jacobian) \ (double t, const double y[], double * dfdy, double dfdt[], \ void * params); cell dim \ dimension; cell params \ * params; } gsl_odeiv_system \ constants related to ODE 1 constant GSL_ODEIV_HADJ_INC 0 constant GSL_ODEIV_HADJ_NIL -1 constant GSL_ODEIV_HADJ_DEC callback gsl_odeiv_func4:1 (int) df int int int callback; callback gsl_odeiv_jac5:1 (int) df int int int int callback; \ --- Stepping Functions libgsl gsl_odeiv_step_alloc ptr int (ptr) gsl_odeiv_step_alloc ( *step_type int -- *step ) libgsl gsl_odeiv_step_reset ptr (int) gsl_odeiv_step_reset ( *step -- r ) libgsl gsl_odeiv_step_free ptr (void) gsl_odeiv_step_free ( *step -- ) libgsl gsl_odeiv_step_name ptr (ptr) gsl_odeiv_step_name ( *step -- *str0 ) libgsl gsl_odeiv_step_order ptr (int) gsl_odeiv_step_order ( *step -- order) libgsl gsl_odeiv_step_apply int int int int int df df int (int) gsl_odeiv_step_apply ( -- ) \ --- Available algorithms libgsl _gsl_odeiv_step_rk2 (int) gsl_odeiv_step_rk2 libgsl _gsl_odeiv_step_rk4 (int) gsl_odeiv_step_rk4 libgsl _gsl_odeiv_step_rkf45 (int) gsl_odeiv_step_rkf45 libgsl _gsl_odeiv_step_rkck (int) gsl_odeiv_step_rkck libgsl _gsl_odeiv_step_rk8pd (int) gsl_odeiv_step_rk8pd libgsl _gsl_odeiv_step_rk2imp (int) gsl_odeiv_step_rk2imp libgsl _gsl_odeiv_step_rk4imp (int) gsl_odeiv_step_rk4imp libgsl _gsl_odeiv_step_bsimp (int) gsl_odeiv_step_bsimp libgsl _gsl_odeiv_step_gear1 (int) gsl_odeiv_step_gear1 libgsl _gsl_odeiv_step_gear2 (int) gsl_odeiv_step_gear2 : gsl_odeiv_step_rk2 [func'] _gsl_odeiv_step_rk2 @ ; : gsl_odeiv_step_rk4 [func'] _gsl_odeiv_step_rk4 @ ; : gsl_odeiv_step_rkf45 [func'] _gsl_odeiv_step_rkf45 @ ; : gsl_odeiv_step_rkck [func'] _gsl_odeiv_step_rkck @ ; : gsl_odeiv_step_rk8pd [func'] _gsl_odeiv_step_rk8pd @ ; : gsl_odeiv_step_rk2imp [func'] _gsl_odeiv_step_rk2imp @ ; : gsl_odeiv_step_rk4imp [func'] _gsl_odeiv_step_rk4imp @ ; : gsl_odeiv_step_bsimp [func'] _gsl_odeiv_step_bsimp @ ; : gsl_odeiv_step_gear1 [func'] _gsl_odeiv_step_gear1 @ ; : gsl_odeiv_step_gear2 [func'] _gsl_odeiv_step_gear2 @ ; \ --- Adaptive Step-size Control libgsl gsl_odeiv_control_standard_new df df df df (ptr) gsl_odeiv_control_standard_new ( a_dydt a_y eps_rel eps_abs -- *control ) libgsl gsl_odeiv_control_y_new df df (int) gsl_odeiv_control_y_new ( eps_abs eps_rel -- *control ) libgsl gsl_odeiv_control_yp_new df df (ptr) gsl_odeiv_control_yp_new ( eps_abs eps_rel -- *control ) libgsl gsl_odeiv_control_free ptr (void) gsl_odeiv_control_free ( *control -- ) libgsl gsl_odeiv_control_name ptr (ptr) gsl_odeiv_control_name ( *c -- *str0 ) \ --- Evolution libgsl gsl_odeiv_evolve_alloc int (int) gsl_odeiv_evolve_alloc ( #dimensions -- evolution_func ) libgsl gsl_odeiv_evolve_apply int int df int int int int int (int) gsl_odeiv_evolve_apply ( -- ) libgsl gsl_odeiv_evolve_reset ptr (int) gsl_odeiv_evolve_reset ( *e -- r ) libgsl gsl_odeiv_evolve_free ptr (void) gsl_odeiv_evolve_free ( *e -- ) \ ----------------------------------------------------------------------------- \ *** Fast Fourier Transform *** \ -- real libgsl gsl_fft_real_wavetable_alloc int (ptr) gsl_fft_real_wavetable_alloc libgsl gsl_fft_real_wavetable_free ptr (void) gsl_fft_real_wavetable_free libgsl gsl_fft_real_workspace_alloc int (ptr) gsl_fft_real_workspace_alloc libgsl gsl_fft_real_workspace_free ptr (void) gsl_fft_real_workspace_free \ in-place libgsl gsl_fft_real_transform ptr int int ptr ptr (int) gsl_fft_real_transform libgsl gsl_fft_real_unpack ptr ptr int int (int) gsl_fft_real_unpack \ -- halfcomplex \ - mixed radix libgsl gsl_fft_hc_wtbl_alloc int (ptr) gsl_fft_halfcomplex_wavetable_alloc libgsl gsl_fft_hc_wtbl_free ptr (void) gsl_fft_halfcomplex_wavetable_free libgsl gsl_fft_hc_backward ptr int int ptr ptr (int) gsl_fft_halfcomplex_backward libgsl gsl_fft_hc_inverse ptr int int ptr ptr (int) gsl_fft_halfcomplex_inverse libgsl gsl_fft_hc_transform ptr int int ptr ptr (int) gsl_fft_halfcomplex_transform libgsl gsl_fft_hc_unpack ptr ptr int int (int) gsl_fft_halfcomplex_unpack \ - radix2 libgsl gsl_fft_hc_r2_unpack ptr ptr int int (int) gsl_fft_halfcomplex_radix2_unpack libgsl gsl_fft_hc_r2_backward ptr int int (int) gsl_fft_halfcomplex_radix2_backward libgsl gsl_fft_hc_r2_inverse ptr int int (int) gsl_fft_halfcomplex_radix2_inverse libgsl gsl_fft_hc_r2_transform ptr int int (int) gsl_fft_halfcomplex_radix2_transform | hashlen 32 vector fftpre( struct{ cell next cell size cell workspace cell r_wavetable cell hc_wavetable } gsl_fft_precomputes | create $buf 255 allot | : 2str dup >r abs s>d <<# #s r> sign #> $buf 0place #>> ; | : s>hash ( n -- key ) 2str $buf hash ; | : (cache-fft) ( n -- addr ) sizeof gsl_fft_precomputes allocate throw >r 0 r@ gsl_fft_precomputes next ! dup r@ gsl_fft_precomputes size ! dup gsl_fft_real_workspace_alloc r@ gsl_fft_precomputes workspace ! dup gsl_fft_real_wavetable_alloc r@ gsl_fft_precomputes r_wavetable ! gsl_fft_hc_wtbl_alloc r@ gsl_fft_precomputes hc_wavetable ! r> ; | : cache-fft ( size -- addr ) dup s>hash fftpre( over )@ 0= if swap (cache-fft) fftpre( rot dup >r )! fftpre( r> )@ else swap (cache-fft) swap fftpre( over )@ over gsl_fft_precomputes next ! fftpre( rot dup >r )! fftpre( r> )@ then ; \ in case not found addr is just the size | : find-fft-cache ( n -- addr 0/1 ) dup s>hash fftpre( swap )@ dup begin while 2dup gsl_fft_precomputes size @ = if nip true exit then gsl_fft_precomputes next @ dup repeat ; legacy on \ Structures struct{ cell name cell max cell min cell size cell set cell get cell get_double } gsl_rng_type struct{ cell type cell state } gsl_rng struct{ cell size cell data } gsl_block struct{ cell size cell stride cell data cell block cell owner } gsl_vector ' gsl_block alias gsl_permutation struct{ cell size1 cell size2 cell tda cell data cell block cell owner } gsl_matrix \ random number generation functions : 0-len dup 1- 0 begin 1+ 2dup + c@ 0= until nip ; : )gsl-rng ( addr i -- *gsl_rng_type ) cells + @ ; \ setting up all available random number generators \ gsl_rng_types_setup value gsl_rng_array( 0 value gsl_rng_default : gsl-free ( -- ) gsl_rng_default gsl_rng_free ; \ the following code was generated automatically, for description \ consult GSL documentation : borosh13 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_borosh13 @ gsl_rng_alloc to gsl_rng_default ; : coveyou ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_coveyou @ gsl_rng_alloc to gsl_rng_default ; : cmrg ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_cmrg @ gsl_rng_alloc to gsl_rng_default ; : fishman18 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_fishman18 @ gsl_rng_alloc to gsl_rng_default ; : fishman20 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_fishman20 @ gsl_rng_alloc to gsl_rng_default ; : fishman2x ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_fishman2x @ gsl_rng_alloc to gsl_rng_default ; : gfsr4 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_gfsr4 @ gsl_rng_alloc to gsl_rng_default ; : knuthran ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_knuthran @ gsl_rng_alloc to gsl_rng_default ; : knuthran2 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_knuthran2 @ gsl_rng_alloc to gsl_rng_default ; : knuthran2002 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_knuthran2002 @ gsl_rng_alloc to gsl_rng_default ; : lecuyer21 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_lecuyer21 @ gsl_rng_alloc to gsl_rng_default ; : minstd ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_minstd @ gsl_rng_alloc to gsl_rng_default ; : mrg ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_mrg @ gsl_rng_alloc to gsl_rng_default ; : mt19937 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_mt19937 @ gsl_rng_alloc to gsl_rng_default ; : mt19937_1999 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_mt19937_1999 @ gsl_rng_alloc to gsl_rng_default ; : mt19937_1998 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_mt19937_1998 @ gsl_rng_alloc to gsl_rng_default ; : r250 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_r250 @ gsl_rng_alloc to gsl_rng_default ; : ran0 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ran0 @ gsl_rng_alloc to gsl_rng_default ; : ran1 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ran1 @ gsl_rng_alloc to gsl_rng_default ; : ran2 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ran2 @ gsl_rng_alloc to gsl_rng_default ; : ran3 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ran3 @ gsl_rng_alloc to gsl_rng_default ; : rand ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_rand @ gsl_rng_alloc to gsl_rng_default ; : rand48 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_rand48 @ gsl_rng_alloc to gsl_rng_default ; : random128_bsd ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random128_bsd @ gsl_rng_alloc to gsl_rng_default ; : random128_glibc2 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random128_glibc2 @ gsl_rng_alloc to gsl_rng_default ; : random128_libc5 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random128_libc5 @ gsl_rng_alloc to gsl_rng_default ; : random256_bsd ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random256_bsd @ gsl_rng_alloc to gsl_rng_default ; : random256_glibc2 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random256_glibc2 @ gsl_rng_alloc to gsl_rng_default ; : random256_libc5 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random256_libc5 @ gsl_rng_alloc to gsl_rng_default ; : random32_bsd ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random32_bsd @ gsl_rng_alloc to gsl_rng_default ; : random32_glibc2 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random32_glibc2 @ gsl_rng_alloc to gsl_rng_default ; : random32_libc5 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random32_libc5 @ gsl_rng_alloc to gsl_rng_default ; : random64_bsd ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random64_bsd @ gsl_rng_alloc to gsl_rng_default ; : random64_glibc2 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random64_glibc2 @ gsl_rng_alloc to gsl_rng_default ; : random64_libc5 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random64_libc5 @ gsl_rng_alloc to gsl_rng_default ; : random8_bsd ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random8_bsd @ gsl_rng_alloc to gsl_rng_default ; : random8_glibc2 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random8_glibc2 @ gsl_rng_alloc to gsl_rng_default ; : random8_libc5 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random8_libc5 @ gsl_rng_alloc to gsl_rng_default ; : random_bsd ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random_bsd @ gsl_rng_alloc to gsl_rng_default ; : random_glibc2 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random_glibc2 @ gsl_rng_alloc to gsl_rng_default ; : random_libc5 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_random_libc5 @ gsl_rng_alloc to gsl_rng_default ; : randu ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_randu @ gsl_rng_alloc to gsl_rng_default ; : ranf ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ranf @ gsl_rng_alloc to gsl_rng_default ; : ranlux ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ranlux @ gsl_rng_alloc to gsl_rng_default ; : ranlux389 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ranlux389 @ gsl_rng_alloc to gsl_rng_default ; : ranlxd1 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ranlxd1 @ gsl_rng_alloc to gsl_rng_default ; : ranlxd2 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ranlxd2 @ gsl_rng_alloc to gsl_rng_default ; : ranlxs0 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ranlxs0 @ gsl_rng_alloc to gsl_rng_default ; : ranlxs1 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ranlxs1 @ gsl_rng_alloc to gsl_rng_default ; : ranlxs2 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ranlxs2 @ gsl_rng_alloc to gsl_rng_default ; : ranmar ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_ranmar @ gsl_rng_alloc to gsl_rng_default ; : slatec ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_slatec @ gsl_rng_alloc to gsl_rng_default ; : taus ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_taus @ gsl_rng_alloc to gsl_rng_default ; : taus2 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_taus2 @ gsl_rng_alloc to gsl_rng_default ; : taus113 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_taus113 @ gsl_rng_alloc to gsl_rng_default ; : transputer ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_transputer @ gsl_rng_alloc to gsl_rng_default ; : tt800 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_tt800 @ gsl_rng_alloc to gsl_rng_default ; : uni ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_uni @ gsl_rng_alloc to gsl_rng_default ; : uni32 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_uni32 @ gsl_rng_alloc to gsl_rng_default ; : vax ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_vax @ gsl_rng_alloc to gsl_rng_default ; : waterman14 ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_waterman14 @ gsl_rng_alloc to gsl_rng_default ; : zuf ( -- ) gsl_rng_default 0<> if gsl-free then [func'] _gsl_rng_zuf @ gsl_rng_alloc to gsl_rng_default ; \ : default ( -- ) gsl_rng_default 0<> if gsl-free then \ [func'] _gsl_rng_default @ gsl_rng_alloc to gsl_rng_default ; \ words for actual generation of random numbers : gsl-randomg ( -- n ) gsl_rng_default gsl_rng_get ; : gsl-randomu /* -- f \in [0,1) */ gsl_rng_default gsl_rng_uniform ; : gsl-randomu+ /* -- f \in (0,1) */ gsl_rng_default gsl_rng_uniform_pos ; : gsl-random-up ( n -- f \in [0,n] ) gsl_rng_default swap gsl_rng_uniform_int ; : gsl-set-seed ( n -- ) gsl_rng_default swap gsl_rng_set ; : gsl-clone ( -- *gsl_rng ) gsl_rng_default gsl_rng_clone ; : gsl-gaussian ( -- df ) gsl_rng_default !1 gsl_ran_gaussian_ratio_method ; : gsl-discrete ( *gsl_ran_discrete -- n ) gsl_rng_default swap gsl_ran_discrete ; \ vectors and matrices 0 constant GSL_VECTOR_TYPE 1 constant GSL_MATRIX_TYPE : gsltypeof ' >body cell + @ ; : fvector ( n -- | -- id addr ) create gsl_vector_calloc , GSL_VECTOR_TYPE , does> @ ; \ allocate a nameless vector : :] ( # -- addr ) gsl_vector_calloc ; \ allocate a nameless matrix : :]] ( # # -- addr ) gsl_matrix_calloc ; : ]@ ( addr i -- df ) gsl_vector_get ; : ]! ( df addr i -- ) gsl_vector_set ; : ]data ( addr -- *data ) gsl_vector data @ ; : ]stride ( addr -- *data ) gsl_vector stride @ ; : ]fill ( df addr -- ) gsl_vector_set_all ; : ]erase ( addr -- ) gsl_vector_set_zero ; : ]+ ( *gsl_vector *gsl_vector -- ) gsl_vector_add drop ; : ]- ( *gsl_vector *gsl_vector -- ) gsl_vector_sub drop ; : ]e*! ( *gsl_vector *gsl_vector -- ) gsl_vector_mul drop ; : ]size ( *gsl_vector -- n ) gsl_vector size @ ; : ]outer* ( *gsl_vector *gsl_vector -- *gsl_matrix ) over ]size over ]size gsl_matrix_calloc dup >r !1 gsl_blas_dger drop r> ; \ no control for division by zero in gsl (is it NAN now?) : ]/ ( *gsl_vector *gsl_vector -- ) gsl_vector_div throw ; : ]clone ( *gsl_vector -- *gsl_vector ) dup gsl_vector size @ gsl_vector_alloc dup -rot swap gsl_vector_memcpy drop ; ' gsl_sort_vector alias ]sort ( v[ -- ) : ]add ( *gsl_vector *gsl_vector -- *gsl_vector ) ]clone dup -rot swap ]+ ; : ]sub ( *gsl_vector *gsl_vector -- *gsl_vector ) swap ]clone dup -rot swap ]- ; : ]mul ( *gsl_vector *gsl_vector -- *gsl_vector ) swap ]clone dup -rot swap ]e*! ; : ]div ( *gsl_vector *gsl_vector -- *gsl_vector ) swap ]clone dup -rot swap ]/ ; \ "*, /, +, -" constant to a vector : ]*c ( df *gsl_vector -- ) gsl_vector_scale drop ; : ]/c ( df *gsl_vector -- ) 1/f ]*c ; : ]+c ( df *gsl_vector -- ) gsl_vector_add_constant drop ; : ]-c ( df *gsl_vector -- ) fnegate ]+c ; \ extrema finding words : ]max ( *gsl_vector -- f:max ) gsl_vector_max ; : ]min ( *gsl_vector -- f:min ) gsl_vector_min ; : ]minmax ( *gsl_vector -- f:min f:max ) >r 0. sp@ >r 0. sp@ >r r> r> r> -rot gsl_vector_minmax fd>f fd>f ; : ]imax ( *gsl_vector -- max_idx ) gsl_vector_max_index ; : ]imin ( *gsl_vector -- min_idx ) gsl_vector_min_index ; : ]iminmax ( *gsl_vector -- min max ) sp@ >r 0 sp@ >r 0 rot 2r> gsl_vector_minmax_index ; : ]copy] ( *gsl_vector_dest *gsl_vector_src -- ) gsl_vector_memcpy drop ; : ]negate !-1.0 ]*c ; : ]ones ( n -- *gsl_vector ) :] dup 1e ]fill ; : ]]slice ( *gsl_matrix x y n m -- *gsl_matrix ) gsl_matrix_alloc_from_matrix ; : ]slice ( *gsl_vector offset length stride -- *gsl_vector ) gsl_vector_alloc_from_vector ; : ]map ( *v -- ) ( word ) ' swap dup ]size 0 do 2dup i ]@ execute dup i ]! loop 2drop ; : (]map) ( *v xt -- ) swap dup ]size 0 do 2dup i ]@ execute dup i ]! loop 2drop ; : ]abs ( *v -- ) ['] fabs (]map) ; : ]indgen ( *v -- ) dup ]size 0 do i s>f dup i ]! loop drop ; : ]null? ( *gsl_vector -- 0/-1 ) gsl_vector_isnull negate ; : ]pos? ( *gsl_vector -- 0/-1 ) gsl_vector_ispos negate ; : ]neg? ( *gsl_vector -- 0/-1 ) gsl_vector_isneg negate ; \ FFT 19jan08sp : ]fft! ( *gsl_vector -- ) dup ]size >r dup ]stride >r ]data r> r> dup find-fft-cache if dup gsl_fft_precomputes r_wavetable @ swap gsl_fft_precomputes workspace @ else drop dup cache-fft dup gsl_fft_precomputes r_wavetable @ swap gsl_fft_precomputes workspace @ then gsl_fft_real_transform throw ; : ]fft ( *gsl_vector -- *gsl_vector ) ]clone dup ]fft! ; : ]ifft! ( *gsl_vector -- ) dup ]size >r dup ]stride >r ]data r> r> dup find-fft-cache if dup gsl_fft_precomputes hc_wavetable @ swap gsl_fft_precomputes workspace @ else drop dup cache-fft dup gsl_fft_precomputes hc_wavetable @ swap gsl_fft_precomputes workspace @ then gsl_fft_hc_inverse throw ; : ]ifft ( *gsl_vector -- *gsl_vector ) ]clone dup ]ifft! ; \ multiply two half complex vectors \ store result in the first - needed for convolution : ]hc*! ( *gsl_vector *gsl_vector -- ) 2dup 0 ]@ 0 ]@ f* over 0 ]! dup ]size dup %1 and not + 1 do 2dup dup i ]@ i 1+ ]@ dup i ]@ i 1+ ]@ z* over dup i 1+ ]! i ]! 2 +loop dup ]size %1 and not if 2dup dup ]size 1- dup >r ]@ dup r@ ]@ f* r> ]! then 2drop ; \ pseudomatrices and vectors : pvector ( *data n -- *gsl_vector ) sizeof gsl_vector allocate throw dup >r dup 1 swap gsl_vector stride ! gsl_vector size ! r@ gsl_vector data ! r@ 0 swap gsl_vector owner ! r> ; : pmatrix! ( *data tda n m *pmatrix -- *gsl_matrix ) dup >r gsl_matrix size2 ! r@ gsl_matrix size1 ! r@ gsl_matrix tda ! r@ gsl_matrix data ! 0 r@ gsl_matrix owner ! r> ; : pmatrix ( *data tda n m -- *gsl_matrix ) sizeof gsl_matrix allocate throw dup 0 swap gsl_matrix owner ! pmatrix! ; \ permutations : fpermutation ( n -- | -- id addr ) create gsl_permutation_calloc , does> @ ; : :} ( n -- *gsl_permutation ) gsl_permutation_calloc ; : }@ ( *gsl_permutation i -- n ) gsl_permutation_get ; : }data ( *gsl_permutation -- *data ) gsl_block data @ ; : }size ( *gsl_permutation -- *data ) gsl_block size @ ; : }free ( *gsl_permutation -- ) gsl_permutation_free ; : }sign ( *gsl_permutation -- 1/-1 ) 1 over dup }size 0 do dup i }@ i <> if swap negate swap then loop drop ; : }>[ ( p{ -- p[ ) dup }size :] dup ]size 0 do 2dup swap i }@ s>f i ]! loop nip ; : ]isort ( *gsl_vector -- *gsl_permutation ) dup ]size :} dup >r swap gsl_sort_vector_index throw r> ; : ]kisort ( *gsl_vector k -- *vector ) dup vector* dup >r -rot swap gsl_sort_vector_smallest_index throw r> ; \ matrices : fmatrix ( n m -- | -- id addr ) create gsl_matrix_calloc , GSL_MATRIX_TYPE , does> @ ; : free_pseudomatrix ( pmatrix/pvector -- ) free throw ; create free_matrix ' free_pseudomatrix , ' gsl_matrix_free , create free_vector ' free_pseudomatrix , ' gsl_vector_free , : ]]free ( *gsl_matrix -- ) dup gsl_matrix owner @ cells free_matrix + @ execute ; : ]free ( addr -- ) dup gsl_vector owner @ cells free_vector + @ execute ; : ]]@ ( *gsl_matrix i j -- df ) gsl_matrix_get ; : ]]*@ ( *gsl_matrix i j -- *[i,j] ) gsl_matrix_ptr ; : ]]! ( *gsl_matrix i j df -- ) gsl_matrix_set ; : ]]fill ( addr df -- ) gsl_matrix_set_all ; : ]]erase ( addr -- ) gsl_matrix_set_zero ; : ]]size1 gsl_matrix size1 @ ; : ]]size2 gsl_matrix size2 @ ; : ]]dim ( *gsl_matrix -- m n ) dup ]]size1 swap ]]size2 ; : ]]dim. ( *gsl_matrix -- ) ]]dim swap . ." x" . cr ; : ]]data ( *gsl_matrix -- addr) gsl_matrix data @ ; : ]]tda gsl_matrix tda @ ; : ]]block gsl_matrix block @ ; : ]]owner gsl_matrix owner @ ; : ]]copy]] ( *gsl_matrix_dest *gsl_matrix_src -- ) gsl_matrix_memcpy drop ; : ]]'copy]] ( *gsl_matrix_dest *gsl_matrix_src -- ) gsl_matrix_transpose_memcpy drop ; \ : ]]row ( *gsl_matrix idx -- *gsl_vector ) gsl_matrix_row ; \ : ]]col ( *gsl_matrix idx -- *gsl_vector ) gsl_matrix_column ; : ]]>] ( *gsl_vector *gsl_matrix i -- ) gsl_matrix_get_col drop ; : ]]>]' ( *gsl_vector *gsl_matrix i -- ) gsl_matrix_get_row drop ; : ]>]] ( *gsl_matrix *gsl_vector i -- ) swap gsl_matrix_set_col drop ; : ]'>]] ( *gsl_matrix *gsl_vector i -- ) swap gsl_matrix_set_row drop ; : ]]max gsl_matrix_max ; : ]]min gsl_matrix_min ; : ]]add! ( *gsl_matrix *gsl_matrix -- ) gsl_matrix_add drop ; : ]]sub! ( *gsl_matrix *gsl_matrix -- ) gsl_matrix_sub drop ; : ]]e*! ( *gsl_matrix *gsl_matrix -- ) gsl_matrix_mul_elements drop ; : ]]*c ( *gsl_matrix df -- ) gsl_matrix_scale drop ; : ]]+c ( df *gsl_matrix -- ) gsl_matrix_add_constant drop ; : ]]clone ( *gsl_matrix -- *gsl_matrix ) dup dup gsl_matrix size1 @ swap gsl_matrix size2 @ gsl_matrix_alloc dup -rot swap gsl_matrix_memcpy drop ; : ]]negate !-1.0 ]]*c ; : ]]+ ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) ]]clone dup -rot swap ]]add! ; : ]]- ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) swap ]]clone dup -rot swap ]]sub! ; : ]]null? ( *gsl_matrix -- 0/-1 ) gsl_matrix_isnull negate ; : ]]pos? ( *gsl_matrix -- 0/-1 ) gsl_matrix_ispos negate ; : ]]neg? ( *gsl_matrix -- 0/-1 ) gsl_matrix_isneg negate ; \ blas \ constants 101 Constant CblasRowMajor 102 Constant CblasColMajor 111 Constant CblasNoTrans 112 Constant CblasTrans 113 Constant CblasConjTrans 121 Constant CblasUpper 122 Constant CblasLower 131 Constant CblasNonUnit 132 Constant CblasUnit 141 Constant CblasLeft 142 Constant CblasRight : action? ( *gsl_matrix *gsl_matrix n n n -- ) dup 0= if drop 2swap 2dup ]]size2 swap ]]size1 swap exit then dup 1 = if drop 2swap 2dup ]]size2 swap ]]size2 swap exit then 2 = if 2swap 2dup ]]size1 swap ]]size1 swap exit then 3 = if 2swap 2dup ]]size1 swap ]]size2 swap exit then ; create samemattable ' noop , ' ]]clone , : samemat ( *gsl_matrix *gsl_matrix -- 1/0 *gsl_matrix ) dup -rot = abs dup -rot cells samemattable + @ execute ; macro : ]]mul ( *gsl_matrix *gsl_matrix n n n -- *gsl_matrix ) !1 !0 action? gsl_matrix_alloc dup >r gsl_blas_dgemm drop r> ; : ]]* ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) 2dup samemat dup rot 2>r nip CblasNoTrans CblasNoTrans 0 ]]mul 2r> if ]]free else drop then ; : ]]'* ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) 2dup samemat dup rot 2>r nip CblasTrans CblasNoTrans 1 ]]mul 2r> if ]]free else drop then ; : ]]*' ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) 2dup samemat dup rot 2>r nip CblasNoTrans CblasTrans 2 ]]mul 2r> if ]]free else drop then ; : ]]'*' ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) 2dup samemat dup rot 2>r nip CblasTrans CblasTrans 3 ]]mul 2r> if ]]free else drop then ; : ]]mul! ( n n *gsl_matrix *gsl_matrix *gsl_matrix -- ) !1 !0 gsl_blas_dgemm drop ; : ]]*! ( *gsl_matrix *gsl_matrix *gsl_matrix -- ) >r CblasNoTrans CblasNoTrans 2swap r> ]]mul! ; : ]]'*! ( *gsl_matrix *gsl_matrix *gsl_matrix -- ) >r CblasTrans CblasNoTrans 2swap r> ]]mul! ; : ]]*'! ( *gsl_matrix *gsl_matrix *gsl_matrix -- ) >r CblasNoTrans CblasTrans 2swap r> ]]mul! ; : ]]*] ( *gsl_matrix *gsl_vector -- *gsl_vector ) over ]]size1 gsl_vector_alloc >r CblasNoTrans -rot r@ !1 !0 gsl_blas_dgemv drop r> ; : ]]'*] ( *gsl_matrix *gsl_vector -- *gsl_vector ) over ]]size1 gsl_vector_alloc >r CblasTrans -rot r@ !1 !0 gsl_blas_dgemv drop r> ; : ]]i ( *gsl_matrix -- ) dup dup ]]size1 swap ]]size2 <> if abort" ERROR: Not a square matrix!" then dup ]]size1 0 do dup i i !1 ]]! loop drop ; : identity ( n -- *gsl_matrix ) dup gsl_matrix_calloc dup ]]i ; : min-identity ( *gsl_matrix -- *gsl_matrix ) dup ]]size1 swap ]]size2 min identity ; : left/right' ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) over ]]size1 over ]]size1 > if swap ]]*' exit else ]]'* exit then ; \ original matrix remains intact : ]]' ( *gsl_matrix -- *gsl_matrix ) dup min-identity dup >r left/right' r> ]]free ; : ]]T! ( *gsl_matrix -- ) gsl_matrix_transpose drop ; : ]]T ( *gsl_matrix -- *gsl_matrix ) dup ]]dim swap gsl_matrix_alloc dup rot gsl_matrix_transpose_memcpy drop ; : ]]2T ( *gsl_matr *gsl_matrix -- ) gsl_matrix_transpose_memcpy drop ; : ]]+! ( *gsl_matrix i j df -- ) >r 2dup r@ ]]@ f+ r> ]]! ; : ]]scale! ( *gsl_matrix i j df -- ) >r 2dup r@ ]]@ f* r> ]]! ; : ]]data_ij ( *gsl_matrix i j -- addr) rot >r swap r@ ]]tda dfloats * swap dfloats + r> ]]data + ; \ Cross product can be either calculated through determinant: : ]x ( *gsl_vector *gsl_vector -- *gsl_vector ) 3 gsl_vector_alloc { x1[ x2[ x3[ | x1[ 2 ]@ fnegate x2[ 1 ]@ f* x1[ 1 ]@ x2[ 2 ]@ f* f+ x3[ 0 ]! x1[ 2 ]@ x2[ 0 ]@ f* x1[ 0 ]@ fnegate x2[ 2 ]@ f* f+ x3[ 1 ]! x1[ 1 ]@ fnegate x2[ 0 ]@ f* x1[ 0 ]@ x2[ 1 ]@ f* f+ x3[ 2 ]! x3[ } ; \ or using algebraic form when first vector in the product is \ rewritten in a matrix form: \ a x b = [C_a] b \ [ 0 -a[2] a[1] ] \ [C_a] = [ a[2] 0 -a[0] ] \ [-a[1] a[0] 0 ] \ a function to convert a vector into such matrix: : ]>[x] ( ] -- ]] ) [IFDEF] debug dup ]size 3 <> abort" Not a 3D vector!" [THEN] 3 3 :]] dup >r swap 2dup 2 ]@ fdup dup fnegate 0 1 ]]! 1 0 ]]! 2dup 1 ]@ fdup dup fnegate 2 0 ]]! 0 2 ]]! 0 ]@ fdup dup fnegate 1 2 ]]! 2 1 ]]! r> ; : ]. ( *gsl_vector *gsl_vector -- f:dot_product ) { x1[ x2[ | 0 0 sp@ x1[ x2[ rot gsl_blas_ddot drop fd>f } ; : ]total ( *gsl_vector -- f:sum ) dup ]size gsl_vector_alloc dup !1 ]fill dup rot ]. ]free ; \ probability normalize - assures sum is unity : ]pnormalize ( *gsl_vector - ) dup ]total 1/f ]*c ; : |]| ( *gsl_vector -- f:norm ) dup ]. fsqrt ; \ assures vector norm is unity : ]normalize ( *gsl_vector - ) dup |]| 1/f ]*c ; : ]distance ( *gsl-vector *gsl-vector -- f ) ]sub dup |]| ]free ; : ]+! ( *gsl_vector i df -- ) 2dup ]@ f+ ]! ; : ]*! ( *gsl_vector i df -- ) 2dup ]@ f* ]! ; : ]]*]m ( *gsl_matrix *gsl_vector -- *gsl_vector ) over ]]size1 gsl_vector_calloc { m[[ x[ y[ | m[[ ]]size1 0 do m[[ ]]size2 0 do m[[ j i ]]@ x[ i ]@ f* y[ j ]+! loop loop y[ } ; : >#rows ( -- ) swap ]]size1 >= abort" number of rows is bigger than available!" ; : >#cols ( -- ) swap ]]size2 >= abort" number of columns is bigger than available!" ; : ]]row ( *gsl_matrix n -- *gsl_vector ) 2dup >#rows sizeof gsl_vector allocate throw dup 1 swap gsl_vector stride ! >r over ]]size2 r@ gsl_vector size ! 0 ]]data_ij r@ gsl_vector data ! 0 r@ gsl_vector owner ! r> ; \ assumes all dimensions are set correctly : ]]row! ( *gsl_vector *gsl_matrix n -- ) rot >r 2dup >#rows 0 ]]data_ij r> gsl_vector data ! ; : ]]col ( *gsl_matrix n -- *gsl_vector ) 2dup >#cols sizeof gsl_vector allocate throw >r over ]]tda r@ gsl_vector stride ! over ]]size1 r@ gsl_vector size ! over ]]block r@ gsl_vector block ! 0 swap ]]data_ij r@ gsl_vector data ! 0 r@ gsl_vector owner ! r> ; : ]]rfill ( f:n *gsl_matrix i -- ) ]]row dup ]fill ]free ; : ]]cfill ( f:n *gsl_matrix i -- ) ]]col dup ]fill ]free ; : ]]submat ( *gsl_matrix n1 n2 m1 m2 -- *gsl_matrix ) { m[[ n1 n2 m1 m2 | sizeof gsl_matrix allocate throw >r n2 n1 - 1+ r@ gsl_matrix size1 ! m2 m1 - 1+ r@ gsl_matrix size2 ! m[[ n1 m1 ]]data_ij r@ gsl_matrix data ! m[[ ]]tda r@ gsl_matrix tda ! 0 r@ gsl_matrix owner ! r> } ; : ?square ( *gsl_matrix -- ) dup ]]size1 swap ]]size2 <> abort" ERROR: Not a square matrix!" ; : ]]diag ( *gsl_matrix n1 n2 -- *gsl_vector ) rot dup ?square -rot sizeof gsl_vector allocate throw { d[ | over - d[ gsl_vector size ! 2dup dup ]]data_ij d[ gsl_vector data ! drop dup ]]tda d[ gsl_vector stride ! ]]block d[ gsl_vector block ! 0 d[ gsl_vector owner ! d[ } ; \ with input matrix replaced by the result : ]]gsl-svd ( *gsl_matrix -- *gsl_matrix *gsl_vector ) dup ]]size2 dup dup gsl_matrix_calloc swap dup gsl_vector_calloc swap gsl_vector_calloc { mV vS vW | mV vS vW gsl_linalg_SV_decomp drop vW ]free mV vS } ; \ seems to be 30% faster : ]]gsl-svdm ( *gsl_matrix -- *gsl_matrix *gsl_vector ) dup ]]size2 dup ( a n n -- ) dup dup gsl_matrix_calloc swap ( a n a n -- ) dup gsl_matrix_calloc rot dup ( a a a n n -- ) gsl_vector_calloc swap gsl_vector_calloc { mX mV vS vW | mX mV vS vW gsl_linalg_SV_decomp_mod drop vW ]free mX ]]free mV vS } ; : ]]alu ( *gsl_matrix -- *gsl_permutation ) ( matrix replaced with its lu ) { a[[ | CblasRowMajor a[[ ]]size1 a[[ ]]size2 a[[ ]]data a[[ ]]size1 dup gsl_permutation_alloc dup >r }data clapack_dgetrf throw r> } ; \ the following version is for pure lapack when no ATLAS is available \ : ]]alu ( *gsl_matrix -- *gsl_permutation ) ( matrix replaced with its lu ) \ 0 0 0 3 ivector* { a[[ v( | \ a[[ ]]dim v( 1 )! v( 0 )! \ v( 0 *) v( 1 *) a[[ ]]data \ v( 1 *) v( 0 )@ gsl_permutation_alloc dup >r }data \ v( 2 *) dgetrf r> v( )free } ; : ]]ainv ( *gsl_matrix *gsl_permutation -- ) \ LU of a matrix replaced with its inverse { a[[ t{ | CblasRowMajor a[[ ]]size2 a[[ ]]data a[[ ]]size1 t{ }data clapack_dgetri throw } ; : ]]ainvert ( *gsl_matrix -- *gsl_matrix ) [IFDEF] отладка dup ?square [THEN] ]]clone dup dup >r ]]alu dup >r ]]ainv r> }free r> ; : ]]det ( *gsl_matrix -- f:determinant ) [IFDEF] отладка dup ?square [THEN] ]]clone dup ]]alu >r 1e0 dup ]]size1 0 do dup i dup ]]@ f* loop ]]free \ compute permutation sign r> }sign s>f f* }free ; \ calculates the work needed for dgesvd ( see man dgesvd ) : lwork ( m n -- c ) 2dup max -rot min 3 * over + swap 5 * max ; \ this svd returns U MxM so eats a lot of memory : ]]asvda ( *gsl_matrix -- *gsl_matrix *gsl_matrix *gsl_vector ) ]]clone { A[[ | A[[ ]]size1 dup gsl_matrix_alloc A[[ ]]size2 dup gsl_matrix_alloc A[[ ]]size1 A[[ ]]size2 min gsl_vector_alloc 8 cells allocate throw { U[[ V[[ W[ p[ | ascii A p[ 0 cells + ! p[ 0 cells + ascii A p[ 1 cells + ! p[ 1 cells + A[[ ]]size1 p[ 2 cells + ! p[ 2 cells + A[[ ]]size2 p[ 3 cells + ! p[ 3 cells + A[[ ]]data p[ 2 cells + W[ ]data U[[ ]]data U[[ ]]size1 p[ 4 cells + ! p[ 4 cells + V[[ ]]data V[[ ]]size1 p[ 5 cells + ! p[ 5 cells + A[[ ]]size1 A[[ ]]size2 lwork dup gsl_vector_alloc dup >r ]data swap p[ 6 cells + ! p[ 6 cells + p[ 7 cells + dgesvd r> ]free p[ free throw A[[ ]]free U[[ V[[ W[ } } ; \ performs A=U*S*V^T \ A = MxN, where M>N, pass it A^T \ returns U^T (MxN), V(NxN) and vector of N eigenvalues : ]]asvdO ( *gsl_matrix -- *gsl_matrix *gsl_matrix *gsl_vector ) { A[[ | A[[ ]]size2 A[[ ]]size1 min dup gsl_matrix_alloc A[[ ]]size1 A[[ ]]size2 min gsl_vector_alloc 8 cells allocate throw { V[[ W[ p[ | ascii O p[ 0 cells + ! p[ 0 cells + ascii S p[ 1 cells + ! p[ 1 cells + A[[ ]]size2 p[ 2 cells + ! p[ 2 cells + A[[ ]]size1 p[ 3 cells + ! p[ 3 cells + A[[ ]]data p[ 2 cells + W[ ]data 0 p[ 2 cells + V[[ ]]data V[[ ]]size2 p[ 5 cells + ! p[ 5 cells + A[[ ]]size2 A[[ ]]size1 lwork dup gsl_vector_alloc dup >r ]data swap p[ 6 cells + ! p[ 6 cells + p[ 7 cells + dgesvd r> ]free p[ free throw A[[ V[[ W[ } } ; \ in place Cholesky decomposition | 4 cell 8 * vector chol( \ a vector that simplifies parameter passing : choldc-errors ( i -- ) dup 0= if drop exit then dup 0< if abs ." Argument " . ." is set incorrectly!" abort then ." Leading minor of order " . ." is not positive definite!" abort ; : ]]choldc! ( *gsl_matrix -- int ) chol( 0 *) swap dup ]]dim chol( 2 )! chol( 1 )! ]]data chol( 1 *) swap chol( 2 *) chol( 3 *) dpotrf chol( 3 )@ ; \ A = U^T U : ]]choldcU! ( *gsl_matrix -- int ) 'U chol( 0 )! ]]choldc! ; \ A = L L^T : ]]choldcL! ( *gsl_matrix -- int ) 'L chol( 0 )! ]]choldc! ; \ Cholesky preserving the original matrix : ]]choldcU ( *gsl_matrix -- *gsl_matrix ) ]]clone dup ]]choldcu! choldc-errors ; : ]]choldcL ( *gsl_matrix -- *gsl_matrix ) ]]clone dup ]]choldcL! choldc-errors ; : ]diag[[ ( *gsl_vector -- *gsl_matrix ) dup ]size dup dup gsl_matrix_calloc swap 0 do 2dup swap i ]@ i i ]]! loop nip ; : ]print\ ( *gsl_vector -- ) dup ]size 0 do dup i ]@ fx. loop drop ; : ]print ( *gsl_vector -- ) precision swap 5 set-precision ]print\ cr set-precision ; : ]]print ( *gsl_matrix -- ) cr precision swap 5 set-precision dup ]]size1 0 do \ i . ." : " dup ]]size2 0 do dup j i ]]@ fs. loop cr loop drop set-precision ; : ]]row-print ( *gsl_matrix i -- ) cr over gsl_matrix size2 @ 0 do 2dup i ]]@ f. loop cr 2drop ; : ]]col-print ( *gsl_matrix i -- ) cr over gsl_matrix size1 @ 0 do 2dup i swap ]]@ f. loop cr 2drop ; : ]]nthrow ( *gsl_matrix n -- addr ) over ]]tda * dfloats swap ]]data + ; : ]]randomize ( *gsl_matrix -- ) dup dup ]]size1 swap ]]size2 * 0 do dup gsl-randomu ]]data i dfloats + df! loop drop ; : ]]randn ( *gsl_matrix -- ) dup dup ]]size1 swap ]]size2 * 0 do dup gsl-gaussian ]]data i dfloats + df! loop drop ; ' ]]randomize alias ]]randu : ]randomize ( *gsl_vector -- ) dup ]size 0 do dup gsl-randomu i ]! loop drop ; : ]randn ( *gsl_vector -- ) \ zero mean unit variance dup ]size 0 do dup gsl-gaussian i ]! loop drop ; ' ]randomize alias ]randu : ]mean ( *gsl_vector -- f ) dup ]stride swap dup ]size swap ]data rot rot gsl_stats_mean ; : ]variance ( *gsl_vector -- f ) dup ]stride swap dup ]size swap ]data rot rot gsl_stats_variance ; : ]sd ( *gsl_vector -- f ) dup ]stride swap dup ]size swap ]data rot rot gsl_stats_sd ; : ]skew ( *gsl_vector -- f ) dup ]stride swap dup ]size swap ]data rot rot gsl_stats_skew ; : ]kurtosis ( *gsl_vector -- f ) dup ]stride swap dup ]size swap ]data rot rot gsl_stats_kurtosis ; \ Add Gaussian noise to a vector. Standard deviation of the noise \ is % of vectors variance : %>r ( -- ) 100e f/ ; \ % to ratio : ]jstdv ( f:ratio m[ -- f:stdv ) dup ]variance fdup f0= if dup ]max f+ then f* drop ; \ for constant vectors : ]rnd-clone ( f:stdv m[ -- c[ ) ]clone dup ]randn dup ]*c ; : ]jitter+ ( f:% v[ -- ) %>r dup ]jstdv dup ]rnd-clone >r r@ ]+ r> ]free ; : ]]gsl-lu ( *gsl_matrix -- *gsl_matrix *gsl_permutation ) 1 sp@ rot ]]clone dup >r dup ]]size1 gsl_permutation_calloc dup >r rot gsl_linalg_LU_decomp drop r> r> swap rot drop ; : ]]gsl-invert ( *gsl_matrix -- *gsl_matrix ) ]]clone dup dup ]]gsl-lu 2dup >r >r rot gsl_linalg_LU_invert drop r> ]]free r> }free ; ' ]]ainvert alias ]]invert ' ]]asvdO alias ]]svd : ]]save ( *gsl_matrix *gsl_matrix_cfa fid -- ) -rot { m[[ name[[ | >r name[[ >name count 1+ nip 0 m[[ ]]size2 m[[ ]]size1 0 sp@ 5 cells r@ write-file throw 2drop 2drop drop name[[ >name count 1+ r@ write-file throw m[[ ]]size1 m[[ ]]size2 * dfloats m[[ ]]T dup s>f ]]data swap r> write-file throw f>s ]]free } ; \ these words do not work with float matrices but are needed for \ scientific calculations, that's why they are in this module : _hmatrix ( n m size -- addr ) rot over * 2 pick * [ 2 cells ] literal + allocate throw dup [ 2 cells ] literal + >r rot over ! [ 1 cells ] literal + ! r> ; : hmatrix ( n m size -- ) create rot over * 2 pick * [ 2 cells ] literal + allocate throw dup , rot over ! [ 1 cells ] literal + ! does> @ [ 2 cells ] literal + ; : }}row-size ( hmatrix -- ) [ 2 cells ] literal - @ ; : freeHmatrix ( hmatrix -- ) [ 2 cells ] literal - free throw ; : }} ( addr i j -- addr[i][j] ) \ word to fetch 2-D array addresses >R >R \ indices to return stack temporarily DUP CELL- CELL- 2@ \ &a[0][0] size m R> * R> + * + ALIGNED ; : h->[[ ( hmatrix -- gsl_matrix ) dup }}row-size 3 swap gsl_matrix_alloc dup ]]size2 0 do 3 0 do 2dup swap i j }} w@ s>f i j ]]! loop loop nip ; \ some sequencing code : arange ( f:start f:end f:step -- x[ ) f-rot fswap fdup f>r f- fover f/ f>s :] fr> dup ]size 0 do dup fover i s>f f* fover f+ i ]! loop fdrop fdrop ; : product ( x[ -- f:P ) !1 dup ]size 0 do dup i ]@ f* loop drop ; \ initializing random number generator to some value in order to have \ it available upon loading of gsl mt19937 \ choose n elements out of m possibilities : rnd-choose ( n m -- s[ ) gsl_rng_default -rot 0e s>f 1e arange dup >r over :] dup >r ]data -rot dup ]data swap ]size 1 dfloats gsl_ran_choose throw r> r> ]free ; : )randperm ( *v( -- ) gsl_rng_default swap dup )size over )type 8 / gsl_ran_shuffle ; : )>] ( v( -- v[ ) dup )size :] over )size 0 do 2dup swap i )@ s>f i ]! loop nip ; previous previous previous previous previous Module;