diff --git a/package/Project.toml b/package/Project.toml
index 7509c9d..4e97225 100644
--- a/package/Project.toml
+++ b/package/Project.toml
@@ -5,12 +5,16 @@ version = "1.0.0-DEV"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
+DelimitedFiles = "1.9.1"
+GZip = "0.6.2"
 LinearAlgebra = "1.11.0"
 Printf = "1.11.0"
 Random = "1.11.0"
diff --git a/package/test/CpuInterpreterTests.jl b/package/test/CpuInterpreterTests.jl
index 02b7be7..19d564c 100644
--- a/package/test/CpuInterpreterTests.jl
+++ b/package/test/CpuInterpreterTests.jl
@@ -1,5 +1,9 @@
 using LinearAlgebra
 using BenchmarkTools
+using DelimitedFiles
+using GZip
+
+include("parser.jl") # to parse expressions from a file
 
 function test_cpu_interpreter(nrows; parallel = false)
     exprs = [
@@ -20,12 +24,12 @@ function test_cpu_interpreter(nrows; parallel = false)
 
     if parallel 
         # t_sec = @elapsed fetch.([Threads.@spawn interpret_cpu(exprs, X, p; repetitions=expr_reps) for i in 1:reps])
-        @btime parallel(exprs, X, p, expr_reps, reps)
-        println("~ $(round(30 * reps * expr_reps * nrows  / 1e9 / t_sec, digits=2)) GFLOPS ($(Threads.nthreads()) threads) ($(round(LinearAlgebra.peakflops(1000, eltype=Float32, ntrials=1) / 1e9, digits=2)) GFLOPS (peak, single-core))")
+        @btime parallel($exprs, $X, $p, $expr_reps, $reps)
+        # println("~ $(round(30 * reps * expr_reps * nrows  / 1e9 / t_sec, digits=2)) GFLOPS ($(Threads.nthreads()) threads) ($(round(LinearAlgebra.peakflops(1000, eltype=Float32, ntrials=1) / 1e9, digits=2)) GFLOPS (peak, single-core))")
     else
         # t_sec = @elapsed for i in 1:reps interpret_cpu(exprs, X, p; repetitions=expr_reps) end
-        @btime single(exprs, X, p, expr_reps, reps)
-        println("~ $(round(30 * reps * expr_reps * nrows  / 1e9 / t_sec, digits=2)) GFLOPS (single-core) ($(round(LinearAlgebra.peakflops(1000, eltype=Float32, ntrials=1) / 1e9, digits=2)) GFLOPS (peak, single-core))")
+        @btime single($exprs, $X, $p, $expr_reps, $reps)
+        # println("~ $(round(30 * reps * expr_reps * nrows  / 1e9 / t_sec, digits=2)) GFLOPS (single-core) ($(round(LinearAlgebra.peakflops(1000, eltype=Float32, ntrials=1) / 1e9, digits=2)) GFLOPS (peak, single-core))")
     end
     true
 end
@@ -45,3 +49,27 @@ end
 @test test_cpu_interpreter(1000, parallel=true) # start julia -t 6 for six threads
 @test test_cpu_interpreter(10000)
 @test test_cpu_interpreter(10000, parallel=true)
+
+
+function test_cpu_interpreter_nikuradse() 
+    data,varnames = readdlm("data/nikuradse_1.csv", ',', header=true);
+    X = convert(Matrix{Float32}, data)
+    
+    exprs = Expr[]
+    parameters = Vector{Vector{Float32}}()
+    varnames = ["x$i" for i in 1:10]
+    paramnames = ["p$i" for i in 1:20]
+    # data/esr_nvar2_len10.txt.gz_9.txt.gz has  ~250_000 exprs
+    # data/esr_nvar2_len10.txt.gz_10.txt.gz has ~800_000 exrps
+    GZip.open("data/esr_nvar2_len10.txt.gz_9.txt.gz") do io 
+        for line in eachline(io)
+            expr, p = parse_infix(line, varnames, paramnames)
+            
+            push!(exprs, expr)
+            push!(parameters, randn(Float32, length(p)))
+        end
+    end
+    
+    
+    interpret_cpu(exprs, X, parameters) # TODO: sufficient to do up to 10 repetitions per expression, 
+end
\ No newline at end of file
diff --git a/package/test/data/esr_nvar2_len10.txt.gz_1.txt.gz b/package/test/data/esr_nvar2_len10.txt.gz_1.txt.gz
new file mode 100644
index 0000000..e627ddb
Binary files /dev/null and b/package/test/data/esr_nvar2_len10.txt.gz_1.txt.gz differ
diff --git a/package/test/data/esr_nvar2_len10.txt.gz_10.txt.gz b/package/test/data/esr_nvar2_len10.txt.gz_10.txt.gz
new file mode 100644
index 0000000..9a19674
Binary files /dev/null and b/package/test/data/esr_nvar2_len10.txt.gz_10.txt.gz differ
diff --git a/package/test/data/esr_nvar2_len10.txt.gz_2.txt.gz b/package/test/data/esr_nvar2_len10.txt.gz_2.txt.gz
new file mode 100644
index 0000000..e7d1019
Binary files /dev/null and b/package/test/data/esr_nvar2_len10.txt.gz_2.txt.gz differ
diff --git a/package/test/data/esr_nvar2_len10.txt.gz_3.txt.gz b/package/test/data/esr_nvar2_len10.txt.gz_3.txt.gz
new file mode 100644
index 0000000..53bf509
Binary files /dev/null and b/package/test/data/esr_nvar2_len10.txt.gz_3.txt.gz differ
diff --git a/package/test/data/esr_nvar2_len10.txt.gz_4.txt.gz b/package/test/data/esr_nvar2_len10.txt.gz_4.txt.gz
new file mode 100644
index 0000000..7c8e16c
Binary files /dev/null and b/package/test/data/esr_nvar2_len10.txt.gz_4.txt.gz differ
diff --git a/package/test/data/esr_nvar2_len10.txt.gz_5.txt.gz b/package/test/data/esr_nvar2_len10.txt.gz_5.txt.gz
new file mode 100644
index 0000000..cc2df90
Binary files /dev/null and b/package/test/data/esr_nvar2_len10.txt.gz_5.txt.gz differ
diff --git a/package/test/data/esr_nvar2_len10.txt.gz_6.txt.gz b/package/test/data/esr_nvar2_len10.txt.gz_6.txt.gz
new file mode 100644
index 0000000..0798c05
Binary files /dev/null and b/package/test/data/esr_nvar2_len10.txt.gz_6.txt.gz differ
diff --git a/package/test/data/esr_nvar2_len10.txt.gz_7.txt.gz b/package/test/data/esr_nvar2_len10.txt.gz_7.txt.gz
new file mode 100644
index 0000000..c96d003
Binary files /dev/null and b/package/test/data/esr_nvar2_len10.txt.gz_7.txt.gz differ
diff --git a/package/test/data/esr_nvar2_len10.txt.gz_8.txt.gz b/package/test/data/esr_nvar2_len10.txt.gz_8.txt.gz
new file mode 100644
index 0000000..584b288
Binary files /dev/null and b/package/test/data/esr_nvar2_len10.txt.gz_8.txt.gz differ
diff --git a/package/test/data/esr_nvar2_len10.txt.gz_9.txt.gz b/package/test/data/esr_nvar2_len10.txt.gz_9.txt.gz
new file mode 100644
index 0000000..2cdc94a
Binary files /dev/null and b/package/test/data/esr_nvar2_len10.txt.gz_9.txt.gz differ
diff --git a/package/test/data/nikuradse_1.csv b/package/test/data/nikuradse_1.csv
new file mode 100644
index 0000000..30cfe6c
--- /dev/null
+++ b/package/test/data/nikuradse_1.csv
@@ -0,0 +1,363 @@
+r_k,log_Re,target
+507,4.114,0.456
+507,4.23,0.438
+507,4.322,0.417
+507,4.362,0.407
+507,4.362,0.403
+507,4.462,0.381
+507,4.491,0.38
+507,4.532,0.366
+507,4.568,0.365
+507,4.591,0.356
+507,4.623,0.347
+507,4.672,0.333
+507,4.69,0.324
+507,4.716,0.32
+507,4.763,0.307
+507,4.806,0.303
+507,4.851,0.292
+507,4.898,0.286
+507,4.94,0.278
+507,4.973,0.274
+507,5.009,0.274
+507,5.025,0.272
+507,5.049,0.27
+507,5.1,0.262
+507,5.143,0.26
+507,5.199,0.255
+507,5.236,0.253
+507,5.27,0.255
+507,5.281,0.253
+507,5.303,0.25
+507,5.326,0.252
+507,5.377,0.255
+507,5.43,0.253
+507,5.493,0.258
+507,5.534,0.26
+507,5.574,0.262
+507,5.608,0.29
+507,5.63,0.272
+507,5.668,0.272
+507,5.709,0.272
+507,5.756,0.278
+507,5.792,0.279
+507,5.833,0.283
+507,5.94,0.286
+507,5.965,0.288
+507,5.929,0.289
+507,5.954,0.288
+507,5.987,0.286
+252,4.21,0.4506
+252,4.279,0.4349
+252,4.465,0.3808
+252,4.507,0.3636
+252,4.549,0.3579
+252,4.597,0.3562
+252,4.644,0.3434
+252,4.778,0.3257
+252,4.82,0.3282
+252,4.916,0.3222
+252,4.987,0.3197
+252,5.057,0.321
+252,5.1,0.3228
+252,5.173,0.3197
+252,5.21,0.3276
+252,5.283,0.3322
+252,5.366,0.3416
+252,5.494,0.3504
+252,5.58,0.3562
+252,5.623,0.3602
+252,5.702,0.3636
+252,4.708,0.3371
+252,5.305,0.3328
+252,5.544,0.3562
+252,5.787,0.3661
+252,4.748,0.3335
+252,4.869,0.3228
+252,4.954,0.321
+252,5.134,0.321
+252,5.255,0.3294
+252,5.415,0.3434
+252,5.58,0.3551
+252,5.748,0.3608
+252,5.845,0.3666
+252,5.881,0.3688
+252,5.924,0.3727
+252,5.967,0.3705
+252,5.991,0.3716
+60,3.653,0.593
+60,3.7,0.577
+60,3.74,0.571
+60,3.785,0.56
+60,3.851,0.544
+60,3.869,0.531
+60,3.909,0.512
+60,3.949,0.512
+60,3.996,0.507
+60,4.057,0.494
+60,4.09,0.49
+60,4.161,0.494
+60,4.236,0.487
+60,4.29,0.487
+60,4.391,0.481
+60,4.412,0.489
+60,4.512,0.49
+60,4.54,0.487
+60,4.553,0.498
+60,4.58,0.493
+60,4.609,0.507
+60,4.654,0.504
+60,4.665,0.507
+60,4.699,0.509
+60,4.74,0.517
+60,4.769,0.52
+60,4.813,0.528
+60,4.849,0.526
+60,4.93,0.543
+60,4.954,0.534
+60,5.034,0.543
+60,5.155,0.543
+60,5.083,0.545
+60,5.185,0.55
+60,5.231,0.537
+60,4.875,0.535
+60,4.924,0.534
+60,4.954,0.542
+60,5.052,0.535
+60,5.033,0.54
+60,5.13,0.545
+60,5.17,0.55
+60,5.196,0.547
+60,5.23,0.568
+60,5.258,0.551
+60,5.283,0.555
+60,5.312,0.551
+60,5.35,0.555
+60,5.408,0.55
+60,5.47,0.555
+60,5.497,0.543
+60,5.515,0.551
+60,5.549,0.55
+60,5.554,0.558
+60,5.575,0.551
+60,5.6,0.55
+60,5.621,0.56
+60,5.625,0.543
+60,5.641,0.543
+60,5.655,0.55
+60,5.659,0.551
+60,5.668,0.56
+60,5.691,0.553
+60,5.714,0.551
+60,5.748,0.558
+60,5.757,0.55
+60,5.789,0.551
+60,5.836,0.547
+60,5.865,0.555
+60,5.914,0.553
+60,5.916,0.55
+60,5.945,0.551
+60,5.962,0.555
+15,3.77,0.696
+15,3.82,0.699
+15,3.855,0.707
+15,3.905,0.712
+15,3.955,0.717
+15,4,0.73
+15,4.041,0.734
+15,4.076,0.736
+15,4.079,0.744
+15,4.114,0.751
+15,4.133,0.74
+15,4.179,0.744
+15,4.196,0.754
+15,4.27,0.76
+15,4.29,0.756
+15,4.314,0.769
+15,4.34,0.763
+15,4.366,0.778
+15,4.386,0.772
+15,4.41,0.772
+15,4.425,0.782
+15,4.466,0.785
+15,4.52,0.78
+15,4.59,0.781
+15,4.63,0.777
+15,4.725,0.78
+15,4.811,0.781
+15,4.865,0.777
+15,4.885,0.776
+15,4.965,0.779
+15,5,0.781
+15,5.042,0.78
+15,5.098,0.781
+15,5.155,0.778
+15,5.179,0.781
+15,5.285,0.779
+15,4.44,0.775
+15,4.5,0.777
+15,4.54,0.778
+15,4.596,0.78
+15,4.685,0.781
+15,4.722,0.777
+15,4.845,0.775
+15,4.869,0.778
+15,4.929,0.78
+15,4.949,0.779
+15,5.002,0.777
+15,5.005,0.775
+15,5.097,0.778
+15,5.139,0.783
+15,5.156,0.784
+15,5.22,0.777
+15,5.236,0.78
+15,5.31,0.778
+15,5.36,0.775
+15,5.41,0.78
+15,5.446,0.78
+15,5.455,0.777
+15,5.515,0.781
+15,5.567,0.778
+15,5.613,0.78
+15,5.69,0.784
+15,5.834,0.781
+15,5.882,0.777
+15,5.959,0.778
+15,6.008,0.78
+15,5.793,0.78
+15,5.857,0.777
+15,5.93,0.778
+15,5.987,0.78
+126,3.63,0.594
+126,3.675,0.588
+126,3.715,0.576
+126,3.76,0.566
+126,3.81,0.552
+126,3.833,0.564
+126,3.895,0.532
+126,3.925,0.515
+126,3.95,0.503
+126,3.965,0.498
+126,4.015,0.491
+126,4.111,0.471
+126,4.196,0.451
+126,4.265,0.435
+126,4.33,0.424
+126,4.386,0.415
+126,4.425,0.412
+126,4.47,0.4
+126,4.496,0.396
+126,4.511,0.4
+126,4.55,0.393
+126,4.62,0.392
+126,4.697,0.391
+126,4.76,0.4
+126,4.82,0.403
+126,4.91,0.408
+126,4.985,0.414
+126,5.057,0.422
+126,5.121,0.424
+126,5.164,0.43
+126,5.591,0.45
+126,5.616,0.453
+126,5.655,0.447
+126,5.675,0.45
+126,5.708,0.445
+126,5.736,0.452
+126,5.756,0.445
+126,5.775,0.445
+126,5.798,0.45
+126,5.831,0.45
+126,5.835,0.446
+126,5.874,0.45
+126,5.894,0.447
+126,5.935,0.45
+126,5.961,0.444
+126,5.97,0.449
+126,5.987,0.447
+126,4.95,0.43
+126,5.049,0.432
+126,5.021,0.415
+126,5.1,0.422
+126,5.13,0.422
+126,5.179,0.43
+126,5.196,0.43
+126,5.225,0.435
+126,5.225,0.43
+126,5.25,0.436
+126,5.274,0.438
+126,5.29,0.438
+126,5.31,0.436
+126,5.33,0.439
+126,5.35,0.439
+126,5.366,0.444
+126,5.393,0.444
+126,5.423,0.446
+126,5.432,0.447
+126,5.455,0.45
+126,5.476,0.452
+126,5.501,0.447
+126,5.525,0.447
+126,5.56,0.45
+30.6,3.672,0.592
+30.6,3.708,0.59
+30.6,3.748,0.592
+30.6,3.763,0.597
+30.6,3.785,0.583
+30.6,3.826,0.585
+30.6,3.869,0.596
+30.6,3.881,0.578
+30.6,3.929,0.578
+30.6,3.935,0.583
+30.6,3.978,0.578
+30.6,4.009,0.585
+30.6,4.049,0.583
+30.6,4.079,0.592
+30.6,4.124,0.59
+30.6,4.13,0.599
+30.6,4.19,0.599
+30.6,4.27,0.609
+30.6,4.29,0.618
+30.6,4.309,0.612
+30.6,4.584,0.639
+30.6,4.653,0.644
+30.6,4.799,0.647
+30.6,4.9,0.656
+30.6,4.965,0.656
+30.6,5.029,0.652
+30.6,5.068,0.65
+30.6,5.134,0.65
+30.6,5.176,0.65
+30.6,4.425,0.637
+30.6,4.44,0.63
+30.6,4.56,0.637
+30.6,4.636,0.647
+30.6,4.74,0.654
+30.6,4.83,0.654
+30.6,4.855,0.661
+30.6,4.99,0.657
+30.6,5.1,0.652
+30.6,5.24,0.657
+30.6,5.275,0.657
+30.6,5.323,0.647
+30.6,5.473,0.657
+30.6,5.655,0.652
+30.6,4.934,0.656
+30.6,5.068,0.657
+30.6,5.17,0.659
+30.6,5.223,0.656
+30.6,5.255,0.652
+30.6,5.342,0.657
+30.6,5.344,0.657
+30.6,5.394,0.659
+30.6,5.428,0.659
+30.6,5.444,0.661
+30.6,5.516,0.657
+30.6,5.541,0.659
+30.6,5.559,0.657
+30.6,5.776,0.659
+30.6,5.81,0.659
+30.6,5.863,0.657
+30.6,5.916,0.659
+30.6,5.962,0.65
+30.6,6,0.659
diff --git a/package/test/parser.jl b/package/test/parser.jl
new file mode 100644
index 0000000..6ec7684
--- /dev/null
+++ b/package/test/parser.jl
@@ -0,0 +1,294 @@
+## Parser for (ESR) expressions in infix format
+
+mutable struct Parser
+    const str::AbstractString               # string to be parsed
+    pos::Int64                        # current position in string
+    sy::Union{AbstractString,Nothing} # current lookahead symbol
+    const pSy::Symbol
+    const xSy::Symbol
+    const varnames::Vector{<:AbstractString}
+    const paramnames::Vector{<:AbstractString}
+    const coeff::Vector{Float64}
+    const numbers_as_parameters::Bool
+    const integers_as_constants::Bool # TODO rename and implement as rationals_as_constants
+    
+    # The kwparam numbers_as_parameters allows to include coefficient values directly in the expression and the values are parsed as parameters
+    # In this mode the suffix 'f' allows to mark constants. E.g. 3 * x ^ 2f  would create the parameterized expression a0*x^2 with 2 a constant value.
+    function Parser(str::AbstractString, varnames::Vector{<:AbstractString}, paramnames::Vector{<:AbstractString}; numbers_as_parameters=false, integers_as_constants=false)
+        if numbers_as_parameters && length(paramnames) > 0 
+            error("the parser does not support paramnames when numbers_as_parameters=true")
+        end
+        if !numbers_as_parameters && integers_as_constants
+            error("Set numbers_as_parameters=true to parse integers_as_constants")
+        end
+
+        p = new(lowercase(str), 1, nothing, :p, :x, varnames, paramnames, Vector{Float64}(), numbers_as_parameters, integers_as_constants)
+        next_symbol!(p)
+        return p;
+    end
+end
+
+# recursive descent parser
+# scanner is also defined in this file
+
+# LL(1) grammar:
+# G(Expr):
+# Expr      = Term { ('+' | '-') Term }
+# Term      = Fact { ('*' | '/') Fact }
+# Fact      = { '+' | '-' }
+#             (ident | number | parameter
+#              | '(' Expr ')'
+#              | ident ParamList                           // function call
+#             ) [ ('**' | '^') Fact ]
+# ParamList = '(' Expr { ',' Expr } ')'
+
+
+
+# scanner
+
+
+function parse_infix(exprStr::AbstractString, varnames::Vector{<:AbstractString}, paramnames::Vector{<:AbstractString}; 
+                     numbers_as_parameters = false, integers_as_constants = false)::Tuple{Expr, Vector{Float64}}
+    parser = Parser(exprStr, varnames, paramnames; 
+                    numbers_as_parameters = numbers_as_parameters, integers_as_constants = integers_as_constants)
+    body = parse_expr!(parser)
+    expr = Expr(:->, Expr(:tuple, :x, :p), body) # :((x,p) -> $body) 
+    (expr, parser.coeff)    
+end
+
+function parse_expr!(p::Parser)
+    t1 = parse_term!(p)
+    while p.sy == "+" || p.sy == "-"
+        if p.sy == "+"
+            next_symbol!(p)
+            t2 = parse_term!(p)
+            t1 = :($t1 + $t2) # add_simpl(t1, t2)
+        else
+            next_symbol!(p)
+            t2 = parse_term!(p)
+            t1 = :($t1 - $t2) # sub_simpl(t1, t2)
+        end
+    end
+    return t1
+end
+
+function parse_term!(p::Parser)
+    f1 = parse_factor!(p)
+    while p.sy == "*" || p.sy == "/"
+        if p.sy == "*"
+            next_symbol!(p)
+            f2 = parse_factor!(p)
+            f1 = :($f1 * $f2) # mul_simpl(f1, f2)
+        else
+            next_symbol!(p)
+            f2 = parse_factor!(p)
+            f1 = :($f1 / $f2) # div_simpl(f1, f2)
+        end
+    end
+    return f1
+end
+
+# Fact      = { '+' | '-' }
+#             (constant | parameter
+#              | '(' Expr ')'
+#              | ident [ ParamList ]        variable or function call
+#             ) [ ('**' | '^') Fact ]
+# ParamList = '(' Expr { ',' Expr } ')'
+
+function parse_factor!(p::Parser)
+    sign = 1.0
+    
+    while p.sy == "+" || p.sy == "-"
+        if p.sy == "-"
+            sign = sign * -1.0
+        end
+        next_symbol!(p)
+    end
+    
+    factor = 1.0
+    
+    if isident(p.sy)
+        ident = p.sy
+        next_symbol!(p)
+        if p.sy == "("
+            parameters = parse_paramlist!(p)
+            
+            if ident == "sqr"
+                # convert sqr(x) call to x**2 (so that we don't have to update the interpreters)
+                factor = Expr(:call, func_symbol("pow"), parameters..., 2.0)
+            else
+                factor = Expr(:call, func_symbol(ident), parameters...)
+            end
+        else
+            idx =  findfirst(p -> p==ident, p.varnames)
+            if !isnothing(idx)
+              factor = Expr(:ref, p.xSy, idx) 
+            elseif !p.numbers_as_parameters  # only if paramnames are given
+               idx =  findfirst(p -> p==ident, p.paramnames)
+
+               # replace parameter variables with access to coefficient vector (initialized to zero)
+               if !isnothing(idx)
+                    factor = Expr(:ref, p.pSy, idx)
+                    push!(p.coeff, 0.0)
+                else 
+                    error("undefined symbol $ident")
+                end
+            else
+                error("undefined variable $ident")
+            end            
+        end
+        
+    elseif isnumber(p.sy)
+        if p.numbers_as_parameters 
+            numStr = p.sy
+            val = parse(Float64, numStr)
+            next_symbol!(p)
+            if p.sy == "f"
+                # constant
+                factor = sign * val    # numbers are parsed without sign (if we parsed a sign above then we can include this in the constant here)
+                sign = 1.0
+                next_symbol!(p)
+            elseif p.integers_as_constants && isinteger(val)
+                # integers are parsed as constants
+                factor = sign * val    # numbers are parsed without sign (if we parsed a sign above then we can include this in the constant here)
+                sign = 1.0
+            else 
+                # parameter
+                factor = new_param!(p, sign * val)
+                sign = 1.0
+            end
+        else 
+            # otherwise all numbers are parsed as constants
+            numStr = p.sy
+            next_symbol!(p)
+            
+            if p.sy == "//"
+                num = parse(Int64, numStr)
+                next_symbol!(p)
+                denom = parse(Int64, p.sy)
+                val = num // denom
+                next_symbol!(p)
+            else
+                val = parse(Float64, numStr)
+            end
+            
+            factor = sign * val
+            sign = 1.0            
+        end
+        
+    elseif p.sy == "("
+        next_symbol!(p)
+        factor = parse_expr!(p)
+        expect_and_next!(p, ")")
+        
+    else
+        error("cannot parse expression")
+    end
+    
+    if p.sy == "**" || p.sy == "^"
+        next_symbol!(p)
+        exponent = parse_factor!(p)
+        factor = :($factor ^ $exponent) # pow_simpl(factor, exponent)
+    end
+
+    if sign == -1
+        :(-$factor)
+    else
+        factor
+    end
+end
+
+function parse_paramlist!(p::Parser)::Vector
+    parameters = Vector()
+    expect_and_next!(p, "(")
+    push!(parameters, parse_expr!(p))
+    while p.sy == ","
+        next_symbol!(p)
+        push!(parameters, parse_expr!(p))
+    end
+    expect_and_next!(p, ")")
+    return parameters
+end
+
+function expect_and_next!(p::Parser, expectedSy::AbstractString)
+    if p.sy != expectedSy
+        error("expected: $(expectedSy) at column $(p.pos)")
+    else
+        next_symbol!(p)
+    end   
+end
+
+function new_param!(p::Parser, val::Float64)::Expr
+    push!(p.coeff, val)
+    return Expr(:ref, p.pSy, length(p.coeff))
+end
+
+
+function isident(s::AbstractString)::Bool
+    return s != "nan" && s != "inf" && !isnothing(match(r"^[_a-zA-Z][_a-zA-Z0-9]*$", s))
+end
+
+function isnumber(s::AbstractString)::Bool
+    return !isnothing(tryparse(Float64, s))
+end
+
+function variable_index(p::Parser, str::AbstractString)
+    return findfirst(s->s==str, p.varNames)
+end
+
+function func_symbol(id::AbstractString)
+    if id == "pow"
+        return :^;
+    else
+        return Symbol(id)
+    end
+end
+
+function next_symbol!(p::Parser)
+    s = p.str
+    pos = p.pos
+    # skip whitespace
+    while pos <= length(s) && isspace(s[pos])
+        pos += 1
+    end
+
+    if pos > length(s)
+        p.sy = nothing
+        p.pos = pos
+        return
+    end
+    
+    if isdigit(s[pos]) # numbers
+        m = match(r"(\d+([.]\d*)?([eE][+-]?\d+)?|[.]\d+([eE][+-]?\d+)?)", s, pos) # match floating point number
+        pos += length(m[1]) # get the whole match
+        p.sy = m[1]
+    elseif isletter(s[pos]) # identifiers
+        idStr = string(s[pos])
+        pos += 1
+        while pos <= length(s) && (isdigit(s[pos]) || isletter(s[pos]) || s[pos] == '_')
+            idStr = idStr * s[pos]
+            pos += 1
+        end
+        p.sy = idStr
+    elseif s[pos] == '*'
+        pos += 1
+        p.sy = "*"
+        if s[pos] == '*'
+            p.sy = "**"
+            pos += 1
+        end
+    elseif s[pos] == '/'
+        pos += 1
+        p.sy = "/"
+        if s[pos] == '/'
+            p.sy = "//"
+            pos += 1
+        end
+    else
+        p.sy = string(s[pos]) # single character symbol
+        pos += 1
+    end
+
+    p.pos = pos
+#    println((p.sy, pos)) # for debugging
+end