# Derivative, dual numbers and Guile

introduction to some basics of Guile

Note:Back in 2019, I wrote «Calcul de Dérivées via les nombres duaux» usingthe`C`

language. The aim was to illustrate how computing helps about concreteengineering problems. The same example is reused here to show Scheme programming.The context of this post is the monthly Café Guix. Join the fun!

In this post, we are going to illustrate some GNU Guile concepts. The programming language is based on Scheme and our aim is to be as functional and declarative as possible.

The jungle of the Lisp family programming language can be confusing, at first.
Let quickly summarized. Lisp, initially named as « *LISt Processor* », is easily
recognisable by its fully parenthesized prefix notation syntax. The first
Lisp had been developed in 1958 and the modern variants inherit many components,
making them all cousins. Scheme, introduced in the 1970s, is one variant and
then various implementations of this programming language had been released;
to name some Racket, MIT/GNU Scheme or Chez Scheme and… GNU Guile. Although
language specifications are written down, all these implementations have their
own specificity.

First thing first, the heart of any Lisp is **S-expression**, defined by an atom
or an expression of the form `(x y …)`

where `x`

and `y`

are S-expression. Some
examples for the usual:

`(list 1 2 "foo")`

is the list containing the number`1`

,`2`

and the string`"foo",`

`(define x (list 1 2 "foo"))`

binds the list to the symbol named`x`

,`(lambda (a) (+ 1 a))`

creates a procedure taking one argument (named`a`

) and then returning the addition by one,`(define f (lambda (a) (+ 1 a)))`

binds the procedure to the symbol named`f`

,- the call of the procedure is simply:
`(f 42)`

.

Instead of,

(define f (lambda (a) (+ 1 a)))

it is possible to use the equivalent sugar,

(define (f a) (+ 1 a))

We have enough… well we will introduce later what we need. Let’s go!

## The problem we try to solve (elementary maths)

This part is not necessary for the programming part. It just helps to have some context about what we are going to implement.

Let a mathematical function \(f\) from the real numbers to the real numbers,

\begin{equation*} \begin{array}{lclcl} f &:& \mathbb{R} &\longrightarrow& \mathbb{R}\\ & & x &\mapsto& y=f(x) \end{array} \end{equation*}and the derivative, noted \(f^\prime\), is defined by,

\begin{equation} \label{org895e05b} f^\prime(x) = \underset{h\rightarrow 0}{\lim} \frac{f(x+h) - f(x)}{h}. \end{equation}Picking one value, say \(a\), it is straightforward to compute the function for this value and the new computed value is \(b=f(a)\). When the function is analytic, symbolic computations are possible using pen and paper or some computer algebra system. From this symbolic formula of \(f^\prime\), the value of the derivative can be then computed to the value \(a\). However, most of the time, \(f\) is not analytic and then how to compute \(f^\prime(a)\)?

For sure, it is possible to reuse the definition \eqref{org895e05b} and approximate it by finite difference (fixing a small enough value for \(h\)). The main issue is the choice of such value \(h\) and the accuracy strongly depends on this choice.

### Differentiation rules

We learn in high school the differentiation rules which reads,

\begin{equation} \label{org03485ea} \begin{array}{lcl} (f+g)^\prime &=& f^\prime + g^\prime\\ (f*g)^\prime &=& f^\prime * g + f* g^\prime\\ (f-g)^\prime &=& f^\prime - g^\prime\\ \displaystyle \left(\frac{f}{g}\right)^\prime &=& \displaystyle \frac{f^\prime * g - f* g^\prime}{g^2}. \end{array} \end{equation}Nothing fancy!

### Dual numbers

Although old, dual numbers are not well-know, from what I have seen. Let define a dual number \(x\) with two real numbers \(v\) and \(d\), such that,

\begin{equation*} x = (v, d) \end{equation*}and similarly \(y=(u, p)\) another dual number. The addition, multiplication, subtraction and division,

\begin{equation} \label{orgc82c208} \begin{array}{lcll} \texttt{add}(x, y) &=& (v+u ,& d+p) \\ \texttt{mul}(x, y) &=& (u*v ,& d*v ~+~ u*p) \\ \texttt{sub}(x, y) &=& (u-v ,& d-p) \\ \texttt{div}(x, y) &=& (\displaystyle \frac{u}{v} ,& \displaystyle \frac{d*v ~-~ u*p}{p^2}). \end{array} \end{equation}Wait! It looks really similar, well even, the same identical with the differentiation rules. Indeed, let denote this number as \((\texttt{value}, \texttt{deriv})\) where,

`value`

behaves as real numbers,`deriv`

behaves as derivatives.

#### Example

Let the function \(f(x)=x^{2}\) and we are able to compute the derivative: \(f^{\prime}(x)=2x\). Let the dual number \(z=(a, 1.)\). Now, we substitute \(z\) in \(f\) applying the arithmetic operation of dual number,

\begin{equation*} \texttt{mul}(z,z) = (a*a, 1.*a + a*1.) = (a*a, 2*a). \end{equation*}Tadam! Neat, isn’t it? Ok, we are going to implement these dual numbers in Guile.

## Dual number implemented with Guile

#### Record

If we think as a `C`

programmer, we would like to use some `struct`

to represent
these dual numbers; something like,

typedef struct { char name[MAXCHAR]; float value; float deriv; } dual;

In Guile, *structures* are named *record*. And a record is the mechanism to
define new disjoint type. Guile provides a syntax for defining new record but
it can be extended adding some facilities when loading the module SRFI-9.
Yeah, the names based on number is impractical–the best is to rely on this
page of the Guile manual.

The `struct`

reads,

(use-modules (srfi srfi-9)) (define-record-type <dual> ;1 name (dual name value deriv) ;2 constructor dual? ;3 predicate (name dual-name) ;4 field and accessor (value dual-value) ;5 idem (deriv dual-deriv)) ;6 idem

where, some explanations,

`<dual>`

is the name of the record; it can be anything and`<>`

is a convention,`(dual name value deriv)`

is the constructor; to create a new dual number, it reads, e.g.,`(define z (dual "whatever" 1 2))`

,`dual?`

is a predicate which tests if it is or not the type of`<dual>`

; for instance,`(dual? z)`

return`#t`

(true) whereas`(dual? 42)`

returns`#f`

(false)`(name dual-name)`

where`name`

stores the value provided in the constructor and`dual-name`

is the accessor of that field`name`

,- Idem for
`value`

, - Idem for
`deriv`

.

scheme@(guile-user)> (define z (dual "whatever" 1 2)) scheme@(guile-user)> z $3 = #<<dual> name: "whatever" value: 1 deriv: 2> scheme@(guile-user)> (dual? z) $4 = #t scheme@(guile-user)> (dual? 42) $5 = #f scheme@(guile-user)> (dual-name z) $6 = "whatever" scheme@(guile-user)> (dual-value z) $7 = 1

#### Format

Now, we would like to show a dual number, for instance, get this output,

scheme@(guile-user)> (dual-show z) "whatever" Value: 1.0E+0 Deriv: 2.0E+0 $8 = #<<dual> name: "whatever" value: 1 deriv: 2>

For one part, it is simply the use of the accessors. For the other part, it
is `format`

– similar to `printf`

from `C`

, though escape sequences are marked with
`~`

instead of `%`

. The procedure `dual-show`

reads,

(define (dual-show z) (format #t "~s~%Value: ~e~%Deriv: ~e~%" (dual-name z) (dual-value z) (dual-deriv z)) z)

We return the argument for returning something. That's arbitrary at this point.

#### Local variable and optional argument

Previously, we used the dual number \((a, 1.)\), so we would like a procedure to build it. From what we know, we would write,

(define (dual-with-one name val) (dual name val 1.))

but instead, let introduce the `let`

binding which allows to declare local
variable bindings,

(define (variable name val) (let ((local 1.0)) (dual name val local)))

Here, `local`

is a *local* variable where the value is set to `1.`

. Note,
conditional as `if`

could be there, for instance,

(let ((local (if name name "c"))) ...)

which means, if `name`

is set, then use this value for the variable `local`

, else
use the string `"c"`

. Hum, we would like to define a new procedure but with
an optional argument. It is provided by the variant `define*`

, for instance,

(define* (constant val #:optional name) (let ((local (if name name "c"))) (dual local val 0)))

where `#:optional name`

means it is… optional. Calling `(constant 2)`

set the
default value `#f`

(false) for `name`

. Let as an exercise to set another default
value than `#f`

, for example 42. A quick remark: atoms with `#:`

as `#:optinonal`

are named keyword.

**Another quick remark:** append a star often means it is a variant; `define*`

is an example, it is a variant of `define`

. For instance, `let*`

is also a
variant of `let`

: it allows to create local variables which depend on previous
declared ones,

```
(let* ((one 1)
(two (+ 1 one))
(three (1+ two)))
(* 2 three))
```

Note that `1+`

is not a typo.

#### Implementation of dual number addition

What is still missing is how to manipulate strings. For instance, when adding
two dual numbers, we would like to append their name around the sign plus
(+). Maybe `string-append`

, let try:

scheme@(guile-user)> (define name1 "x") scheme@(guile-user)> (define name2 "y") scheme@(guile-user)> (string-append "(" name1 "+" name2 ")") $9 = "(x+y)"

Neat! Well, thus for implementing the addition of dual numbers, we only need to access to all the fields and construct a new dual number,

(define (dual-add z1 z2) (let ((name1 (dual-name z1)) (value1 (dual-value z1)) (deriv1 (dual-deriv z1)) (name2 (dual-name z2)) (value2 (dual-value z2)) (deriv2 (dual-deriv z2))) (dual (string-append "(" name1 "+" name2 ")") (+ value1 value2) (+ deriv1 deriv2))))

Well, something does not appear smooth. It could be nice to easy the unwrap way.

#### Pattern matching

Pattern matching is a powerful feature by functional programming languages.
It allows to match an object against several patterns and extract the elements
that make it up. It requires the module `(ice-9 match)`

. Let take some
examples to fix the ideas. The procedure `example`

returns,

scheme@(guile-user)> (example (list 'first 2 'third)) $10 = third scheme@(guile-user)> (example (list 'first 'second 3)) $11 = 9 scheme@(guile-user)> (example (list 1 2 3)) $12 = 5 scheme@(guile-user)> (example (list 1 2 3 5)) $13 = (1 2 3 5) scheme@(guile-user)> (example "that") $14 = "that"

How do you implement it? Using many conditionals? Nope, using pattern matching! The syntax reads,

(define (example object) (match object (('first 2 something) ;1 something) ((_ 'second n) ;2 (* n n)) ((1 b c) ;3 (+ b c)) (x ;4 x)))

where `object`

is extracted following patterns:

- if
`object`

matches the pattern`('first 2 something)`

where`something`

is the name of a local variable, then return the value of that variable, if

`object`

matches the pattern`(_ 'second number)`

where`_`

(underscore) means anything and`number`

is the name of a local variable, then return the multiplication with this number. It would be more correct to also check, when unwrapping that the third element is a number (predicate`number?`

),((_ 'second (? number? n)) (* n n))

- etc.

As you see, it is really powerful because on one hand, it extracts the elements and on the other hands it checks the pattern. Ok, another example. Lists are fundamental and primitive operations in Lisp are car and cdr which returns the first element of a list and the rest except the first.

scheme@(guile-user)> (car (list 'first 2 'third)) $15 = first scheme@(guile-user)> (cdr (list 'first 2 'third)) $16 = (2 third)

The pattern matching allow *ellipsis* (`...`

). It means that the previous pattern
may be matched zero or more times in a list. Hence, a procedure using `match`

returning the first element.

(define (car-match lst) (match lst ((head tail ...) head))) ;here

Instead of `head`

(line `;here`

), try to replace by `tail`

(and rename the
procedure, e.g., `cdr-match`

).

scheme@(guile-user)> (cdr-match (list 'first 2 'third)) $17 = (2 third)

So far, so good! Let implement the subtraction.

#### Implementation of dual number subtraction

Once the pattern matching documentation careful read, instead of several access to the record field, now it appears straightforward to pattern match using the record. It reads,

(define (dual-sub z1 z2) (match z1 (($ <dual> name1 value1 deriv1) (match z2 (($ <dual> name2 value2 deriv2) (dual (string-append "(" name1 "-" name2 ")") (- value1 value2) (- deriv1 deriv2)))))))

Much more elegant than `dual-add`

, isn’t it? The two remaining arithmetic
operations (multiplication and division) are let as an exercise.

#### Implementation of dual number exponential

Since the exponential is an analytic function and we know its derivative (\((e^{u})^{\prime}=u^{\prime}e^{u}\)), the efficient implementation is to exploit such, for instance,

(define (dual-exp z) (match z (($ <dual> name value deriv) (let ((exp-value (exp value))) (dual (format #f "~~e^(~a)" (dual-name z)) exp-value (* deriv exp-value))))))

Well, we would like to practise Guile, let implement the exponential using only the four arithmetic operations. The exponential function is defined by power series with the formula,

\begin{equation*} e^{x} = \sum_{k=0}^{\infty} \frac{x^{k}}{k!} = 1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \dots \end{equation*}
We will not compute until infinity, say only the first twenty terms. The
first thinking is to iterate using a for-loop. However, it does not respect
the functional style–instead recursion is the functional style. Well, we
would say using plain words: if \(k\) is smaller than `max-iter`

(say 20), then
update and apply again (recursion-loop) this test, else return the result.
Somehow, it would translate to,

(if (<= k max-iter) (loop (+1 k) expx xk kfac) (dual (format #f "~~e^(~a)" (dual-name x)) (dual-value expx) (dual-deriv expx)))

where `expx`

, `xk`

and `kfac`

are not well-defined yet. Knowing \(x^{k}\) or \(k!\),
the update is provided by,

and thus, knowing previous summation (`expx`

), we just compute the \(k\) -th term
(`term_k`

) and add it. Somehow, this `loop`

procedure could be defined with,

(define (loop k expx xk kfac) (let* ((kk (constant k)) ;convert integer to dual (xk (dual-mul xk x)) (kfac (dual-mul kfac kk)) (term_k (dual-div xk kfac)) (expx (dual-add expx term_k))) (if (<= k max-iter) (loop (1+ k) expx xk kfac) (dual (format #f "~~e^(~a)" (dual-name x)) (dual-value expx) (dual-deriv expx)))))

And what remains is the initialization and/or initial call. It can be
confusing at first but this `loop`

procedure can be defined via `let`

(see
iteration mechanism).

(define (dual-exponential x) (let ((max-iter 20)) (let loop ((k 1) (expx (constant 1)) (xk (constant 1)) (kfac (constant 1))) (let* ((kk (constant k)) (xk (dual-mul xk x)) (kfac (dual-mul kfac kk)) (term_k (dual-div xk kfac)) (expx (dual-add expx term_k))) (if (<= k max-iter) (loop (1+ k) expx xk kfac) (dual (format #f "~~e^(~a)" (dual-name x)) (dual-value expx) (dual-deriv expx)))))))

Confusing… yeah! Functional style appears unexpected at first; it is true for OCaml, Haskell or similar when coming from imperative style as Python or similar.

Because we have fun, let demo some other nice features.

## Extra

#### High-order function

In functional style, procedure (function) is a first-class citizen. Scheme makes it easy to compose them. Consider two procedures,

(define (f x) (+ x 2)) (define (g x) (* x 2))

and composing them reads,

(define (h x) (f (g x)))

or simply,

(define fog (compose f g))

And `(fog 5)`

returns 12 as expected. Now, consider the new procedure `f`

taking
a `value`

argument instead of plain 2.

(define (f x value) (+ x value))

And we would like to be able to vary this `value`

in the composition `fog`

.
First, note that the composition does not match. Well, to understand
correctly, we are going to write using good ol’ `lambda`

,

(define fog (lambda (value) (lambda (x) (f (g x) value))))

Hum, what does it mean? Applying an argument `(fog 2)`

returns a new procedure
with one argument.

scheme@(guile-user)> (fog 2) $18 = #<procedure 7fae6c192180 at <unknown port>:542:4 (x)> scheme@(guile-user)> (define fog-2 (fog 2)) scheme@(guile-user)> (fog-2 5) $19 = 12

The module `(use-modules (srfi srfi-26))`

provides a syntax for conveniently
specializing selected parameters of a function; especially `cut`

.

(define fog (lambda (value) (compose (cut f <> value) g)))

where the slot `<>`

is filled by an argument from the created procedure.

Another example. From a number (potential floating-point), let convert to an integer. To do so, we convert the number as a string, then split the string against the dot character and last convert back the string to a number. The code reads,

(define number->integer (compose string->number car (cut string-split <> #\.) number->string))

And it works, `(number->integer 1.2)`

returns 1.

#### Lists: `map`

, `fold`

As we said earlier, the `list`

data structure plays a central role. Let
manipulate it. It requires `(use-modules (srfi srfi-1))`

.

The first essential component is `map`

. Assuming a procedure and a list, the
`map`

procedure apply the procedure to each element of the list and returns a
list. Roughly, naming `f`

the procedure and considering the list `(list 1 2 3)`

,
the result of `(map f (list 1 2 3))`

will be `(list (f 1) (f 2) (f 3))`

.

The procedure `map`

is built-in but now we have all the tools to have an
implementation, for instance, renamed `nap`

here,

(define (nap f lst) (match lst ((x xs ...) (cons (f x) (nap f xs))) (_ '())))

Give a look to this introduction about monads for some questions about composition.

Another essential component is `fold`

(also termed reduce, accumulate,
aggregate, compress, or inject). Assuming we would like to sum all the
elements of the list `(1 2 3 4 5`

), the structure of the operations reads,

+ / \ 5 + / \ 4 + / \ 3 + / \ 2 + / \ 1 0

which is exactly what `fold`

does. For instance, the code reads,

(define (sum-list lst) (fold (lambda (current previous) ;reducer (+ current previous)) 0 ;init lst)) ;list

and the result is indeed 15.

scheme@(guile-user)> (sum-list (iota 5 1)) $20 = 15

where `iota`

returns a sequence of 5 elements starting with 1.

It becomes straightforward to implement a procedure that returns ```
(-3 -2 -1 0
1 2 3)
```

. It just reads,

(define (make-real-numbers n) (let ((positives (iota (1+ n)))) (fold (lambda (element result) (cons (- element) result)) (cdr positives) positives)))

where `cons`

appends one element to an existing list. Based on that, we can
also make dual numbers.

(define (make-variables n) (let* ((varname (lambda (i) (string-append "x_" (number->string i))))) (map (lambda (i) (variable (varname i) i)) (make-numbers n))))

and the outputs looks like,

scheme@(guile-user)> ,pp (make-variables 3) $21 = (#<<dual> name: "x_-3" value: -3 deriv: 1.0> #<<dual> name: "x_-2" value: -2 deriv: 1.0> #<<dual> name: "x_-1" value: -1 deriv: 1.0> #<<dual> name: "x_0" value: 0 deriv: 1.0> #<<dual> name: "x_1" value: 1 deriv: 1.0> #<<dual> name: "x_2" value: 2 deriv: 1.0> #<<dual> name: "x_3" value: 3 deriv: 1.0>)

Awesome! One last?

scheme@(guile-user)> (for-each (compose dual-show dual-exponential2) (make-variables 2)) "~e^(x_-2)" Value: 1.3533528323660934E-1 Deriv: 1.353352832366504E-1 "~e^(x_-1)" Value: 3.6787944117144245E-1 Deriv: 3.6787944117144245E-1 "~e^(x_0)" Value: 1.0E+0 Deriv: 1.0E+0 "~e^(x_1)" Value: 2.7182818284590455E+0 Deriv: 2.7182818284590455E+0 "~e^(x_2)" Value: 7.389056098930645E+0 Deriv: 7.389056098930604E+0

with,

(define (dual-exponential2 x) (define result (let ((max-iter 20)) (fold (lambda (k result) (match result ((xk kfac expx) (let* ((kk (constant k)) (xk (dual-mul xk x)) (kfac (dual-mul kfac kk)) (term_k (dual-div xk kfac))) (list xk kfac (dual-add expx term_k)))))) (list (constant 1.) (constant 1.) (constant 1.)) (iota (1+ max-iter) 1)))) (match result ((_ _ ($ <dual> _ value deriv)) (dual (format #f "~~e^(~a)" (dual-name x)) value deriv))))

## Take the pill

**Manual**Guile « Hello Scheme »**Mini-tuto**« A Scheme Primer »**Mini-tuto**« Unlock Lisp / Scheme's magic: beginner to Scheme-in-Scheme in one hour » (video 1h)**Talk**« A tour of the Guix source tree » (video 40min)**Book**« How to Design Programs » (famous course about Racket, another Scheme; syntactic sugar varies)**Book**« Structure and Interpretation of Computer Programs (famous course from MIT)**Post**Whitespace to Lisp (wisp) « Going from Python to Guile Scheme »

Join the fun, join Guix!