ctx->entries->entry-2016123118

numerical mathematical functions in scheme (1)

Let's imagine a function. For example f(x) = x^2. Symbolically we're done. We defined an operation that takes values from one space and maps them to another. But if we want to work with it numerically we better be able to pass it some of these values and get back the transformed ones. So we might want a function that generates values first, on a linear space. Let's start with an extremely naive implementation. We define a recursive function that takes the limits of the space we want to generate and the step to increase. if we go past the right limit we stop.

(define genlist
        (lambda (start end step)
                (if (> start end) 
                    '()
                    (cons start (genlist (+ start step) end step)))))

>(genlist 0 1 0.1)
(0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0)
>(genlist 0 1 0.01)
(0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4 0.41 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.6 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.7 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8 0.81 0.820000000000001 0.830000000000001 0.840000000000001 0.850000000000001 0.860000000000001 0.870000000000001 0.880000000000001 0.890000000000001 0.900000000000001 0.910000000000001 0.920000000000001 0.930000000000001 0.940000000000001 0.950000000000001 0.960000000000001 0.970000000000001 0.980000000000001 0.990000000000001)

`

We see that due to the floating point nature of the numbers we're using accurancy is not controlled. We can think of a way to control it. We are going to need a function to round-off to the number of decimal points we want, and a function to detect the current number of digits. The function provided here employs a string trick. it is also possible to detect it numerically by succesful divisions or multiplications by 10 , checking the modulo of the operation with 10 and tracking the number of operations with an accumulator.

(define (round-off z n)
  (let ((power (expt 10 n)))
    (/ (round (* power z)) power)))
>(round-off 12.124124124 3)
12.124

(define (number-of-digits x)
  (let ((parts (string-split (number->string x) ".")))
    (if (= 2 (length parts))
        (+ (string-length (car parts)) (string-length (cadr parts)))
        (string-length (car parts)))))

>(number-of-digits 0.10111)
6

Now we can change our genlist function to maintain the same nunber of digits in the generated number as the step.

(define genlist
        (lambda (start end step)
                (if (> start end) '()
                    (cons start (genlist (round-off (+ start step) (number-of-digits step)) end step)))))
>(genlist 0 1 0.01)
(0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4 0.41 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.6 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.7 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1.0)

` Ok now we need a function to time the execution of our functions. They won't be extremely performant because of the use of plain lists to represent our space. We will generate 10k , 100k and 1m numbers to see how it scales

(define timeit 
   (lambda (f) 
           (let ((start 0) 
                 (end 0)) 
             (set! start (current-milliseconds)) 
             (f) 
             (set! end (current-milliseconds)) 
             (- end start))))
>(len (genlist 0 100 0.01))
10001
>(timeit (lambda () (genlist 0 100 0.01)))
137.0
>(timeit (lambda () (genlist 0 1000 0.01)))
1235.0
>(timeit (lambda () (genlist 0 10000 0.01)))
14763.0

That's quite bad. But we can use vectors to the rescue.

(use vector-lib) ;;chicken scheme specific

(define genlist2
    (lambda (start end step)
          (let* ((startcnt start)
                 (nd (number-of-digits step))
                 (times (inexact->exact (ceiling (/ (- end start) step))))
                (incs (lambda (s) 
                         (let ((ret startcnt)) 
                           (set! startcnt (round-off (+ s startcnt) nd)) 
                           (if (< ret end)
                               ret)))))
       (vector-unfold (lambda (x) (incs step)) times))))

>(timeit (lambda () (genlist2 0 1000 0.1)))
39.0
>(timeit (lambda () (genlist2 0 10000 0.1)))
390.0
>(timeit (lambda () (genlist2 0 100000 0.1)))
4086.0

That's an one order of magnitude improvement in generating linear spaces with a lot of elements. On future posts we will investigate how we map functions on those spaces and plot them.

Get in touch: dsp ..at. 2f30 .dot.. org / irc.freenode.org/#2f30