1\ GNU Scientific Library interface Mon Sep 12 14:40:15 MDT 2005 2\ Copyright (C) 2007, Sergey Plis 3\ 4\ This program is free software; you can redistribute it and/or modify 5\ it under the terms of the GNU General Public License as published by 6\ the Free Software Foundation; either version 2 of the License, or 7\ (at your option) any later version. 8\ 9\ This program is distributed in the hope that it will be useful, 10\ but WITHOUT ANY WARRANTY; without even the implied warranty of 11\ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12\ GNU General Public License for more details. 13 14\needs float import float 15warning off 16\needs locals| include locals.fs 17warning on 18\needs atlas include atlas.fs 19\needs callback include callback.fs 20\needs vectors include vectors.fs 21\needs complex include complex.fb 22 23Module GSL 24 25\ stole the hash function from hash.fs of bigforth 26| &14 Value Hashbits 27| 1 Hashbits << Value Hashlen 28 29Label (hash ( SI:string -- AX:key ) :R DX push 30 .b lods $1F # AX and AX CX mov DX DX xor CX 1 # shr 31 b IF SI ) AH mov SI inc THEN CX dec 32 0>= IF BEGIN .w SI ) DX mov 2 # SI add CX dec 33 DX AX *2 I) AX lea 0< UNTIL THEN 34 & Hashbits A#) CX mov AX DX mov AX shr DX AX add 35 & Hashlen A#) CX mov CX dec CX AX and DX pop ret 36| Code Hash ( string -- key ) 37 R: SI push AX SI mov (hash rel) call SI pop 38 Next end-code 39 40also dos also complex also float also atlas also vectors 41 42s" libptcblas.so" getlib 0<> 43[IF] 44 library libblas libptcblas.so depends libatlas 45[ELSE] 46 library libblas libcblas.so depends libatlas 47[THEN] 48 49library libgsl libgsl.so.0 depends libblas 50 51legacy off 52 53\ some functions 54libgsl gsl_log1p df (fp) gsl_log1p ( df -- df ) 55libgsl gsl_acosh df (fp) gsl_acosh ( df -- df ) 56 57\ error handling Wed Sep 21 23:04:06 MDT 2005 58libgsl gsl_set_error_handler ptr (int) gsl_set_error_handler 59( function -- function ) 60libgsl gsl_strerror int (ptr) gsl_strerror 61 62callback 4:0 (void) int int int int callback; 63: cstr-fstr ( addr -- addr len ) 64 0 65 begin 2dup + c@ 0 = not while 66 1+ 67 repeat ; 68 69| : .bold-red ." [1;31;40m" ; 70| : .red ." [2;31;40m" ; 71| : .reset ." [0;37;40m" ; 72| : cb-test 73 cr 74 \ .bold-red 75 ." GSL ERROR: " cr 76 \ .reset cr 77 10 spaces gsl_strerror cstr-fstr type cr 78 drop \ ." at line: " . cr 79 drop \ ." of file: " cstr-fstr type cr 80 10 spaces cstr-fstr type cr 81 \ .red 82 -1 abort" failed at" ; 83' cb-test 4:0 c_plus 84 85\ 1 2 c_plus 2:1call . 86| variable old_handler 87c_plus gsl_set_error_handler old_handler ! 88 890 Constant GSL_SUCCESS 90 91\ random number generation Mon Sep 12 22:06:01 MDT 2005 92 93libgsl gsl_rng_types_setup (ptr) gsl_rng_types_setup ( -- *gsl_rng_type) 94libgsl gsl_rng_env_setup (ptr) gsl_rng_env_setup ( -- *gsl_rng) 95libgsl gsl_rng_alloc int (int) gsl_rng_alloc ( *gsl_rng_type -- *gsl_rng ) 96libgsl gsl_rng_name int (int) gsl_rng_name ( *gsl_rng -- string ) 97libgsl gsl_rng_set int int (void) gsl_rng_set ( *gsl_rng int -- ) 98libgsl gsl_rng_uniform int (fp) gsl_rng_uniform ( *gsl_rng -- df ) 99libgsl gsl_rng_uniform_pos int (fp) gsl_rng_uniform_pos ( *gsl_rng -- df ) 100libgsl gsl_rng_uniform_int int int (int) gsl_rng_uniform_int ( *gsl_rng n --n ) 101libgsl gsl_rng_get int (int) gsl_rng_get ( *gsl_rng -- int ) 102libgsl gsl_rng_max int (int) gsl_rng_max ( *gsl_rng -- int ) 103libgsl gsl_rng_min int (int) gsl_rng_min ( *gsl_rng -- int ) 104libgsl gsl_rng_clone int (int) gsl_rng_clone ( *gsl_rng -- *gsl_rng ) 105libgsl gsl_rng_free int (int) gsl_rng_free ( *gsl_rng -- ) 106 107 108 109\ random number distributions Tue Sep 13 00:44:35 MDT 2005 110\ Gaussian 111libgsl gsl_ran_gaussian df int (fp) gsl_ran_gaussian ( *gsl_rng df -- df ) 112libgsl gsl_ran_gaussian_ratio_method df int (fp) gsl_ran_gaussian_ratio_method ( *gsl_rng df -- df ) 113libgsl gsl_ran_gaussian_pdf df df (fp) gsl_ran_gaussian_pdf ( df df -- df ) 114\ sigma = 1 115libgsl gsl_ran_ugaussian int (fp) gsl_ran_ugaussian ( *gsl_rng -- df ) 116libgsl gsl_ran_ugaussian_ratio_method int (fp) gsl_ran_ugaussian_ratio_method ( *gsl_rng -- df ) 117libgsl gsl_ran_ugaussian_pdf df (fp) gsl_ran_ugaussian_pdf ( df df -- df ) 118libgsl gsl_ran_discrete_preproc int int (int) gsl_ran_discrete_preproc ( int int -- int ) 119libgsl gsl_ran_discrete int int (int) gsl_ran_discrete 120libgsl gsl_ran_discrete_free int (void) gsl_ran_discrete_free 121libgsl gsl_ran_shuffle int int ptr ptr (void) gsl_ran_shuffle 122\ cdf P(x) = \int_{-\infty}^{x} p(x)dx Q(x) = \int_{x}^{\infty} p(x)dx 123libgsl gsl_cdf_gaussian_P df df (fp) gsl_cdf_gaussian_P ( df df -- df ) 124libgsl gsl_cdf_gaussian_Q df df (fp) gsl_cdf_gaussian_Q ( df df -- df ) 125libgsl gsl_cdf_gaussian_Pinv df df (fp) gsl_cdf_gaussian_Pinv ( df df -- df ) 126libgsl gsl_cdf_gaussian_Qinv df df (fp) gsl_cdf_gaussian_Qinv ( df df -- df ) 127\ sigma = 1 cdf 128libgsl gsl_cdf_ugaussian_P df (fp) gsl_cdf_ugaussian_P ( df -- df ) 129libgsl gsl_cdf_ugaussian_Q df (fp) gsl_cdf_ugaussian_Q ( df -- df ) 130libgsl gsl_cdf_ugaussian_Pinv df (fp) gsl_cdf_ugaussian_Pinv ( df -- df ) 131libgsl gsl_cdf_ugaussian_Qinv df (fp) gsl_cdf_ugaussian_Qinv ( df -- df ) 132 133 134\ statistics Tue Sep 13 01:17:35 MDT 2005 135libgsl gsl_stats_mean int int int (fp) gsl_stats_mean ( array{ step size -- df ) 136libgsl gsl_stats_variance int int int (fp) gsl_stats_variance ( array{ step size -- df ) 137libgsl gsl_stats_variance_m df int int int (fp) gsl_stats_variance_m ( df array{ step size -- df ) 138libgsl gsl_stats_sd int int int (fp) gsl_stats_sd ( array{ step size -- df ) 139libgsl gsl_stats_sd_m df int int int (fp) gsl_stats_sd_m ( df array{ step size -- df ) 140libgsl gsl_stats_skew int int int (fp) gsl_stats_skew ( array{ step size -- df ) 141libgsl gsl_stats_kurtosis int int int (fp) gsl_stats_kurtosis ( array{ step size -- df ) 142libgsl gsl_stats_lag1_autocorrelation int int int (fp) gsl_stats_lag1_autocorrelation 143( array{ step size -- df ) 144libgsl gsl_stats_max int int int (fp) gsl_stats_max ( array{ step size -- df ) 145libgsl gsl_stats_min int int int (fp) gsl_stats_min ( array{ step size -- df ) 146libgsl gsl_stats_max_index int int int (int) gsl_stats_max_index ( array{ step size -- n ) 147libgsl gsl_stats_min_index int int int (int) gsl_stats_min_index ( array{ step size -- n ) 148 149\ vectors and matrices Wed Sep 14 00:15:36 MDT 2005 150 151\ Vectors 152libgsl gsl_block_alloc int (int) gsl_block_alloc ( n -- addr ) 153libgsl gsl_block_calloc int (int) gsl_block_calloc ( n -- addr ) 154libgsl gsl_block_free int (int) gsl_block_free ( n -- addr ) 155 156libgsl gsl_vector_alloc int (int) gsl_vector_alloc ( n -- addr ) 157libgsl gsl_vector_calloc int (int) gsl_vector_calloc ( n -- addr ) 158libgsl gsl_vector_alloc_from_vector int int int ptr (int) gsl_vector_alloc_from_vector 159libgsl gsl_vector_free int (void) gsl_vector_free ( addr -- ) 160libgsl gsl_vector_get int int (fp) gsl_vector_get ( addr i -- df ) 161libgsl gsl_vector_set df int int (void/fp) gsl_vector_set ( df addr i -- ) 162libgsl gsl_vector_set_all df int (void) gsl_vector_set_all ( df addr -- ) 163libgsl gsl_vector_set_zero int (void) gsl_vector_set_zero ( addr -- ) 164libgsl gsl_vector_memcpy int int (int) gsl_vector_memcpy ( dest_addr src_addr -- n ) 165 166libgsl gsl_vector_add int int (int) gsl_vector_add ( addr addr -- n ) 167libgsl gsl_vector_sub int int (int) gsl_vector_sub ( addr addr -- n ) 168libgsl gsl_vector_mul int int (int) gsl_vector_mul ( addr addr -- n ) 169libgsl gsl_vector_div int int (int) gsl_vector_div ( addr addr -- n ) 170libgsl gsl_vector_scale df int (int) gsl_vector_scale ( df addr -- n ) 171libgsl gsl_vector_add_constant df int (int) gsl_vector_add_constant ( df addr -- n ) 172libgsl gsl_vector_max int (fp) gsl_vector_max ( addr -- df ) 173libgsl gsl_vector_min int (fp) gsl_vector_min ( addr -- df ) 174libgsl gsl_vector_max_index int (fp) gsl_vector_max_index ( addr -- df ) 175libgsl gsl_vector_min_index int (fp) gsl_vector_min_index ( addr -- df ) 176libgsl gsl_vector_subvector int int int (int) gsl_vector_subvector 177\ Vector properties 178libgsl gsl_vector_isnull ptr (int) gsl_vector_isnull 179libgsl gsl_vector_ispos ptr (int) gsl_vector_ispos 180libgsl gsl_vector_isneg ptr (int) gsl_vector_isneg 181 182\ permutations 183libgsl gsl_permutation_alloc int (int) gsl_permutation_alloc ( n -- *gsl_prm) 184libgsl gsl_permutation_calloc int (int) gsl_permutation_calloc ( n -- *gsl_prm) 185libgsl gsl_permutation_init int (void) gsl_permutation_init ( *gsl_prm -- ) 186libgsl gsl_permutation_free int (void) gsl_permutation_free ( *gsl_prm -- ) 187libgsl gsl_permutation_get int int (int) gsl_permutation_get ( *gsl_prm i -- n) 188 189\ Matrices 190\ Allocation 191libgsl gsl_matrix_alloc int int (int) gsl_matrix_alloc 192libgsl gsl_matrix_calloc int int (int) gsl_matrix_calloc 193libgsl gsl_matrix_alloc_from_block [ 5 ] ints (int) gsl_matrix_alloc_from_block 194libgsl gsl_matrix_alloc_from_matrix [ 5 ] ints (int) gsl_matrix_alloc_from_matrix 195libgsl gsl_matrix_free ( *gsl_matrix -- ) int (void) gsl_matrix_free 196\ Accessing matrix elements 197libgsl gsl_matrix_get int int int (fp) gsl_matrix_get ( *m i j -- df ) 198libgsl gsl_matrix_set df int int int (void) gsl_matrix_set ( df *m i j -- ) 199libgsl gsl_matrix_ptr int int int (int) gsl_matrix_ptr ( *m i j -- *[i,j] ) 200\ Initializing matrix elements 201libgsl gsl_matrix_set_all df int (void) gsl_matrix_set_all ( *m df -- n ) 202libgsl gsl_matrix_set_zero df int (void) gsl_matrix_set_zero ( *m df -- n ) 203libgsl gsl_matrix_set_identity df int (void) gsl_matrix_set_identity ( *m df -- n ) 204\ Reading and writing matrices 205libgsl gsl_matrix_fwrite ptr ptr (int) gsl_matrix_fwrite 206libgsl gsl_matrix_fread ptr ptr (int) gsl_matrix_fread 207libgsl gsl_matrix_fprintf ptr ptr ptr (int) gsl_matrix_fprintf 208libgsl gsl_matrix_fscanf ptr ptr (int) gsl_matrix_fscanf 209\ Copying matrices 210libgsl gsl_matrix_memcpy int int (int) gsl_matrix_memcpy ( *m *m -- n ) 211libgsl gsl_matrix_swap int int (int) gsl_matrix_swap ( *m *m -- n ) 212\ Copying Rows and columns 213libgsl gsl_matrix_get_row int int int (int) gsl_matrix_get_row 214libgsl gsl_matrix_set_row int int int (int) gsl_matrix_set_row 215libgsl gsl_matrix_get_col int int int (int) gsl_matrix_get_col 216libgsl gsl_matrix_set_col int int int (int) gsl_matrix_set_col 217\ Exchanging rows and columns 218libgsl gsl_matrix_swap_rows int int ptr (int) gsl_matrix_swap_rows 219libgsl gsl_matrix_swap_columns int int ptr (int) gsl_matrix_swap_columns 220libgsl gsl_matrix_swap_rowcol int int ptr (int) gsl_matrix_swap_rowcol 221libgsl gsl_matrix_transpose_memcpy int int (int) gsl_matrix_transpose_memcpy 222libgsl gsl_matrix_transpose int (int) gsl_matrix_transpose 223\ Matrix operations 224libgsl gsl_matrix_add int int (int) gsl_matrix_add 225libgsl gsl_matrix_sub int int (int) gsl_matrix_sub 226libgsl gsl_matrix_mul_elements int int (int) gsl_matrix_mul_elements 227libgsl gsl_matrix_div_elements int int (int) gsl_matrix_div_elements 228libgsl gsl_matrix_scale df int (int) gsl_matrix_scale 229libgsl gsl_matrix_add_constant df int (int) gsl_matrix_add_constant 230\ Finding maximum and minimum elements of matrices 231libgsl gsl_matrix_max ptr (fp) gsl_matrix_max 232libgsl gsl_matrix_min ptr (fp) gsl_matrix_min 233libgsl gsl_matrix_minmax ptr ptr ptr (void) gsl_matrix_minmax 234libgsl gsl_matrix_min_index ptr ptr ptr (void) gsl_matrix_min_index 235libgsl gsl_matrix_max_index ptr ptr ptr (void) gsl_matrix_max_index 236libgsl gsl_matrix_minmax_index ptr ptr ptr ptr ptr (void) gsl_matrix_minmax_index 237\ Matrix properties 238libgsl gsl_matrix_isnull ptr (int) gsl_matrix_isnull 239libgsl gsl_matrix_ispos ptr (int) gsl_matrix_ispos 240libgsl gsl_matrix_isneg ptr (int) gsl_matrix_isneg 241\ libgsl gsl_matrix_isnonneg ptr (int) gsl_matrix_isnonneg 242 243 244libgsl gsl_matrix_submatrix int int int int int (int) gsl_matrix_submatrix ( *gsl_matrix k1 k2 n1 n2 -- n ) 245libgsl gsl_matrix_row int int (int) gsl_matrix_row ( *gsl_matrix idx -- *gsl_vector ) 246libgsl gsl_matrix_column int int (int) gsl_matrix_column ( *gsl_matrix idx -- *gsl_vector ) 247libgsl gsl_matrix_diagonal int (int) gsl_matrix_diagonal ( *gsl_matrix -- *gsl_vector ) 248 249 250\ BLAS Wed Sep 14 16:10:34 MDT 2005 251\ libblas cblas_dgemm int int df int int int 252\ int df int int int int int int (void/fp) cblas_dgemm 253libblas cblas_dgemv int int int int df int 254int df int int int int (void/fp) cblas_dgemv 255libgsl gsl_blas_ddot int int int (int) gsl_blas_ddot 256( *gsl_vector *gsl_vector df -- n ) 257libgsl gsl_blas_dgemm int df int int df int int (int/fp) gsl_blas_dgemm 258libgsl gsl_blas_dger int int int df (int/fp) gsl_blas_dger 259( alpha *gsl_vector *gsl_vector *gsl_matrix -- n ) ( A=\alpha x y^T+A ) 260libgsl gsl_blas_dgemv int df int int df int (int/fp) gsl_blas_dgemv 261( n alpha *gsl_matrix *gsl_vector beta *gsl_vector -- n ) 262 263\ Linear ALgebra Wed Sep 14 13:39:22 MDT 2005 264 265libgsl gsl_linalg_LU_decomp int int int (int) gsl_linalg_LU_decomp 266( *gsl_matrix *gsl_permutation *variable -- n ) 267libgsl gsl_linalg_LU_invert int int int (int) gsl_linalg_LU_invert 268( *gsl_matrix *gsl_permutation *gsl_matrix -- n ) 269libgsl gsl_linalg_SV_decomp int int int int (int) gsl_linalg_SV_decomp 270( *gsl_matrix *gsl_matrix *gsl_vector *gsl_vector -- n ) 271libgsl gsl_linalg_SV_decomp_mod int int int int int (int) gsl_linalg_SV_decomp_mod 272( *gsl_matrix *gsl_matrix *gsl_matrix *gsl_vector *gsl_vector -- n ) 273 274\ ----------------------------------------------------------------------------- 275\ *** Ordinary Differential Equations *** 276\ --- ODE system 277struct{ 278 cell func \ (* function) 279 \ (double t, const double y[], double dydt[], void * params); 280 cell jac \ (* jacobian) 281 \ (double t, const double y[], double * dfdy, double dfdt[], 282 \ void * params); 283 cell dim \ dimension; 284 cell params \ * params; 285} gsl_odeiv_system 286\ constants related to ODE 287 1 constant GSL_ODEIV_HADJ_INC 288 0 constant GSL_ODEIV_HADJ_NIL 289-1 constant GSL_ODEIV_HADJ_DEC 290 291callback gsl_odeiv_func4:1 (int) df int int int callback; 292callback gsl_odeiv_jac5:1 (int) df int int int int callback; 293 294\ --- Stepping Functions 295libgsl gsl_odeiv_step_alloc ptr int (ptr) gsl_odeiv_step_alloc 296( *step_type int -- *step ) 297libgsl gsl_odeiv_step_reset ptr (int) gsl_odeiv_step_reset ( *step -- r ) 298libgsl gsl_odeiv_step_free ptr (void) gsl_odeiv_step_free ( *step -- ) 299libgsl gsl_odeiv_step_name ptr (ptr) gsl_odeiv_step_name ( *step -- *str0 ) 300libgsl gsl_odeiv_step_order ptr (int) gsl_odeiv_step_order ( *step -- order) 301libgsl gsl_odeiv_step_apply int int int int int df df int (int) gsl_odeiv_step_apply 302( -- ) 303\ --- Available algorithms 304libgsl _gsl_odeiv_step_rk2 (int) gsl_odeiv_step_rk2 305libgsl _gsl_odeiv_step_rk4 (int) gsl_odeiv_step_rk4 306libgsl _gsl_odeiv_step_rkf45 (int) gsl_odeiv_step_rkf45 307libgsl _gsl_odeiv_step_rkck (int) gsl_odeiv_step_rkck 308libgsl _gsl_odeiv_step_rk8pd (int) gsl_odeiv_step_rk8pd 309libgsl _gsl_odeiv_step_rk2imp (int) gsl_odeiv_step_rk2imp 310libgsl _gsl_odeiv_step_rk4imp (int) gsl_odeiv_step_rk4imp 311libgsl _gsl_odeiv_step_bsimp (int) gsl_odeiv_step_bsimp 312libgsl _gsl_odeiv_step_gear1 (int) gsl_odeiv_step_gear1 313libgsl _gsl_odeiv_step_gear2 (int) gsl_odeiv_step_gear2 314 315: gsl_odeiv_step_rk2 [func'] _gsl_odeiv_step_rk2 @ ; 316: gsl_odeiv_step_rk4 [func'] _gsl_odeiv_step_rk4 @ ; 317: gsl_odeiv_step_rkf45 [func'] _gsl_odeiv_step_rkf45 @ ; 318: gsl_odeiv_step_rkck [func'] _gsl_odeiv_step_rkck @ ; 319: gsl_odeiv_step_rk8pd [func'] _gsl_odeiv_step_rk8pd @ ; 320: gsl_odeiv_step_rk2imp [func'] _gsl_odeiv_step_rk2imp @ ; 321: gsl_odeiv_step_rk4imp [func'] _gsl_odeiv_step_rk4imp @ ; 322: gsl_odeiv_step_bsimp [func'] _gsl_odeiv_step_bsimp @ ; 323: gsl_odeiv_step_gear1 [func'] _gsl_odeiv_step_gear1 @ ; 324: gsl_odeiv_step_gear2 [func'] _gsl_odeiv_step_gear2 @ ; 325 326\ --- Adaptive Step-size Control 327libgsl gsl_odeiv_control_standard_new df df df df (ptr) gsl_odeiv_control_standard_new ( a_dydt a_y eps_rel eps_abs -- *control ) 328libgsl gsl_odeiv_control_y_new df df (int) gsl_odeiv_control_y_new 329( eps_abs eps_rel -- *control ) 330libgsl gsl_odeiv_control_yp_new df df (ptr) gsl_odeiv_control_yp_new 331( eps_abs eps_rel -- *control ) 332libgsl gsl_odeiv_control_free ptr (void) gsl_odeiv_control_free ( *control -- ) 333libgsl gsl_odeiv_control_name ptr (ptr) gsl_odeiv_control_name ( *c -- *str0 ) 334 335\ --- Evolution 336libgsl gsl_odeiv_evolve_alloc int (int) gsl_odeiv_evolve_alloc 337( #dimensions -- evolution_func ) 338libgsl gsl_odeiv_evolve_apply int int df int int int int int (int) gsl_odeiv_evolve_apply 339( -- ) 340libgsl gsl_odeiv_evolve_reset ptr (int) gsl_odeiv_evolve_reset ( *e -- r ) 341libgsl gsl_odeiv_evolve_free ptr (void) gsl_odeiv_evolve_free ( *e -- ) 342\ ----------------------------------------------------------------------------- 343\ *** Fast Fourier Transform *** 344\ -- real 345libgsl gsl_fft_real_wavetable_alloc int (ptr) gsl_fft_real_wavetable_alloc 346libgsl gsl_fft_real_wavetable_free ptr (void) gsl_fft_real_wavetable_free 347libgsl gsl_fft_real_workspace_alloc int (ptr) gsl_fft_real_workspace_alloc 348libgsl gsl_fft_real_workspace_free ptr (void) gsl_fft_real_workspace_free 349\ in-place 350libgsl gsl_fft_real_transform ptr int int ptr ptr (int) gsl_fft_real_transform 351libgsl gsl_fft_real_unpack ptr ptr int int (int) gsl_fft_real_unpack 352 353\ -- halfcomplex 354\ - mixed radix 355libgsl gsl_fft_hc_wtbl_alloc int (ptr) gsl_fft_halfcomplex_wavetable_alloc 356libgsl gsl_fft_hc_wtbl_free ptr (void) gsl_fft_halfcomplex_wavetable_free 357libgsl gsl_fft_hc_backward ptr int int ptr ptr (int) gsl_fft_halfcomplex_backward 358libgsl gsl_fft_hc_inverse ptr int int ptr ptr (int) gsl_fft_halfcomplex_inverse 359libgsl gsl_fft_hc_transform ptr int int ptr ptr (int) gsl_fft_halfcomplex_transform 360libgsl gsl_fft_hc_unpack ptr ptr int int (int) gsl_fft_halfcomplex_unpack 361\ - radix2 362libgsl gsl_fft_hc_r2_unpack ptr ptr int int (int) gsl_fft_halfcomplex_radix2_unpack 363libgsl gsl_fft_hc_r2_backward ptr int int (int) gsl_fft_halfcomplex_radix2_backward 364libgsl gsl_fft_hc_r2_inverse ptr int int (int) gsl_fft_halfcomplex_radix2_inverse 365libgsl gsl_fft_hc_r2_transform ptr int int (int) gsl_fft_halfcomplex_radix2_transform 366 367 368| hashlen 32 vector fftpre( 369struct{ 370 cell next 371 cell size 372 cell workspace 373 cell r_wavetable 374 cell hc_wavetable 375} gsl_fft_precomputes 376| create $buf 255 allot 377| : 2str dup >r abs s>d <# #s r> sign #> $buf 0place ; 378| : s>hash ( n -- key ) 2str $buf hash ; 379| : (cache-fft) ( n -- addr ) 380 sizeof gsl_fft_precomputes allocate throw >r 381 0 r@ gsl_fft_precomputes next ! 382 dup r@ gsl_fft_precomputes size ! 383 dup gsl_fft_real_workspace_alloc r@ gsl_fft_precomputes workspace ! 384 dup gsl_fft_real_wavetable_alloc r@ gsl_fft_precomputes r_wavetable ! 385 gsl_fft_hc_wtbl_alloc r@ gsl_fft_precomputes hc_wavetable ! 386 r> ; 387| : cache-fft ( size -- addr ) 388 dup s>hash 389 fftpre( over )@ 0= if 390 swap (cache-fft) 391 fftpre( rot dup >r )! 392 fftpre( r> )@ 393 else 394 swap (cache-fft) 395 swap fftpre( over )@ 396 over gsl_fft_precomputes next ! 397 fftpre( rot dup >r )! 398 fftpre( r> )@ 399 then ; 400\ in case not found addr is just the size 401| : find-fft-cache ( n -- addr 0/1 ) 402 dup s>hash fftpre( swap )@ dup 403 begin while 404 2dup gsl_fft_precomputes size @ = 405 if nip true exit then 406 gsl_fft_precomputes next @ dup 407 repeat ; 408 409legacy on 410 411\ Structures 412 413struct{ 414 cell name 415 cell max 416 cell min 417 cell size 418 cell set 419 cell get 420 cell get_double 421} gsl_rng_type 422 423struct{ 424 cell type 425 cell state 426} gsl_rng 427 428struct{ 429 cell size 430 cell data 431} gsl_block 432 433struct{ 434 cell size 435 cell stride 436 cell data 437 cell block 438 cell owner 439} gsl_vector 440 441' gsl_block alias gsl_permutation 442 443struct{ 444 cell size1 445 cell size2 446 cell tda 447 cell data 448 cell block 449 cell owner 450} gsl_matrix 451 452\ random number generation functions 453: 0-len dup 1- 0 begin 1+ 2dup + c@ 0= until nip ; 454: )gsl-rng ( addr i -- *gsl_rng_type ) 455 cells + @ ; 456 457\ setting up all available random number generators 458gsl_rng_types_setup value gsl_rng_array( 4590 value gsl_rng_default 460: gsl-free ( -- ) 461 gsl_rng_default gsl_rng_free ; 462 463: borosh13 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 464 gsl_rng_array( 0 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 465: cmrg ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 466 gsl_rng_array( 1 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 467: coveyou ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 468 gsl_rng_array( 2 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 469: fishman18 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 470 gsl_rng_array( 3 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 471: fishman20 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 472 gsl_rng_array( 4 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 473: fishman2x ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 474 gsl_rng_array( 5 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 475: gfsr4 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 476 gsl_rng_array( 6 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 477: knuthran ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 478 gsl_rng_array( 7 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 479: knuthran2 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 480 gsl_rng_array( 8 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 481: lecuyer21 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 482 gsl_rng_array( 9 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 483: minstd ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 484 gsl_rng_array( 10 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 485: mrg ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 486 gsl_rng_array( 11 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 487: mt19937 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 488 gsl_rng_array( 12 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 489: mt19937_1999 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 490 gsl_rng_array( 13 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 491: mt19937_1998 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 492 gsl_rng_array( 14 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 493: r250 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 494 gsl_rng_array( 15 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 495: ran0 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 496 gsl_rng_array( 16 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 497: ran1 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 498 gsl_rng_array( 17 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 499: ran2 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 500 gsl_rng_array( 18 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 501: ran3 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 502 gsl_rng_array( 19 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 503: rand ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 504 gsl_rng_array( 20 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 505: rand48 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 506 gsl_rng_array( 21 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 507: random128-bsd ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 508 gsl_rng_array( 22 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 509: random128-glibc2 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 510 gsl_rng_array( 23 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 511: random128-libc5 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 512 gsl_rng_array( 24 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 513: random256-bsd ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 514 gsl_rng_array( 25 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 515: random256-glibc2 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 516 gsl_rng_array( 26 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 517: random256-libc5 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 518 gsl_rng_array( 27 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 519: random32-bsd ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 520 gsl_rng_array( 28 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 521: random32-glibc2 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 522 gsl_rng_array( 29 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 523: random32-libc5 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 524 gsl_rng_array( 30 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 525: random64-bsd ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 526 gsl_rng_array( 31 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 527: random64-glibc2 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 528 gsl_rng_array( 32 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 529: random64-libc5 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 530 gsl_rng_array( 33 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 531: random8-bsd ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 532 gsl_rng_array( 34 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 533: random8-glibc2 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 534 gsl_rng_array( 35 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 535: random8-libc5 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 536 gsl_rng_array( 36 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 537: random-bsd ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 538 gsl_rng_array( 37 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 539: random-glibc2 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 540 gsl_rng_array( 38 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 541: random-libc5 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 542 gsl_rng_array( 39 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 543: randu ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 544 gsl_rng_array( 40 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 545: ranf ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 546 gsl_rng_array( 41 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 547: ranlux ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 548 gsl_rng_array( 42 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 549: ranlux389 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 550 gsl_rng_array( 43 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 551: ranlxd1 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 552 gsl_rng_array( 44 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 553: ranlxd2 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 554 gsl_rng_array( 45 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 555: ranlxs0 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 556 gsl_rng_array( 46 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 557: ranlxs1 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 558 gsl_rng_array( 47 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 559: ranlxs2 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 560 gsl_rng_array( 48 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 561: ranmar ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 562 gsl_rng_array( 49 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 563: slatec ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 564 gsl_rng_array( 50 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 565: taus ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 566 gsl_rng_array( 51 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 567: taus2 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 568 gsl_rng_array( 52 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 569: taus113 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 570 gsl_rng_array( 53 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 571: transputer ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 572 gsl_rng_array( 54 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 573: tt800 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 574 gsl_rng_array( 55 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 575: uni ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 576 gsl_rng_array( 56 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 577: uni32 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 578 gsl_rng_array( 57 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 579: vax ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 580 gsl_rng_array( 58 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 581: waterman14 ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 582 gsl_rng_array( 59 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 583: zuf ( -- *gsl_rng ) gsl_rng_default 0<> if gsl-free then 584 gsl_rng_array( 60 )gsl-rng gsl_rng_alloc to gsl_rng_default ; 585 586\ words for actual generation of random numbers 587: gsl-randomg ( -- n ) 588 gsl_rng_default gsl_rng_get ; 589: gsl-randomu /* -- f \in [0,1) */ 590 gsl_rng_default gsl_rng_uniform ; 591: gsl-randomu+ /* -- f \in (0,1) */ 592 gsl_rng_default gsl_rng_uniform_pos ; 593: gsl-random-up ( n -- f \in [0,n] ) 594 gsl_rng_default swap gsl_rng_uniform_int ; 595: gsl-set-seed ( n -- ) 596 gsl_rng_default swap gsl_rng_set ; 597: gsl-clone ( -- *gsl_rng ) 598 gsl_rng_default gsl_rng_clone ; 599: gsl-gaussian ( -- df ) 600 gsl_rng_default !1 gsl_ran_gaussian ; 601: gsl-discrete ( *gsl_ran_discrete -- n ) 602 gsl_rng_default swap gsl_ran_discrete ; 603 604\ vectors and matrices 6050 constant GSL_VECTOR_TYPE 6061 constant GSL_MATRIX_TYPE 607: gsltypeof ' >body cell + @ ; 608 609: fvector ( n -- | -- id addr ) 610 create 611 gsl_vector_calloc , 612 GSL_VECTOR_TYPE , 613 does> @ ; 614 615\ allocate a nameless vector 616: :] ( # -- addr ) gsl_vector_calloc ; 617\ allocate a nameless matrix 618: :]] ( # # -- addr ) gsl_matrix_calloc ; 619 620: ]@ ( addr i -- df ) gsl_vector_get ; 621: ]! ( df addr i -- ) gsl_vector_set ; 622: ]data ( addr -- *data ) gsl_vector data @ ; 623: ]stride ( addr -- *data ) gsl_vector stride @ ; 624: ]fill ( df addr -- ) gsl_vector_set_all ; 625: ]erase ( addr -- ) gsl_vector_set_zero ; 626: ]+ ( *gsl_vector *gsl_vector -- ) gsl_vector_add drop ; 627: ]- ( *gsl_vector *gsl_vector -- ) gsl_vector_sub drop ; 628: ]e*! ( *gsl_vector *gsl_vector -- ) gsl_vector_mul drop ; 629: ]size ( *gsl_vector -- n ) gsl_vector size @ ; 630: ]outer* ( *gsl_vector *gsl_vector -- *gsl_matrix ) 631 over ]size over ]size gsl_matrix_calloc dup >r !1 632 gsl_blas_dger drop r> ; 633\ no control for divizion by zero (I get segfaults) 634: ]/ ( *gsl_vector *gsl_vector -- ) gsl_vector_div throw ; 635: ]clone ( *gsl_vector -- *gsl_vector ) 636 dup gsl_vector size @ gsl_vector_alloc 637 dup -rot swap gsl_vector_memcpy drop ; 638 639: ]add ( *gsl_vector *gsl_vector -- *gsl_vector ) 640 ]clone dup -rot swap ]+ ; 641: ]sub ( *gsl_vector *gsl_vector -- *gsl_vector ) 642 swap ]clone dup -rot swap ]- ; 643: ]mul ( *gsl_vector *gsl_vector -- *gsl_vector ) 644 swap ]clone dup -rot swap ]e*! ; 645: ]div ( *gsl_vector *gsl_vector -- *gsl_vector ) 646 swap ]clone dup -rot swap ]/ ; 647 648: ]*c ( df *gsl_vector -- ) gsl_vector_scale drop ; 649: ]+c ( df *gsl_vector -- ) gsl_vector_add_constant drop ; 650: ]max ( *gsl_vector -- ) gsl_vector_max ; 651: ]min ( *gsl_vector -- ) gsl_vector_min ; 652: ]imax ( *gsl_vector -- ) gsl_vector_max_index ; 653: ]imin ( *gsl_vector -- ) gsl_vector_min_index ; 654: ]copy] ( *gsl_vector_dest *gsl_vector_src -- ) gsl_vector_memcpy drop ; 655: ]negate !-1.0 ]*c ; 656: ]ones ( n -- *gsl_vector ) :] dup 1e ]fill ; 657 658: ]]slice ( *gsl_matrix x y n m -- *gsl_matrix ) 659 gsl_matrix_alloc_from_matrix ; 660: ]slice ( *gsl_vector offset length stride -- *gsl_vector ) 661 gsl_vector_alloc_from_vector ; 662 663: ]null? ( *gsl_vector -- 0/-1 ) gsl_vector_isnull negate ; 664: ]pos? ( *gsl_vector -- 0/-1 ) gsl_vector_ispos negate ; 665: ]neg? ( *gsl_vector -- 0/-1 ) gsl_vector_isneg negate ; 666 667\ FFT 19jan08sp 668: ]fft! ( *gsl_vector -- ) 669 dup ]size >r dup ]stride >r ]data r> r> 670 dup find-fft-cache if 671 dup gsl_fft_precomputes r_wavetable @ 672 swap gsl_fft_precomputes workspace @ 673 else 674 drop 675 dup cache-fft 676 dup gsl_fft_precomputes r_wavetable @ 677 swap gsl_fft_precomputes workspace @ 678 then 679 gsl_fft_real_transform throw ; 680: ]fft ( *gsl_vector -- *gsl_vector ) 681 ]clone dup ]fft! ; 682: ]ifft! ( *gsl_vector -- ) 683 dup ]size >r dup ]stride >r ]data r> r> 684 dup find-fft-cache if 685 dup gsl_fft_precomputes hc_wavetable @ 686 swap gsl_fft_precomputes workspace @ 687 else 688 drop 689 dup cache-fft 690 dup gsl_fft_precomputes hc_wavetable @ 691 swap gsl_fft_precomputes workspace @ 692 then 693 gsl_fft_hc_inverse throw ; 694: ]ifft ( *gsl_vector -- *gsl_vector ) 695 ]clone dup ]ifft! ; 696\ multiply two half complex vectors 697\ store result in the first 698: ]hc*! ( *gsl_vector *gsl_vector -- ) 699 2dup 0 ]@ 0 ]@ f* over 0 ]! 700 dup ]size dup %1 and not + 1 do 701 2dup 702 dup i ]@ i 1+ ]@ 703 dup i ]@ i 1+ ]@ z* 704 over dup i 1+ ]! i ]! 705 2 +loop 706 dup ]size %1 and not if 707 2dup 708 dup ]size 1- dup >r ]@ dup r@ ]@ f* r> ]! 709 then 710 2drop ; 711 712\ pseudomatrices and vectors 713: pvector ( *data n -- *gsl_vector ) 714 sizeof gsl_vector allocate throw 715 dup >r dup 1 swap gsl_vector stride ! 716 gsl_vector size ! r@ 717 gsl_vector data ! r@ 718 0 swap gsl_vector owner ! r> ; 719 720: pmatrix! ( *data tda n m *pmatrix -- *gsl_matrix ) 721 dup >r gsl_matrix size2 ! 722 r@ gsl_matrix size1 ! 723 r@ gsl_matrix tda ! 724 r@ gsl_matrix data ! 725 0 r@ gsl_matrix owner ! 726 r> ; 727 728: pmatrix ( *data tda n m -- *gsl_matrix ) 729 sizeof gsl_matrix allocate throw 730 dup 0 swap gsl_matrix owner ! 731 pmatrix! ; 732 733\ permutations 734 735: fpermutation ( n -- | -- id addr ) 736 create 737 gsl_permutation_calloc , 738 does> @ ; 739: }@ ( *gsl_permutation i -- n ) gsl_permutation_get ; 740: }data ( *gsl_permutation -- *data ) gsl_block data @ ; 741: }size ( *gsl_permutation -- *data ) gsl_block size @ ; 742: }free ( *gsl_permutation -- ) gsl_permutation_free ; 743: }sign ( *gsl_permutation -- 1/-1 ) 744 1 over dup }size 0 do 745 dup i }@ i <> if swap negate swap then 746 loop drop ; 747 748\ matrices 749 750: fmatrix ( n m -- | -- id addr ) 751 create 752 gsl_matrix_calloc , 753 GSL_MATRIX_TYPE , 754 does> @ ; 755 756: free_pseudomatrix ( pmatrix/pvector -- ) free throw ; 757 758create free_matrix ' free_pseudomatrix , ' gsl_matrix_free , 759create free_vector ' free_pseudomatrix , ' gsl_vector_free , 760 761: ]]free ( *gsl_matrix -- ) 762 dup gsl_matrix owner @ 763 cells free_matrix + @ execute ; 764: ]free ( addr -- ) 765 dup gsl_vector owner @ 766 cells free_vector + @ execute ; 767: ]]@ ( *gsl_matrix i j -- df ) gsl_matrix_get ; 768: ]]*@ ( *gsl_matrix i j -- *[i,j] ) gsl_matrix_ptr ; 769: ]]! ( *gsl_matrix i j df -- ) gsl_matrix_set ; 770: ]]fill ( addr df -- ) gsl_matrix_set_all ; 771: ]]size1 gsl_matrix size1 @ ; 772: ]]size2 gsl_matrix size2 @ ; 773: ]]dim ( *gsl_matrix -- m n ) dup ]]size1 swap ]]size2 ; 774: ]]dim. ( *gsl_matrix -- ) ]]dim swap . ." x" . cr ; 775: ]]data ( *gsl_matrix -- addr) gsl_matrix data @ ; 776: ]]tda gsl_matrix tda @ ; 777: ]]block gsl_matrix block @ ; 778: ]]owner gsl_matrix owner @ ; 779: ]]copy]] ( *gsl_matrix_dest *gsl_matrix_src -- ) gsl_matrix_memcpy drop ; 780: ]]'copy]] ( *gsl_matrix_dest *gsl_matrix_src -- ) gsl_matrix_transpose_memcpy drop ; 781\ : ]]row ( *gsl_matrix idx -- *gsl_vector ) gsl_matrix_row ; 782\ : ]]col ( *gsl_matrix idx -- *gsl_vector ) gsl_matrix_column ; 783: ]]>] ( *gsl_vector *gsl_matrix i -- ) gsl_matrix_get_col drop ; 784: ]]>]' ( *gsl_vector *gsl_matrix i -- ) gsl_matrix_get_row drop ; 785: ]>]] ( *gsl_matrix *gsl_vector i -- ) swap gsl_matrix_set_col drop ; 786: ]'>]] ( *gsl_matrix *gsl_vector i -- ) swap gsl_matrix_set_row drop ; 787 788: ]]max gsl_matrix_max ; 789: ]]min gsl_matrix_min ; 790: ]]add! ( *gsl_matrix *gsl_matrix -- ) 791 gsl_matrix_add drop ; 792: ]]sub! ( *gsl_matrix *gsl_matrix -- ) 793 gsl_matrix_sub drop ; 794: ]]e*! ( *gsl_matrix *gsl_matrix -- ) 795 gsl_matrix_mul_elements drop ; 796: ]]*c ( *gsl_matrix df -- ) 797 gsl_matrix_scale drop ; 798: ]]+c ( df *gsl_matrix -- ) gsl_matrix_add_constant drop ; 799: ]]clone ( *gsl_matrix -- *gsl_matrix ) 800 dup dup gsl_matrix size1 @ swap gsl_matrix size2 @ 801 gsl_matrix_alloc 802 dup -rot swap gsl_matrix_memcpy drop ; 803: ]]negate !-1.0 ]]*c ; 804 805: ]]+ ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) 806 ]]clone dup -rot swap ]]add! ; 807 808: ]]- ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) 809 swap ]]clone dup -rot swap ]]sub! ; 810: ]]null? ( *gsl_matrix -- 0/-1 ) gsl_matrix_isnull negate ; 811: ]]pos? ( *gsl_matrix -- 0/-1 ) gsl_matrix_ispos negate ; 812: ]]neg? ( *gsl_matrix -- 0/-1 ) gsl_matrix_isneg negate ; 813 814\ blas 815 816\ constants 817101 Constant CblasRowMajor 818102 Constant CblasColMajor 819111 Constant CblasNoTrans 820112 Constant CblasTrans 821113 Constant CblasConjTrans 822121 Constant CblasUpper 823122 Constant CblasLower 824131 Constant CblasNonUnit 825132 Constant CblasUnit 826141 Constant CblasLeft 827142 Constant CblasRight 828 829: action? ( *gsl_matrix *gsl_matrix n n n -- ) 830 dup 0= if 831 drop 832 2swap 2dup 833 ]]size2 swap ]]size1 swap 834 exit 835 then 836 dup 1 = if 837 drop 838 2swap 2dup 839 ]]size2 swap ]]size2 swap 840 exit 841 then 842 2 = if 843 2swap 2dup 844 ]]size1 swap ]]size1 swap 845 exit 846 then 847 3 = if 848 2swap 2dup 849 ]]size1 swap ]]size2 swap 850 exit 851 then ; 852 853create samemattable ' noop , ' ]]clone , 854: samemat ( *gsl_matrix *gsl_matrix -- 1/0 *gsl_matrix ) 855 dup -rot = abs dup -rot cells samemattable + @ execute ; macro 856 857: ]]mul ( *gsl_matrix *gsl_matrix n n n -- *gsl_matrix ) 858 !1 !0 action? 859 gsl_matrix_alloc dup >r 860 gsl_blas_dgemm drop r> ; 861: ]]* ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) 862 2dup samemat dup rot 2>r nip 863 CblasNoTrans CblasNoTrans 0 ]]mul 864 2r> if ]]free else drop then ; 865: ]]'* ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) 866 2dup samemat dup rot 2>r nip 867 CblasTrans CblasNoTrans 1 ]]mul 868 2r> if ]]free else drop then ; 869: ]]*' ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) 870 2dup samemat dup rot 2>r nip 871 CblasNoTrans CblasTrans 2 ]]mul 872 2r> if ]]free else drop then ; 873: ]]'*' ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) 874 2dup samemat dup rot 2>r nip 875 CblasTrans CblasTrans 3 ]]mul 876 2r> if ]]free else drop then ; 877 878: ]]mul! ( n n *gsl_matrix *gsl_matrix *gsl_matrix -- ) 879 !1 !0 gsl_blas_dgemm drop ; 880: ]]*! ( *gsl_matrix *gsl_matrix *gsl_matrix -- ) 881 >r CblasNoTrans CblasNoTrans 2swap r> ]]mul! ; 882: ]]'*! ( *gsl_matrix *gsl_matrix *gsl_matrix -- ) 883 >r CblasTrans CblasNoTrans 2swap r> ]]mul! ; 884: ]]*'! ( *gsl_matrix *gsl_matrix *gsl_matrix -- ) 885 >r CblasNoTrans CblasTrans 2swap r> ]]mul! ; 886 887 888: ]]*] ( *gsl_matrix *gsl_vector -- *gsl_vector ) 889 over ]]size1 gsl_vector_alloc >r 890 CblasNoTrans -rot r@ !1 !0 gsl_blas_dgemv drop r> ; 891: ]]'*] ( *gsl_matrix *gsl_vector -- *gsl_vector ) 892 over ]]size1 gsl_vector_alloc >r 893 CblasTrans -rot r@ !1 !0 gsl_blas_dgemv drop r> ; 894 895: ]]i ( *gsl_matrix -- ) 896 dup dup ]]size1 swap ]]size2 <> if 897 abort" ERROR: Not a square matrix!" 898 then 899 dup ]]size1 0 do 900 dup i i !1 ]]! 901 loop drop ; 902: identity ( n -- *gsl_matrix ) 903 dup gsl_matrix_calloc dup ]]i ; 904: min-identity ( *gsl_matrix -- *gsl_matrix ) 905 dup ]]size1 swap ]]size2 min identity ; 906: left/right' ( *gsl_matrix *gsl_matrix -- *gsl_matrix ) 907 over ]]size1 over ]]size1 > if 908 swap ]]*' exit 909 else 910 ]]'* exit 911 then ; 912 913\ original matrix remains intact 914: ]]' ( *gsl_matrix -- *gsl_matrix ) 915 dup min-identity dup >r 916 left/right' 917 r> ]]free ; 918: ]]T! ( *gsl_matrix -- ) 919 gsl_matrix_transpose drop ; 920 921: ]]T ( *gsl_matrix -- *gsl_matrix ) 922 dup ]]dim swap gsl_matrix_alloc dup rot gsl_matrix_transpose_memcpy drop ; 923 924: ]]2T ( *gsl_matr *gsl_matrix -- ) 925 gsl_matrix_transpose_memcpy drop ; 926 927: ]]+! ( *gsl_matrix i j df -- ) >r 2dup r@ ]]@ f+ r> ]]! ; 928: ]]scale! ( *gsl_matrix i j df -- ) >r 2dup r@ ]]@ f* r> ]]! ; 929: ]]data_ij ( *gsl_matrix i j -- addr) 930 rot >r swap r@ ]]tda dfloats * swap dfloats + r> ]]data + ; 931\ Cross product can be either calculated through determinant: 932: ]x ( *gsl_vector *gsl_vector -- *gsl_vector ) 933 3 gsl_vector_alloc 934 { x1[ x2[ x3[ | 935 x1[ 2 ]@ fnegate x2[ 1 ]@ f* x1[ 1 ]@ x2[ 2 ]@ f* f+ x3[ 0 ]! 936 x1[ 2 ]@ x2[ 0 ]@ f* x1[ 0 ]@ fnegate x2[ 2 ]@ f* f+ x3[ 1 ]! 937 x1[ 1 ]@ fnegate x2[ 0 ]@ f* x1[ 0 ]@ x2[ 1 ]@ f* f+ x3[ 2 ]! 938 x3[ } ; 939\ or using algebraic form when first vector in the product is 940\ rewritten in a matrix form: 941\ a x b = [C_a] b 942\ [ 0 -a[2] a[1] ] 943\ [C_a] = [ a[2] 0 -a[0] ] 944\ [-a[1] a[0] 0 ] 945\ a function to convert a vector into such matrix: 946: ]>[x] ( ] -- ]] ) 947 [IFDEF] debug 948 dup ]size 3 <> abort" Not a 3D vector!" 949 [THEN] 950 3 3 :]] dup >r 951 swap 2dup 2 ]@ fdup dup fnegate 0 1 ]]! 1 0 ]]! 952 2dup 1 ]@ fdup dup fnegate 2 0 ]]! 0 2 ]]! 953 0 ]@ fdup dup fnegate 1 2 ]]! 2 1 ]]! r> ; 954: ]. ( *gsl_vector *gsl_vector -- f:dot_product ) 955 { x1[ x2[ | 956 0 0 sp@ x1[ x2[ rot gsl_blas_ddot drop fd>f } ; 957: ]total ( *gsl_vector -- f:sum ) 958 dup ]size gsl_vector_alloc dup !1 ]fill dup rot ]. ]free ; 959\ probability normalize - assures sum is unity 960: ]pnormalize ( *gsl_vector - ) 961 dup ]total 1/f ]*c ; 962: |]| ( *gsl_vector -- f:norm ) dup ]. fsqrt ; 963\ assures vector norm is unity 964: ]normalize ( *gsl_vector - ) 965 dup |]| 1/f ]*c ; 966: ]distance ( *gsl-vector *gsl-vector -- f ) 967 ]sub dup |]| ]free ; 968: ]+! ( *gsl_vector i df -- ) 969 2dup ]@ f+ ]! ; 970: ]*! ( *gsl_vector i df -- ) 971 2dup ]@ f* ]! ; 972 973: ]]*]m ( *gsl_matrix *gsl_vector -- *gsl_vector ) 974 over ]]size1 gsl_vector_calloc 975 { m[[ x[ y[ | 976 m[[ ]]size1 0 do 977 m[[ ]]size2 0 do 978 m[[ j i ]]@ x[ i ]@ f* y[ j ]+! 979 loop 980 loop y[ } ; 981 982: >#rows ( -- ) 983 swap ]]size1 >= abort" number of rows is bigger than available!" ; 984: >#cols ( -- ) 985 swap ]]size2 >= abort" number of columns is bigger than available!" ; 986 987: ]]row ( *gsl_matrix n -- *gsl_vector ) 988 2dup >#rows 989 sizeof gsl_vector allocate throw 990 dup 1 swap gsl_vector stride ! >r 991 over ]]size2 r@ gsl_vector size ! 992 0 ]]data_ij r@ gsl_vector data ! 993 0 r@ gsl_vector owner ! r> ; 994\ assumes all dimensions are set correctly 995: ]]row! ( *gsl_vector *gsl_matrix n -- ) 996 rot >r 2dup >#rows 0 ]]data_ij r> gsl_vector data ! ; 997: ]]col ( *gsl_matrix n -- *gsl_vector ) 998 2dup >#cols 999 sizeof gsl_vector allocate throw >r 1000 over ]]tda r@ gsl_vector stride ! 1001 over ]]size1 r@ gsl_vector size ! 1002 over ]]block r@ gsl_vector block ! 1003 0 swap ]]data_ij r@ gsl_vector data ! 1004 0 r@ gsl_vector owner ! r> ; 1005: ]]rfill ( f:n *gsl_matrix i -- ) ]]row dup ]fill ]free ; 1006: ]]cfill ( f:n *gsl_matrix i -- ) ]]col dup ]fill ]free ; 1007 1008 1009: ]]submat ( *gsl_matrix n1 n2 m1 m2 -- *gsl_matrix ) 1010 { m[[ n1 n2 m1 m2 | 1011 sizeof gsl_matrix allocate throw >r 1012 n2 n1 - 1+ r@ gsl_matrix size1 ! 1013 m2 m1 - 1+ r@ gsl_matrix size2 ! 1014 m[[ n1 m1 ]]data_ij r@ gsl_matrix data ! 1015 m[[ ]]tda r@ gsl_matrix tda ! 1016 0 r@ gsl_matrix owner ! r> } ; 1017 1018: ?square ( *gsl_matrix -- ) 1019 dup ]]size1 swap ]]size2 <> abort" ERROR: Not a square matrix!" ; 1020: ]]diag ( *gsl_matrix n1 n2 -- *gsl_vector ) 1021 rot dup ?square -rot 1022 sizeof gsl_vector allocate throw { d[ | 1023 over - d[ gsl_vector size ! 1024 2dup dup ]]data_ij d[ gsl_vector data ! drop 1025 dup ]]tda d[ gsl_vector stride ! 1026 ]]block d[ gsl_vector block ! 1027 0 d[ gsl_vector owner ! 1028 d[ } ; 1029 1030\ with input matrix replaced by the result 1031: ]]gsl-svd ( *gsl_matrix -- *gsl_matrix *gsl_vector ) 1032 dup ]]size2 dup dup gsl_matrix_calloc 1033 swap dup gsl_vector_calloc swap 1034 gsl_vector_calloc 1035 { mV vS vW | 1036 mV vS vW gsl_linalg_SV_decomp drop 1037 vW ]free 1038 mV vS } ; 1039\ seems to be 30% faster 1040: ]]gsl-svdm ( *gsl_matrix -- *gsl_matrix *gsl_vector ) 1041 dup ]]size2 dup ( a n n -- ) 1042 dup dup gsl_matrix_calloc swap ( a n a n -- ) 1043 dup gsl_matrix_calloc rot dup ( a a a n n -- ) 1044 gsl_vector_calloc swap 1045 gsl_vector_calloc 1046 { mX mV vS vW | 1047 mX mV vS vW gsl_linalg_SV_decomp_mod drop 1048 vW ]free mX ]]free 1049 mV vS } ; 1050 1051 1052: ]]alu ( *gsl_matrix -- *gsl_permutation ) ( matrix replaced with its lu ) 1053 { a[[ | 1054 CblasRowMajor a[[ ]]size1 a[[ ]]size2 a[[ ]]data a[[ ]]size1 dup 1055 gsl_permutation_alloc dup >r }data 1056 clapack_dgetrf throw r> } ; 1057: ]]ainv ( *gsl_matrix *gsl_permutation -- ) 1058 \ LU of a matrix replaced with its inverse 1059 { a[[ t{ | 1060 CblasRowMajor a[[ ]]size2 a[[ ]]data a[[ ]]size1 t{ }data 1061 clapack_dgetri throw } ; 1062: ]]ainvert ( *gsl_matrix -- *gsl_matrix ) 1063 [IFDEF] отладка 1064 dup ?square 1065 [THEN] 1066 ]]clone dup dup >r ]]alu dup >r ]]ainv r> }free r> ; 1067: ]]det ( *gsl_matrix -- f:determinant ) 1068 [IFDEF] отладка 1069 dup ?square 1070 [THEN] 1071 ]]clone dup ]]alu >r 1e0 1072 dup ]]size1 0 do dup i dup ]]@ f* loop ]]free 1073 \ compute permutation sign 1074 r> }sign s>f f* }free ; 1075\ calculates the work needed for dgesvd_ ( see man dgesvd ) 1076: lwork ( m n -- c ) 1077 2dup max -rot min 3 * over + swap 5 * max ; 1078\ this svd returns U MxM so eats a lot of memory 1079: ]]asvda ( *gsl_matrix -- *gsl_matrix *gsl_matrix *gsl_vector ) 1080 ]]clone { A[[ | 1081 A[[ ]]size1 dup gsl_matrix_alloc 1082 A[[ ]]size2 dup gsl_matrix_alloc 1083 A[[ ]]size1 A[[ ]]size2 min gsl_vector_alloc 1084 8 cells allocate throw 1085 { U[[ V[[ W[ p[ | 1086 ascii A p[ 0 cells + ! p[ 0 cells + 1087 ascii A p[ 1 cells + ! p[ 1 cells + 1088 A[[ ]]size1 p[ 2 cells + ! p[ 2 cells + 1089 A[[ ]]size2 p[ 3 cells + ! p[ 3 cells + 1090 A[[ ]]data 1091 p[ 2 cells + 1092 W[ ]data 1093 U[[ ]]data 1094 U[[ ]]size1 p[ 4 cells + ! p[ 4 cells + 1095 V[[ ]]data 1096 V[[ ]]size1 p[ 5 cells + ! p[ 5 cells + 1097 A[[ ]]size1 A[[ ]]size2 lwork 1098 dup gsl_vector_alloc dup >r 1099 ]data swap p[ 6 cells + ! p[ 6 cells + 1100 p[ 7 cells + 1101 dgesvd_ 1102 r> ]free p[ free throw A[[ ]]free 1103 U[[ V[[ W[ } } ; 1104 1105\ performs A=U*S*V^T 1106\ A = MxN, where M>N, pass it A^T 1107\ returns U^T (MxN), V(NxN) and vector of N eigenvalues 1108: ]]asvdO ( *gsl_matrix -- *gsl_matrix *gsl_matrix *gsl_vector ) 1109 { A[[ | 1110 A[[ ]]size2 A[[ ]]size1 min dup gsl_matrix_alloc 1111 A[[ ]]size1 A[[ ]]size2 min gsl_vector_alloc 1112 8 cells allocate throw 1113 { V[[ W[ p[ | 1114 ascii O p[ 0 cells + ! p[ 0 cells + 1115 ascii S p[ 1 cells + ! p[ 1 cells + 1116 A[[ ]]size2 p[ 2 cells + ! p[ 2 cells + 1117 A[[ ]]size1 p[ 3 cells + ! p[ 3 cells + 1118 A[[ ]]data 1119 p[ 2 cells + 1120 W[ ]data 1121 0 1122 p[ 2 cells + 1123 V[[ ]]data 1124 V[[ ]]size2 p[ 5 cells + ! p[ 5 cells + 1125 A[[ ]]size2 A[[ ]]size1 lwork 1126 dup gsl_vector_alloc dup >r 1127 ]data swap p[ 6 cells + ! p[ 6 cells + 1128 p[ 7 cells + 1129 dgesvd_ 1130 r> ]free p[ free throw 1131 A[[ V[[ W[ } } ; 1132 1133 1134: ]diag[[ ( *gsl_vector -- *gsl_matrix ) 1135 dup ]size dup dup gsl_matrix_calloc swap 1136 0 do 1137 2dup swap i ]@ i i ]]! 1138 loop nip ; 1139 1140: ]print\ ( *gsl_vector -- ) 1141 dup ]size 0 do dup i ]@ fx. loop drop ; 1142: ]print ( *gsl_vector -- ) ]print\ cr ; 1143: ]]print ( *gsl_matrix -- ) 1144 cr precision swap 1145 5 set-precision 1146 dup ]]size1 0 do 1147 \ i . ." : " 1148 dup ]]size2 0 do 1149 dup 1150 j i ]]@ fs. 1151 loop 1152 cr 1153 loop 1154 drop set-precision ; 1155: ]]row-print ( *gsl_matrix i -- ) 1156 cr 1157 over gsl_matrix size2 @ 0 do 1158 2dup 1159 i ]]@ f. 1160 loop 1161 cr 2drop ; 1162 1163: ]]col-print ( *gsl_matrix i -- ) 1164 cr 1165 over gsl_matrix size1 @ 0 do 1166 2dup 1167 i swap ]]@ f. 1168 loop 1169 cr 2drop ; 1170 1171: ]]nthrow ( *gsl_matrix n -- addr ) 1172 over ]]tda * dfloats swap ]]data + ; 1173 1174: ]]randomize ( *gsl_matrix -- ) 1175 dup dup ]]size1 swap ]]size2 * 0 do 1176 dup 1177 gsl-randomu 1178 ]]data i dfloats + df! 1179 loop drop ; 1180: ]randomize ( *gsl_vector -- ) 1181 dup ]size 0 do 1182 dup 1183 gsl-randomu 1184 i ]! 1185 loop drop ; 1186: ]mean ( *gsl_vector -- f ) 1187 dup ]stride swap dup ]size swap ]data 1188 rot rot gsl_stats_mean ; 1189 1190: ]variance ( *gsl_vector -- f ) 1191 dup ]stride swap dup ]size swap ]data 1192 rot rot gsl_stats_variance ; 1193 1194: ]sd ( *gsl_vector -- f ) 1195 dup ]stride swap dup ]size swap ]data 1196 rot rot gsl_stats_sd ; 1197 1198: ]skew ( *gsl_vector -- f ) 1199 dup ]stride swap dup ]size swap ]data 1200 rot rot gsl_stats_skew ; 1201 1202: ]kurtosis ( *gsl_vector -- f ) 1203 dup ]stride swap dup ]size swap ]data 1204 rot rot gsl_stats_kurtosis ; 1205 1206: ]]gsl-lu ( *gsl_matrix -- *gsl_matrix *gsl_permutation ) 1207 1 sp@ rot ]]clone dup >r dup ]]size1 gsl_permutation_calloc dup >r rot 1208 gsl_linalg_LU_decomp drop r> r> swap rot drop ; 1209 1210: ]]gsl-invert ( *gsl_matrix -- *gsl_matrix ) 1211 ]]clone dup dup ]]gsl-lu 2dup >r >r rot 1212 gsl_linalg_LU_invert drop r> ]]free r> }free ; 1213 1214' ]]ainvert alias ]]invert 1215' ]]asvdO alias ]]svd 1216 1217: ]]save ( *gsl_matrix *gsl_matrix_cfa fid -- ) 1218 -rot { m[[ name[[ | 1219 >r 1220 name[[ >name count 1+ nip 0 m[[ ]]size2 m[[ ]]size1 0 1221 sp@ 5 cells r@ write-file throw 1222 2drop 2drop drop 1223 name[[ >name count 1+ r@ write-file throw 1224 m[[ ]]size1 m[[ ]]size2 * dfloats m[[ ]]T dup s>f ]]data swap 1225 r> write-file throw f>s ]]free } ; 1226 1227\ these words do not work with float matrices but are needed for 1228\ scientific calculations, that's why they are in this module 1229 1230: _hmatrix ( n m size -- addr ) 1231 rot over * 2 pick * [ 2 cells ] literal + 1232 allocate throw dup [ 2 cells ] literal + >r 1233 rot over ! [ 1 cells ] literal + ! r> ; 1234: hmatrix ( n m size -- ) 1235 create 1236 rot over * 2 pick * [ 2 cells ] literal + allocate throw dup , 1237 rot over ! [ 1 cells ] literal + ! 1238 does> @ [ 2 cells ] literal + ; 1239: }}row-size ( hmatrix -- ) [ 2 cells ] literal - @ ; 1240: freeHmatrix ( hmatrix -- ) [ 2 cells ] literal - free throw ; 1241: }} ( addr i j -- addr[i][j] ) \ word to fetch 2-D array addresses 1242 >R >R \ indices to return stack temporarily 1243 DUP CELL- CELL- 2@ \ &a[0][0] size m 1244 R> * R> + * 1245 + 1246 ALIGNED ; 1247: h->[[ ( hmatrix -- gsl_matrix ) 1248 dup }}row-size 3 swap gsl_matrix_alloc 1249 dup ]]size2 0 do 1250 3 0 do 1251 2dup swap i j }} w@ s>f i j ]]! 1252 loop 1253 loop nip ; 1254 1255\ some sequencing code 1256: arange ( f:start f:end f:step -- x[ ) 1257 f-rot 1258 fswap fdup f>r f- fover f/ f>s :] fr> 1259 dup ]size 0 do 1260 dup fover i s>f f* fover f+ i ]! 1261 loop ; 1262: product ( x[ -- f:P ) 1263 !1 dup ]size 0 do 1264 dup i ]@ f* 1265 loop drop ; 1266 1267\ initializing random number generator to some value in order to have 1268\ it available upon loading of gsl 1269mt19937 1270: )randperm ( *v( -- ) 1271 gsl_rng_default swap 1272 dup )size over )type 8 / gsl_ran_shuffle ; 1273 1274previous previous previous previous previous 1275 1276Module; 1277