gsl.fs 57.1 KB
Newer Older
bp's avatar
bp committed
1 2 3 4 5 6 7 8 9 10 11 12 13
\ 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.

bp's avatar
bp committed
14
\needs float    import float
bp's avatar
bp committed
15 16 17
warning off
\needs locals|  include locals.fs
warning on
bp's avatar
bp committed
18
\needs atlas    include atlas.fs
bp's avatar
bp committed
19
\needs callback include callback.fs
bp's avatar
bp committed
20 21
\needs vectors  include vectors.fs
\needs complex  include complex.fb
bp's avatar
bp committed
22 23 24

Module GSL

bp's avatar
bp committed
25 26 27
\ stole the hash function from hash.fs of bigforth
| &14 Value Hashbits
| 1 Hashbits << Value Hashlen
bp's avatar
bp committed
28

bp's avatar
bp committed
29 30 31 32 33 34 35 36 37 38
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
bp's avatar
bp committed
39

bp's avatar
bp committed
40 41 42 43 44 45 46 47 48 49
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
bp's avatar
bp committed
50 51 52

legacy off

bp's avatar
bp committed
53 54 55 56 57 58 59 60 61 62 63 64 65
\ 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
bp's avatar
bp committed
66
            1+
bp's avatar
bp committed
67 68
    repeat ;

bp's avatar
bp committed
69 70 71 72
| : .bold-red ." " ;
| : .red ." " ;
| : .reset    ." " ;
| : cb-test
bp's avatar
bp committed
73
    cr
bp's avatar
bp committed
74 75 76
    \ .bold-red
    ." GSL ERROR: " cr
    \ .reset cr
bp's avatar
bp committed
77 78 79 80
    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
bp's avatar
bp committed
81 82
    \ .red
    -1 abort" failed at" ;
bp's avatar
bp committed
83 84 85
' cb-test 4:0 c_plus

\ 1 2 c_plus 2:1call .
bp's avatar
bp committed
86
| variable old_handler
bp's avatar
bp committed
87 88
c_plus gsl_set_error_handler old_handler !

bp's avatar
bp committed
89 90
0 Constant GSL_SUCCESS

bp's avatar
bp committed
91 92 93 94 95 96 97 98 99 100 101 102
\ 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

bp's avatar
bp committed
103 104 105 106 107 108 109 110 111
\ 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 )
bp's avatar
bp committed
112
libgsl gsl_rng_uniform_int int int (int) gsl_rng_uniform_int ( *gsl_rng n --n )
bp's avatar
bp committed
113 114 115 116
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 )
bp's avatar
bp committed
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181
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
bp's avatar
bp committed
182 183 184 185 186 187 188 189 190 191 192


\ 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 )
bp's avatar
bp committed
193 194 195
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
bp's avatar
bp committed
196 197 198 199
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

bp's avatar
bp committed
200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228
\ 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

bp's avatar
bp committed
229 230
\ Vectors 
libgsl gsl_block_alloc  int (int) gsl_block_alloc ( n -- addr )
bp's avatar
bp committed
231
libgsl gsl_block_calloc int (int) gsl_block_calloc ( n -- addr )
bp's avatar
bp committed
232
libgsl gsl_block_free   int (int) gsl_block_free ( n -- addr )
bp's avatar
bp committed
233

bp's avatar
bp committed
234 235 236 237 238
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 )
bp's avatar
bp committed
239
libgsl gsl_vector_set df int int (void/fp) gsl_vector_set ( df addr i --  )
bp's avatar
bp committed
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
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
bp's avatar
bp committed
256 257

libgsl gsl_vector_subvector int int int (int) gsl_vector_subvector
bp's avatar
bp committed
258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301
\ 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 )
bp's avatar
bp committed
302 303 304 305 306
\ 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
bp's avatar
bp committed
307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337
\ 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 )
bp's avatar
bp committed
338

bp's avatar
bp committed
339 340 341 342 343

\ 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
344
int df int int int int (void/fp) cblas_dgemv
bp's avatar
bp committed
345 346
libgsl gsl_blas_ddot int int int (int) gsl_blas_ddot
( *gsl_vector *gsl_vector df -- n )
bp's avatar
bp committed
347
libgsl gsl_blas_dgemm int df int int df int int (int/fp) gsl_blas_dgemm
bp's avatar
bp committed
348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363
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 )

bp's avatar
bp committed
364
\ -----------------------------------------------------------------------------
bp's avatar
bp committed
365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392
\                  *** 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
( -- )
bp's avatar
bp committed
393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414
\ --- 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 @ ;
bp's avatar
bp committed
415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431

\ --- 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 --  )
bp's avatar
bp committed
432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466
\ -----------------------------------------------------------------------------
\                     *** 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
bp's avatar
bp committed
467
| : 2str dup >r abs s>d <<# #s r> sign #> $buf 0place #>> ;
bp's avatar
bp committed
468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497
| : 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 ;
bp's avatar
bp committed
498

bp's avatar
bp committed
499
legacy on
bp's avatar
bp committed
500

bp's avatar
bp committed
501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547
\ 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
bp's avatar
bp committed
548
\ gsl_rng_types_setup  value gsl_rng_array(
bp's avatar
bp committed
549 550 551 552
0 value gsl_rng_default
: gsl-free ( -- )
    gsl_rng_default gsl_rng_free ;

bp's avatar
bp committed
553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680
\ 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 ;
bp's avatar
bp committed
681 682 683 684 685 686 687 688 689 690 691 692 693 694 695

\ 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 )
bp's avatar
bp committed
696
    gsl_rng_default !1 gsl_ran_gaussian_ratio_method ;
bp's avatar
bp committed
697 698
: gsl-discrete ( *gsl_ran_discrete -- n )
    gsl_rng_default swap gsl_ran_discrete ;
bp's avatar
bp committed
699 700

\ vectors and matrices
bp's avatar
bp committed
701 702 703
0 constant GSL_VECTOR_TYPE
1 constant GSL_MATRIX_TYPE
: gsltypeof ' >body cell + @ ;
bp's avatar
bp committed
704 705

: fvector ( n -- | -- id addr )
bp's avatar
bp committed
706
    create    
bp's avatar
bp committed
707
    gsl_vector_calloc ,
bp's avatar
bp committed
708
    GSL_VECTOR_TYPE ,
bp's avatar
bp committed
709 710
  does> @ ;

bp's avatar
bp committed
711 712 713 714 715
\ allocate a nameless vector
: :] ( # -- addr ) gsl_vector_calloc ;
\ allocate a nameless matrix
: :]] ( # # -- addr ) gsl_matrix_calloc ;

bp's avatar
bp committed
716 717 718
: ]@ ( addr i -- df ) gsl_vector_get ;
: ]! ( df addr i -- ) gsl_vector_set ;
: ]data ( addr -- *data ) gsl_vector data @ ;
bp's avatar
bp committed
719
: ]stride ( addr -- *data ) gsl_vector stride @ ;
bp's avatar
bp committed
720 721 722 723 724 725 726 727 728
: ]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> ;
bp's avatar
bp committed
729
\ no control for division by zero in gsl (is it NAN now?)
bp's avatar
bp committed
730 731 732 733
: ]/ ( *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 ;
bp's avatar
bp committed
734
' gsl_sort_vector alias ]sort ( v[ -- )
bp's avatar
bp committed
735 736 737 738 739 740 741 742 743 744

: ]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 ]/ ;

bp's avatar
bp committed
745
\ "*, /, +, -" constant to a vector
bp's avatar
bp committed
746
: ]*c ( df *gsl_vector -- ) gsl_vector_scale drop ;
bp's avatar
bp committed
747 748 749 750 751 752 753 754 755 756 757 758 759 760
: ]/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 ;

bp's avatar
bp committed
761 762
: ]copy] ( *gsl_vector_dest *gsl_vector_src -- ) gsl_vector_memcpy drop ;
: ]negate !-1.0 ]*c ;
bp's avatar
bp committed
763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837
: ]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
bp's avatar
bp committed
838 839 840 841
: pvector ( *data n -- *gsl_vector )
    sizeof gsl_vector allocate throw
    dup >r dup 1 swap gsl_vector stride !
    gsl_vector size ! r@
bp's avatar
bp committed
842 843
    gsl_vector data ! r@
    0 swap gsl_vector owner ! r> ;
bp's avatar
bp committed
844

bp's avatar
bp committed
845 846 847 848 849
: 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 !
bp's avatar
bp committed
850
    0 r@ gsl_matrix owner !
bp's avatar
bp committed
851 852 853 854
    r> ;

: pmatrix ( *data tda n m -- *gsl_matrix )
    sizeof gsl_matrix allocate throw
bp's avatar
bp committed
855
    dup 0 swap gsl_matrix owner !
bp's avatar
bp committed
856 857 858 859 860 861 862 863
    pmatrix! ;

\ permutations

: fpermutation ( n -- | -- id addr )
    create
    gsl_permutation_calloc ,
  does> @ ;
bp's avatar
bp committed
864 865
: :}    ( n -- *gsl_permutation ) gsl_permutation_calloc ;
: }@    ( *gsl_permutation i -- n ) gsl_permutation_get ;
bp's avatar
bp committed
866
: }data ( *gsl_permutation -- *data ) gsl_block data @ ;
bp's avatar
bp committed
867
: }size ( *gsl_permutation -- *data ) gsl_block size @ ;
bp's avatar
bp committed
868
: }free ( *gsl_permutation -- ) gsl_permutation_free ;
bp's avatar
bp committed
869 870 871 872 873 874 875 876 877 878 879 880
: }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> ;
bp's avatar
bp committed
881

bp's avatar
bp committed
882
\ matrices
bp's avatar
bp committed
883

bp's avatar
bp committed
884 885 886
: fmatrix ( n m -- | -- id addr )
    create
    gsl_matrix_calloc ,
bp's avatar
bp committed
887
    GSL_MATRIX_TYPE ,
bp's avatar
bp committed
888 889
  does> @ ;

bp's avatar
bp committed
890 891
: free_pseudomatrix ( pmatrix/pvector -- ) free throw ;

bp's avatar
bp committed
892 893 894 895 896 897 898 899 900 901 902 903 904
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 ;
bp's avatar
bp committed
905
: ]]erase ( addr -- ) gsl_matrix_set_zero ;
bp's avatar
bp committed
906 907 908
: ]]size1 gsl_matrix size1 @ ;
: ]]size2 gsl_matrix size2 @ ;
: ]]dim ( *gsl_matrix -- m n ) dup ]]size1 swap ]]size2 ;
bp's avatar
bp committed
909
: ]]dim. ( *gsl_matrix -- ) ]]dim swap . ." x" . cr ;
bp's avatar
bp committed
910 911 912
: ]]data ( *gsl_matrix -- addr) gsl_matrix data @ ;
: ]]tda gsl_matrix tda @ ;
: ]]block gsl_matrix block @ ;
bp's avatar
bp committed
913
: ]]owner gsl_matrix owner @ ;
bp's avatar
bp committed
914
: ]]copy]] ( *gsl_matrix_dest *gsl_matrix_src -- ) gsl_matrix_memcpy drop ;
bp's avatar
bp committed
915
: ]]'copy]] ( *gsl_matrix_dest *gsl_matrix_src -- ) gsl_matrix_transpose_memcpy drop ;
bp's avatar
bp committed
916 917
\ : ]]row ( *gsl_matrix idx -- *gsl_vector ) gsl_matrix_row ;
\ : ]]col ( *gsl_matrix idx -- *gsl_vector ) gsl_matrix_column ;
bp's avatar
bp committed
918 919 920 921
: ]]>]  ( *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 ;
bp's avatar
bp committed
922 923 924 925 926 927 928 929 930 931 932

: ]]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 ;
bp's avatar
bp committed
933
: ]]+c ( df *gsl_matrix -- ) gsl_matrix_add_constant drop ;
bp's avatar
bp committed
934 935 936 937 938 939 940 941 942 943 944
: ]]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! ;
bp's avatar
bp committed
945 946 947
: ]]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 ;
bp's avatar
bp committed
948

bp's avatar
bp committed
949 950 951 952 953
\ blas

\ constants
101 Constant CblasRowMajor
102 Constant CblasColMajor
bp's avatar
bp committed
954
111 Constant CblasNoTrans
bp's avatar
bp committed
955 956 957 958 959 960 961 962
112 Constant CblasTrans
113 Constant CblasConjTrans
121 Constant CblasUpper
122 Constant CblasLower
131 Constant CblasNonUnit
132 Constant CblasUnit
141 Constant CblasLeft
142 Constant CblasRight
bp's avatar
bp committed
963

bp's avatar
bp committed
964 965
: action? (  *gsl_matrix *gsl_matrix n n n -- )
    dup 0= if
bp's avatar
bp committed
966 967 968 969
        drop
        2swap 2dup
        ]]size2 swap ]]size1 swap
        exit
bp's avatar
bp committed
970 971
    then
    dup 1 = if
bp's avatar
bp committed
972 973 974 975
        drop
        2swap 2dup
        ]]size2 swap ]]size2 swap
        exit
bp's avatar
bp committed
976 977
    then
    2 = if
bp's avatar
bp committed
978 979 980
        2swap 2dup
        ]]size1 swap ]]size1 swap
        exit
bp's avatar
bp committed
981 982
    then
    3 = if
bp's avatar
bp committed
983 984 985
        2swap 2dup
        ]]size1 swap ]]size2 swap
        exit
bp's avatar
bp committed
986 987 988 989 990 991
    then ;

create samemattable ' noop , ' ]]clone ,
: samemat (  *gsl_matrix *gsl_matrix -- 1/0 *gsl_matrix )
    dup -rot = abs dup -rot cells samemattable + @ execute ; macro

bp's avatar
bp committed
992 993
: ]]mul (  *gsl_matrix *gsl_matrix n n n -- *gsl_matrix )
    !1 !0 action?
bp's avatar
bp committed
994 995
    gsl_matrix_alloc dup >r
    gsl_blas_dgemm drop r> ;
bp's avatar
bp committed
996
: ]]* (  *gsl_matrix *gsl_matrix -- *gsl_matrix )
bp's avatar
bp committed
997 998 999
    2dup samemat dup rot 2>r nip
    CblasNoTrans CblasNoTrans 0 ]]mul
    2r> if ]]free else drop then ;
bp's avatar
bp committed
1000
: ]]'* (  *gsl_matrix *gsl_matrix -- *gsl_matrix )
bp's avatar
bp committed
1001 1002 1003
    2dup samemat dup rot 2>r nip
    CblasTrans CblasNoTrans 1 ]]mul
    2r> if ]]free else drop then ;    
bp's avatar
bp committed
1004
: ]]*' (  *gsl_matrix *gsl_matrix -- *gsl_matrix )
bp's avatar
bp committed
1005 1006 1007
    2dup samemat dup rot 2>r nip
    CblasNoTrans CblasTrans 2 ]]mul
    2r> if ]]free else drop then ;        
bp's avatar
bp committed
1008
: ]]'*' (  *gsl_matrix *gsl_matrix -- *gsl_matrix )
bp's avatar
bp committed
1009 1010 1011
    2dup samemat dup rot 2>r nip
    CblasTrans CblasTrans 3 ]]mul
    2r> if ]]free else drop then ;        
bp's avatar
bp committed
1012

bp's avatar
bp committed
1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023
: ]]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 )
bp's avatar
bp committed
1024
    over ]]size1 gsl_vector_alloc >r
bp's avatar
bp committed
1025 1026
    CblasNoTrans -rot r@ !1 !0 gsl_blas_dgemv drop r> ;
: ]]'*] ( *gsl_matrix *gsl_vector -- *gsl_vector )
bp's avatar
bp committed
1027
    over ]]size1 gsl_vector_alloc >r
bp's avatar
bp committed
1028 1029 1030 1031
    CblasTrans -rot r@ !1 !0 gsl_blas_dgemv drop r> ;

: ]]i ( *gsl_matrix -- )
    dup dup ]]size1 swap ]]size2 <> if
bp's avatar
bp committed
1032
        abort" ERROR: Not a square matrix!"
bp's avatar
bp committed
1033 1034
    then
    dup ]]size1 0 do
bp's avatar
bp committed
1035
        dup i i !1 ]]! 
bp's avatar
bp committed
1036 1037 1038 1039 1040 1041 1042
    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
bp's avatar
bp committed
1043
        swap ]]*' exit
bp's avatar
bp committed
1044
    else
bp's avatar
bp committed
1045
        ]]'* exit
bp's avatar
bp committed
1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065
    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 + ;
bp's avatar
bp committed
1066
\ Cross product can be either calculated through determinant:
bp's avatar
bp committed
1067 1068 1069 1070 1071 1072 1073
: ]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[ } ;
bp's avatar
bp committed
1074 1075 1076 1077 1078 1079 1080 1081 1082
\ 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
bp's avatar
bp committed
1083
        dup ]size 3 <> abort" Not a 3D vector!"
bp's avatar
bp committed
1084 1085 1086 1087 1088
    [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> ;
bp's avatar
bp committed
1089 1090 1091 1092
: ]. ( *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 )
bp's avatar
bp committed
1093
    dup ]size gsl_vector_alloc dup !1 ]fill dup rot ]. ]free ;
bp's avatar
bp committed
1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111
\ 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
bp's avatar
bp committed
1112 1113 1114
        m[[ ]]size2 0 do
            m[[ j i ]]@ x[ i ]@ f* y[ j ]+! 
        loop
bp's avatar
bp committed
1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128
    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> ;
bp's avatar
bp committed
1129 1130 1131
\ assumes all dimensions are set correctly
: ]]row! ( *gsl_vector *gsl_matrix n -- )
    rot >r 2dup >#rows 0 ]]data_ij r> gsl_vector data ! ;
bp's avatar
bp committed
1132 1133 1134 1135 1136 1137 1138 1139
: ]]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> ;
bp's avatar
bp committed
1140 1141 1142 1143
: ]]rfill ( f:n *gsl_matrix i -- ) ]]row dup ]fill ]free ;
: ]]cfill ( f:n *gsl_matrix i -- ) ]]col dup ]fill ]free ;


bp's avatar
bp committed
1144 1145 1146
: ]]submat ( *gsl_matrix n1 n2 m1 m2 -- *gsl_matrix )
    { m[[ n1 n2 m1 m2 |
    sizeof gsl_matrix allocate throw >r
bp's avatar
bp committed
1147 1148 1149 1150 1151
    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> } ;
bp's avatar
bp committed
1152
    
bp's avatar
bp committed
1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163
: ?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[ } ;
bp's avatar
bp committed
1164

bp's avatar
bp committed
1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185
\ 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 } ;

bp's avatar
bp committed
1186

bp's avatar
bp committed
1187
: ]]alu ( *gsl_matrix -- *gsl_permutation ) ( matrix replaced with its lu )
bp's avatar
bp committed
1188
    { a[[ | 
bp's avatar
bp committed
1189 1190 1191
    CblasRowMajor a[[ ]]size1 a[[ ]]size2 a[[ ]]data a[[ ]]size1 dup
    gsl_permutation_alloc dup >r }data
    clapack_dgetrf throw r> } ;
bp's avatar
bp committed
1192 1193 1194 1195 1196 1197 1198 1199
\ 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 } ;
    
bp's avatar
bp committed
1200
: ]]ainv ( *gsl_matrix *gsl_permutation -- )
bp's avatar
bp committed
1201
    \ LU of a matrix replaced with its inverse 
bp's avatar
bp committed
1202 1203 1204 1205
    { a[[ t{ |
    CblasRowMajor a[[ ]]size2 a[[ ]]data a[[ ]]size1 t{ }data
    clapack_dgetri throw } ;
: ]]ainvert ( *gsl_matrix -- *gsl_matrix )
bp's avatar
bp committed
1206 1207 1208
    [IFDEF] отладка
        dup ?square
    [THEN]
bp's avatar
bp committed
1209
    ]]clone dup dup >r ]]alu dup >r ]]ainv r> }free r> ;
bp's avatar
bp committed
1210 1211 1212 1213 1214 1215 1216 1217 1218
: ]]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 )
bp's avatar
bp committed
1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243
: 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 +
bp's avatar
bp committed
1244
    dgesvd
bp's avatar
bp committed
1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271
    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 +
bp's avatar
bp committed
1272
    dgesvd
bp's avatar
bp committed
1273 1274 1275
    r> ]free p[ free throw
    A[[ V[[ W[ } } ;

bp's avatar
bp committed
1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295
\ 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 ;

bp's avatar
bp committed
1296 1297 1298 1299

: ]diag[[ ( *gsl_vector -- *gsl_matrix )
    dup ]size dup dup gsl_matrix_calloc swap
    0 do
bp's avatar
bp committed
1300
        2dup swap i ]@ i i ]]!
bp's avatar
bp committed
1301 1302
    loop nip ;

bp's avatar
bp committed
1303 1304
: ]print\ ( *gsl_vector -- )
    dup ]size 0 do dup i ]@ fx. loop drop ;
bp's avatar
bp committed
1305
: ]print ( *gsl_vector -- )
bp's avatar
bp committed
1306 1307 1308
    precision swap
    5 set-precision
    ]print\ cr set-precision ;
bp's avatar
bp committed
1309
: ]]print ( *gsl_matrix -- )
bp's avatar
bp committed
1310 1311
    cr precision swap
    5 set-precision
bp's avatar
bp committed
1312
    dup ]]size1 0 do
bp's avatar
bp committed
1313 1314 1315
        \ i . ." :  "
        dup ]]size2 0 do
            dup
bp's avatar
bp committed
1316
            j i ]]@ fs.
bp's avatar
bp committed
1317 1318
        loop
        cr
bp's avatar
bp committed
1319
    loop
bp's avatar
bp committed
1320 1321
    drop set-precision ;

bp's avatar
bp committed
1322 1323 1324
: ]]row-print ( *gsl_matrix i -- )
    cr
    over gsl_matrix size2 @ 0 do
bp's avatar
bp committed
1325 1326
        2dup
         i ]]@ f.
bp's avatar
bp committed
1327 1328 1329 1330 1331 1332
    loop
    cr 2drop ;

: ]]col-print ( *gsl_matrix i -- )
    cr
    over gsl_matrix size1 @ 0 do
bp's avatar
bp committed
1333 1334
        2dup
        i swap ]]@ f.
bp's avatar
bp committed
1335 1336 1337 1338 1339 1340 1341 1342
    loop
    cr 2drop ;

: ]]nthrow ( *gsl_matrix n -- addr )
    over ]]tda * dfloats swap ]]data + ;

: ]]randomize ( *gsl_matrix -- )
    dup dup ]]size1 swap ]]size2 * 0 do
bp's avatar
bp committed
1343 1344 1345
            dup
            gsl-randomu
            ]]data i dfloats + df!
bp's avatar
bp committed
1346
    loop drop ;
bp's avatar
bp committed
1347 1348 1349 1350 1351 1352 1353 1354
: ]]randn ( *gsl_matrix -- )
    dup dup ]]size1 swap ]]size2 * 0 do
            dup
            gsl-gaussian
            ]]data i dfloats + df!
    loop drop ;
' ]]randomize alias ]]randu

bp's avatar
bp committed
1355 1356
: ]randomize ( *gsl_vector -- )
    dup ]size 0 do
bp's avatar
bp committed
1357 1358 1359
            dup
            gsl-randomu
            i ]!
bp's avatar
bp committed
1360
    loop drop ;
bp's avatar
bp committed
1361 1362 1363 1364
: ]randn ( *gsl_vector -- ) \ zero mean unit variance
    dup ]size 0 do dup gsl-gaussian i ]! loop drop ;
' ]randomize alias ]randu

bp's avatar
bp committed
1365
: ]mean ( *gsl_vector -- f )
bp's avatar
bp committed
1366 1367
    dup ]stride swap dup ]size swap ]data
    rot rot gsl_stats_mean ;
bp's avatar
bp committed
1368 1369

: ]variance ( *gsl_vector -- f )
bp's avatar
bp committed
1370 1371
    dup ]stride swap dup ]size swap ]data
    rot rot gsl_stats_variance ;
bp's avatar
bp committed
1372 1373

: ]sd ( *gsl_vector -- f )
bp's avatar
bp committed
1374 1375
    dup ]stride swap dup ]size swap ]data    
    rot rot gsl_stats_sd ;
bp's avatar
bp committed
1376 1377

: ]skew ( *gsl_vector -- f )
bp's avatar
bp committed
1378 1379
    dup ]stride swap dup ]size swap ]data    
    rot rot gsl_stats_skew ;
bp's avatar
bp committed
1380 1381

: ]kurtosis ( *gsl_vector -- f )
bp's avatar
bp committed
1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393
    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 ;
bp's avatar
bp committed
1394

bp's avatar
bp committed
1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435
: ]]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 ;
bp's avatar
bp committed
1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454
: 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 ;
bp's avatar
bp committed
1455

bp's avatar
bp committed
1456 1457 1458
\ initializing random number generator to some value in order to have
\ it available upon loading of gsl
mt19937
bp's avatar
bp committed
1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471
\ 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
bp's avatar
bp committed
1472 1473

Module;