Linking (lisp) libraries to FriCAS?

Alasdair

unread,

Oct 22, 2015, 11:05:18 PM10/22/15

to FriCAS - computer algebra system

I was looking at some lisp numeric libraries: gsll at https://common-lisp.net/project/gsll/ (which is a wrapper for the Gnu Scientific Library GSL), as well as mjrcalc at http://www.mitchr.me/SS/mjrcalc/. I imagined that as these were Lisp libraries, there would be some way of integrating them into FriCAS - or at least making them available so that they can be used from within a FriCAS session. However, my attempts (such as they were) came to nothing: and I wasn't sure whether to )compile, )read, or do something else.

Is there some documented standard method of integrating external libraries with FriCAS?

Thanks,
Alasdair

Waldek Hebisch

unread,

Oct 23, 2015, 3:25:02 AM10/23/15

to fricas...@googlegroups.com

Concerning documentation, I wrote some time ago about interfacing
to one of Lapack function.

Lisp functions are a bit easier than other languages, but not
much. To use a library from FriCAS, one need to first decide
in interface. That is how FriCAS types correspond to library
types. Normally this involves writing FriCAS domain/package
with apropritate declarations. One example may be plotting
subsystem: at Spad level there are viewport domains which
export somewhat high-level interface and implement operations
calling low-level functions. The low-level parts are done
by several Lisp functions which are called using '$Lisp'
qualification. The interface uses several constans defined
as Spad macros -- the constants on Spad side must agree
with constants used at other side (that is in C code).

At lower level one needs to ensure that the code is available.
Standard FriCAS code is either already loaded or FriCAS knows
how to load it on demand. With arbitrary Lisp code one
can use ')read' (which compiles files and then loads it).
I am affraid we do not have a special command to load
Lisp without compiling, be one can do this via ')lisp'
command like:

)lisp (load "prog.fasl")

Note that FriCAS may run only using AXIOMsys executable,
which contains Lisp compiler but may miss large part
of Lisp environment. For example, you had problem
because Lisp compiler inside FriCAS did not know
where to find Lisp libraries bundled with sbcl.

Currently for non-Lisp libraries one have to create
apropriate wrappers at Lisp level -- basically specify
type of C routine so that Lisp knowns how to pass
arguments. In sbcl non-Lisp dynamic libraries are
easy to use: to load them one needs to use
'|quiet_load_alien|' routine giving it path to the
library. After loadnig wrappers become usable.
When it comes to Lisp libraries each Lisp is
slightly different. In sbcl one can load Lisp
code (which either looses execution time if
file is interpreted or looses time compiling it),
load a fasl (which contains already compiled code
but still needs nontrivial time to load) or create
a new image. Images start fast, but they take
some effect to create, in particular one needs
to make sure that they are properly initialised.
Also. each image contains a copy of AXIOMsys
executable so multiple images could take
a lot of space. Concerning fasl-s, there
are some provisions to create aggregate fasl-s
but usually each file gives a sparate fasl,
so in case of multifile library one have to
load several files.

There is extra difficulty with using Lisp libraries:
modern Lisp code frequently uses so called keyword
arguments. In FriCAS using keyword arguments is
awkward, so one may be forced to create wrappers
at Lisp level which present interface free of
keyword arguments.

To summarize, if you have Lisp library in a single
file that uses FriCAS compatible data representation,
than use may be as easy as

)read "file.lisp"

DO_SOMETHING(42)$Lisp

But interesting things may require extra interface code.

BTW: Are there some specific routines that you need or
do want just to have the libraries handy?

--
Waldek Hebisch

Alasdair McAndrew

unread,

Oct 23, 2015, 5:34:49 AM10/23/15

to fricas...@googlegroups.com

Well, I was just experimenting with one of the numerical integration routines from GSLL, the Lisp wrapper library for GSL: Gnu Scientific Laboratory. In plain lisp, I can, for example enter:

*(ql:quickload "gsll")
*(gsll:integration-qng (lambda (x) (exp (- (* x x)))) 0.0 1.0)

which produces

0.7468241328124271
8.291413475940725e-15
21

Here I'm using the SLIME Emacs mode, and quicklisp to manage the Lisp libraries.

In FriCAS, I can:

)read /home/amca/quicklisp/setup.lisp

)lisp (ql:quickload "gsll")

but the command

)lisp (gsll:integration-qng (lambda (x) (exp (- (* x x)))) 0.0 1.0)

produces the error

  NIL is not a valid identifier to use in FriCAS.

which I don't understand.

I notice in Maxima (which is also based on Lisp), there is a "quadpack.lisp" file which provides an interface to the QUADPACK library.

However it seems from your detailed discussion that interfaces are in fact more difficult than I expected. Ah well. If it was easy I might just be able to manage a hacked work-around, but it looks like a job for a better programmer than me!

--

http://www.facebook.com/alasdair.mcandrew

https://plus.google.com/+AlasdairMcAndrew/posts

https://www.linkedin.com/pub/alasdair-mcandrew/a/178/108

Kurt Pagani

unread,

Oct 23, 2015, 1:40:08 PM10/23/15

to fricas...@googlegroups.com

Supplemental to Waldek's statements you might have a look at the
following examples where some nontrivial Lisp packages have been
integrated to FricAS:

https://bitbucket.org/kfp/snark

Espsecially the file snark.lisp and snark.spad.

Also https://bitbucket.org/kfp/fricas_ajax
and https://bitbucket.org/kfp/fricas_wax are examples for how to easily
include some libraries (Hunchentoot webserver in this case).

Last but not least there is now a binary version of iSPAD for the
FriCAS/SBCL where you can see how to link AxiomSYS to your favourite
library (see install.sh and sbcl.lisp).
https://github.com/nilqed/fricas_jupyter
https://hub.docker.com/r/nilqed/fricas_jupyter/

Kurt

Kurt Pagani

unread,

Oct 23, 2015, 7:03:11 PM10/23/15

to fricas...@googlegroups.com

Hello Alasdair

What I have completely forgotten to mention in the last post:

In the fricas/contrib folder there is a file load-fricas.lisp which
allows you to use FriCAS as a cl library itself. This way you can mix in
any kind of library as long there is no interference to the "boot"
package (I'm aware of only very few ones).

BTW: I found that
(load "../lisp/primitives.lisp")
is missing in the file, so if you encounter problems it might be this
what caused an error. (see below). I used it occasionally in the
cl-jupyter frontend.

--
Kurt

#|
To use FriCAS as a Lisp library build FriCAS as usual and
keep the build tree. Change path below to point to build
tree. After that doing:

(load "load-fricas.lisp")

will load and initialize FriCAS.

Examples:
Using FriCAS command interpreter:

(in-package "BOOT")
(|parseAndInterpret| "x^2")

2
(2) x
Type: Polynomial Integer
((Polynomial (Integer)) WRAPPED 1 x (2 0 . 1))

Directly calling math functions:

(defvar *mult-i-sup*
(|getFunctionFromDomain|
'*
'(|SparseUnivariatePolynomial| (|Integer|))
'((|SparseUnivariatePolynomial| (|Integer|))
(|SparseUnivariatePolynomial| (|Integer|)))))
*MULT-I-SUP*

(SPADCALL '((2 . 1) (0 . 1)) '((5 . 1) (0 . 1)) *mult-i-sup*)
((7 . 1) (5 . 1) (2 . 1) (0 . 1))

Notes:
- all intersting functionality is in package called "BOOT".
- at mathematical level FriCAS is case-sensitive, so at Lisp level one
has to use bars.
- the simplest interface is |parseAndInterpret| which takes a string
as input and produces a Lisp form repesenting printed output. As
side effect |parseAndInterpret| prints the result.
- at deeper lever FriCAS functions are overloaded, so to call correct
function one has to first use |getFunctionFromDomain| to get
function which matches to given argument types. Above I want to
multiply two sparse univarate polynomials with integer coefficients.
Since lookup may be expensive the caller is adviced to cache result
of the lookup.
- FriCAS functions use special calling convention, so one has to use
SPADCALL macro to call them. Actually, |getFunctionFromDomain|
returns a pair consistion of a function and an extra argument.
SPADCALL takes care of decomposing the pair and appending the
extra argument to the argument list.

Currently FriCAS sets a few system (global) variables, for example
*read-default-float-format* is set to 'double-float -- in principle
FriCAS settings may interfere with other programs.

|#
(let ((*default-pathname-defaults*
#P"~/Development/ax-build/src/interp/"))
(load "../lisp/fricas-package.lisp")
(load "../lisp/fricas-config.lisp")
(load "../lisp/fricas-lisp")
(load "../lisp/primitives.lisp") ;;;+kfp
(load "makeint.lisp"))
(in-package "BOOT")
(fricas-init)

(in-package "BOOT")
(|parseAndInterpret| "x^2")
kfp@helix:~/Development/fricas/contrib$

Am 23.10.2015 um 09:24 schrieb Waldek Hebisch:

Alasdair McAndrew

unread,

Oct 23, 2015, 9:13:36 PM10/23/15

to fricas...@googlegroups.com

Just as an experiment I'm taking the quadpack.lisp file from the Maxima distribution, to see if I can install it and have it working under FriCAS. However, as a total Lisp newbie, I'm continuously stumped by trivial errors. For example, quadpack.lisp uses the f2cl-lib package. I can install this with quicklisp, but I don't know (yet) how to make that package available from within FriCAS. For that matter, I don't know how to make the quicklisp package ql available: even the line

(load "~/quicklisp/setup.lisp")

which should work, doesn't appear to help - lines beginning with "(ql: " are immediately flagged as errors and the ")read" compile is aborted.

I'll keep on fiddling... (while Rome burns).

-Alasdair

Kurt Pagani

unread,

Oct 24, 2015, 12:11:27 AM10/24/15

to fricas...@googlegroups.com

I'm trying to explain it as good as I can (experts will correct me):

I see three different methods:

1. start your Lisp REPL (pref. SBCL) and try to load QuickLisp, FriCAS:

- (load "~/quicklisp/setup.lisp")
- (load "load-fricas.lisp")
- (ql:quickload ... whatever you want to via QuickLisp)
- (load "other_lisp_files_of_interest")
- (main)

where "main" is a function which combines the various objects. You said
you were "a total Lisp newbie", then this method might not be the first
choice.

2. Instead of using a Lisp REPL one can use the one inside FriCAS,
however, except from loading "load-fricas.lisp" nothing changes with
respect to method 1, at least as long as you will stay within the Lisp
REPL.

$ fricas -nosman
)fin
...
|spad| ; back to fricas

3. Write a Lisp file containing all the "load" commands (except
"load-fricas.lisp" of course) and any packages and functions necessary
to "export" to the Fricas shell. Afterwards you have to write a SPAD
package which "imports" your Lisp functions (those you want to have in
FriCAS). I've found an example for QuickLisp (it's rather unfinished, so
I don't recommend to compile it): the SPAD section below is the file
"clx.spad" while in the second section below (after LISP) you can see
the file clx.lisp.

Now compare for instance the function qlQuickLoad in clx.spad:

--
qlQuickLoad : String -> SExpression
++ qlQuickLoad is Lisp ql:quickload
...
qlQuickLoad(s) == QL_-QUICKLOAD(s)$Lisp
--

with the function ql-quickload in clx.lisp:

;;; =================
;;; QuickLisp section
;;; =================

(defun ql-quickload (s)
(if (find-package "QUICKLISP")
(ql:quickload s)
(clx-error *pkg-err-ql*)))

Now freely interpreted we have the following situation:

If you call (in FriCAS) the functtion

(1)-> qlQuickLoad("myLib")

then in the underground (Lisp) the expression

(ql-quickload "myLib") is evaluated.

Note that Common-Lisp is case insensitive (unless you use |...|), so
that in SPAD the function "ql-quickload" should be written as
QL_-QUICKLOAD, where the underscore is necessary to escape the special
character "-".

I admit it's confusing at times. But comparing corresponding functions
in the different sections you will get a feeling for it.

In order to use "quadpack.lisp" within FriCAS you have just to do the
same as described above; choose the functions you want to see in FriCAS
and map them to the corresponding Lisp expressions in a SPAD package (or
domain). Of course, Maxima and FriCAS objects may be slightly different
internally, but in this case I think "quadpack" is quite generic since
it is a translation from Fortran (as you mentioned) and the
implementation should be pretty straighforward.

A last remark: the Maxima function "quad_qag" for instance should be
named as "quadQag" (recommended but not mandatory ;)

Kurt

**************************

----------
-- SPAD --
----------

$ cat clx.spad
)lisp (load "clx")
)abbrev package CLX clx
++ Date Created: Sun Dec 21 03:29:45 CET 2014
++ License: BSD (same as Axiom)
++ Date Last Updated:
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:
++
++ Description:
++ CLX is a support package for Quicklisp, MathJax, tex2png and some
++ other utilities.
++
clx() : Exports == Implementation where

Exports == with

%Exec : String -> SExpression
++ \spad{\%Exec} executes a command s in the shell
%writeMathJax : String -> SExpression
++ \spad{\%writeMathJax} write TeX string to /tmp/mjax.html
%displayMathJax : String -> SExpression
++ \spad{\%displayMathJax} write&display TeX string to /tmp/mjax.html

qlQuickLoad : String -> SExpression
++ qlQuickLoad is Lisp ql:quickload
qlUnInstall : String -> SExpression
++ qlUnInstall is Lisp ql:uninstall
qlSystemApropos : String -> SExpression
++ qlSystemApropos is Lisp ql:system-apropos
qlUpdateAllDists : () -> SExpression
++ qlUpdateAllDists is Lisp ql:update-all-dists
qlUpdateClient : () -> SExpression
++ qlUpdateClient is Lisp ql:update-client
qlWhoDependsOn : String -> SExpression
++ qlWhoDependsOn is Lisp ql:who-depends-on

mjax : Any -> SExpression
++ mjax display any expression (if possible) in the browser

evalToSlide : String -> SExpression
evalToString : String -> SExpression

Implementation == add

%Exec(s) == CLX_-EXEC(s)$Lisp
%writeMathJax(s) == CLX_-WRITE_-MATHJAX(s)$Lisp
%displayMathJax(s) == CLX_-DISPLAY_-MATHJAX(s)$Lisp

qlQuickLoad(s) == QL_-QUICKLOAD(s)$Lisp
qlUnInstall(s) == QL_-UNINSTALL(s)$Lisp
qlSystemApropos(s) == QL_-SYSTEM_-APROPOS(s)$Lisp
qlUpdateAllDists() == QL_-UPDATE_-ALL_-DISTS()$Lisp
qlUpdateClient() == QL_-UPDATE_-CLIENT()$Lisp
qlWhoDependsOn(s) == QL_-WHO_-DEPENDS_-ON(s)$Lisp

mjax(x:Any):SExpression ==
t:String:=first tex(x::OutputForm::TexFormat)
s:String:= concat ["$$", t::String, "$$"]
CLX_-DISPLAY_-MATHJAX(s)$Lisp

--evalToString(s) == CLX_-EVAL_-TO_-STRING(s)$Lisp -- blocking ???
--evalToSlide(s) == CLX_-EVAL_-TO_-SLIDE(s)$Lisp
--r:String:=
--split(r,char "@")

;;;;;;;;;;;;
;;; LISP ;;;
;;;;;;;;;;;;

$ cat clx.lisp
(in-package :boot)

(defconstant *pkg-err-asdf* "ASDF package missing.")
(defconstant *pkg-err-ql* "QuickLisp package missing.")

(defun clx-error (s)
"Print error message"
(print s))

(defun clx-exec (s)
"Function: run-shell-command control-string &rest args [oboslete]"
(if (find-package "ASDF")
(asdf:run-shell-command s)
(clx-error *pkg-err-asdf*)))

(defun clx-concstr (list)
;; concatenate a list of strings ; recall ~% = newline"
(if (listp list)
(with-output-to-string (s)
(dolist (item list)
(if (stringp item)
(format s "~a~%" item))))))

(defun clx-getenv (name &optional default)
"http://cl-cookbook.sourceforge.net/os.html; get env variables"
;; e.g. (clx-getenv "OS") => "Windows_NT"
#+CMU
(let ((x (assoc name ext:*environment-list* :test #'string=)))
(if x (cdr x) default))
#-CMU
(or
#+Allegro (sys:getenv name)
#+CLISP (ext:getenv name)
#+ECL (si:getenv name)
#+SBCL (sb-unix::posix-getenv name)
#+LISPWORKS (lispworks:environment-variable name)
default))

;;; =================
;;; QuickLisp section
;;; =================

(defun ql-quickload (s)
(if (find-package "QUICKLISP")
(ql:quickload s)
(clx-error *pkg-err-ql*)))

(defun ql-uninstall (s)
(if (find-package "QUICKLISP")
(ql:uninstall s)
(clx-error *pkg-err-ql*)))

(defun ql-system-apropos (s)
(if (find-package "QUICKLISP")
(ql:system-apropos s)
(clx-error *pkg-err-ql*)))

(defun ql-update-all-dists ()
(if (find-package "QUICKLISP")
(ql:update-all-dists)
(clx-error *pkg-err-ql*)))

(defun ql-update-client ()
(if (find-package "QUICKLISP")
(ql:update-client)
(clx-error *pkg-err-ql*)))

(defun ql-who-depends-on (s)
(if (find-package "QUICKLISP")
(ql:who-depends-on s)
(clx-error *pkg-err-ql*)))

;;; ===============
;;; MathJax section
;;; ===============

(defconstant *url-show* "firefox -url file://C:/cygwin/tmp/~A")

(defconstant *mjax-template* "<!DOCTYPE html>
<html>
<head>
<title>FriCAS Output</title>
<script type='text/x-mathjax-config'>
MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\$','\$']]}});
</script>
<script type='text/javascript'

src='http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'>
</script>
</head>
<body>
~A
</body>
</html>")

(defun clx-write-mathjax (s)
(with-open-file (*standard-output* "/tmp/mjax.html"
:direction :output
:if-exists :supersede)
(format t *mjax-template* s)))

(defun clx-display-mathjax (s)
(progn (clx-write-mathjax s)
(clx-exec (format nil *url-show* "mjax.html"))))

;;; =============
;;; LaTeX section
;;; =============

(defconstant *do-latex* "latex")

(defconstant *beamer-template*
"\\documentclass[pdf]{beamer}
%\\usetheme{Warsaw}
\\usetheme{PaloAlto}
\\title[]{FriCAS}
\\subtitle[subttitle]{eval}
\\author[]{FriCAS}
\\institute[]{Cygwin}
\\date{\\today}
\\begin{document}
\\begin{frame}
\\titlepage
\\end{frame}
~A
\\end{document}")

(defconstant *frame-template* "\\begin{frame} $$~A$$ \\end{frame}")

(defun clx-eval-to-string (s) (|parseAndEvalToString| s))
(defun clx-eval-to-latex (s) (clx-eval-to-string (format nil "tex(~A)" s)))

(defun clx-make-slide (s)
(format nil *frame-template* (car (clx-eval-to-latex s))))

(defun clx-write-slide (s)
(with-open-file (*standard-output* "/tmp/slide.tex"
:direction :output
:if-exists :supersede)
(format t *beamer-template* (clx-make-slide s))))

(defun clx-display-slide (s)
(progn (clx-write-slide s)
(clx-exec (format nil *do-latex* "slide.tex"))))

Alasdair McAndrew

unread,

Oct 24, 2015, 4:28:09 AM10/24/15

to fricas...@googlegroups.com

The more I read about Lisp the more difficult it seems. It would probably be just as easy for me to write a integrate.input file from scratch containing the Gauss-Kronrod code, rather than wrestling with a language of which I have little understanding, and no competence!

Bill Page

unread,

Oct 24, 2015, 10:28:29 AM10/24/15

to fricas-devel

Alasdair,

For entirely selfish reasons :) I would like to encourage you to
continue work on this since I also would like to be able use some of
the standard numeric libraries in FriCAS. But I am also sympathetic to
your reaction concerning the details in Lisp. One thing that might
help is that both FriCAS and OpenAxiom follow the original design of
Axiom which includes an intermediate language called BOOT. BOOT was
intended to make Lisp-level programming easier for people who mostly
program in SPAD since the syntax is similar to SPAD without static
types. Some (most?) of what Kurt wrote in Lisp below and all of the
interface code that he did not include could be written in a
combination of BOOT and SPAD.

In any case I am motivated to try to help with this.

Cheers,
Bill Page

> You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.

Kurt Pagani

unread,

Oct 24, 2015, 5:16:34 PM10/24/15

to fricas...@googlegroups.com

I'm endorsing what Bill suggested. If you are going to manage just one
function the rest of it will be routine.

This might interest you:
https://github.com/matlisp/matlisp

Especially
https://github.com/matlisp/matlisp/blob/t2/src/packages/quadpack/quadpack.lisp
which seems to me quite easier to implement.

Alasdair McAndrew

unread,

Oct 24, 2015, 9:22:06 PM10/24/15

to fricas...@googlegroups.com

I'll give it my best shot... time permitting! Here are the issues I face so far:

QUADPACK is designed for use with IEEE double precision arithmetic. Thus the nodes and weights for numeric integration can be given as fixed values. However, for FriCAS it would be proper, I think, to allow arbitrary precision. The Gaussian nodes can be obtained as the roots of Legendre polynomials, but the interpolating Kronrod nodes are more difficult; one would have to use a method first developed by Dirk Laurie: http://www.ams.org/mcom/1997-66-219/S0025-5718-97-00861-2/S0025-5718-97-00861-2.pdf and since extended by Calvetti et al: http://www.ams.org/mcom/2000-69-231/S0025-5718-00-01174-1/S0025-5718-00-01174-1.pdf

Many computations require eigenvalues of floating point matrices. We could probably do this in FriCAS by solving the characteristic equation - I don't think FriCAS has a inbuilt method for floating point eigenvalues. However when I tried to create a characteristic equation just now FriCAS complained that

"Union(Polynomial(Float),"failed") cannot be coerced to mode Polynomial(Float)"

The QUADPACK routines were initially written in Fortran, then ported to C for use in GSL, as well as to Lisp (GSLL), Python (Numpy), and lots of other languages. So it's a matter of picking one which is easiest to work with!

Maxima (which is also written in Lisp) has a QUADPACK wrapper: there is a function quadpack.lisp which allows access from with Maxima to all the QUADPACK routines. Unfortunately you can't browse Maxima source online, so I've attached it here for your inspection. This file requires access to other files, just as the numeric routines of QUADPACK itself make use of various "include" files (if in C) and associated libraries.

It is probably just a matter, in the first instance, of putting everything where FriCAS can find it, and stitching everything together so that all the necessary libraries can be found as needed. However, the devil is in the details, and my knowledge of how Lisp works could be engraved in 10pt bold font on the head of a small pin.

I would be happy to start by being able to enter, in FriCAS, for example:

)lisp (gsll:integration-qng (lambda (x) (exp (- (* x x)))) 0.0 1.0)

but at the moment, although I can run that command from with SBCL, I can't run it from FriCAS - there's a Lisp problem in making quicklisp libraries and GSLL available from within FriCAS. This is where my ignorance becomes an obstacle.

If I was to choose one function it would be one of QNG (non-adaptive integrator using Gauss-Kronrod-Peterson rules of increasing precision), or QAG (an adaptive integrator using Gauss-Kronrod integration).

Righto, back to real work of collating student results...

quadpack.lisp

Waldek Hebisch

unread,

Oct 24, 2015, 11:00:54 PM10/24/15

to fricas...@googlegroups.com

Alasdair McAndrew wrote:
>
> The more I read about Lisp the more difficult it seems. It would probably
> be just as easy for me to write a integrate.input file from scratch
> containing the Gauss-Kronrod code,

Any reason for specifically wanting Gauss-Kronrod code? On
my todo list was writing adaptive Gauss routine. So I did
one using 9 point rule. Hopefuly it should be good enough
to get 14 digit accuracy with reasonable number of evaluation
points for large class of functions. This is preliminary
version, improvements welcome.

BTW: I believe that adaptive quadratures, recursing on
subintervals are much better than nonadaptive. This one
is quite simplistic, in particular it can handle singularity
in sqrt, but fails on 1/sqrt.

-------------<cut here>-----------------

)abbrev package GAUSSI GaussianQuadrature
GaussianQuadrature : Exports == Implementation where
FT ==> DoubleFloat
Exports ==> with
adaptive : (FT -> FT, FT, FT, FT,
Integer, (FT -> FT, FT, FT) -> FT) -> FT
++ adaptive(f, a, b, eps, n, g) is general adaptive
++ quadrature on interval (a, b) using g as base
++ (nonadaptive) rule and stopping when error is less than
++ eps. adaptive signals error when number of recursion
++ levels exceeds n.
++ Note: agruments to g are f, midpoint of interval and
++ half of length
gauss9 : (FT -> FT, FT, FT) -> FT
++ gauss9(f, mid, h) computes 0 point Gaussian quadrature
++ of f on the interval with midpoint mid and lenght 2*h
Implementation ==> add

x1 := 0.3242534234_0380892904::DoubleFloat
x2 := 0.6133714327_0059039731::DoubleFloat
x3 := 0.8360311073_266357943::DoubleFloat
x4 := 0.9681602395_0762608983::DoubleFloat
w0 := 0.3302393550_0125976316::DoubleFloat
w1 := 0.3123470770_4000284007::DoubleFloat
w2 := 0.2606106964_0293546232::DoubleFloat
w3 := 0.1806481606_9485740405_8::DoubleFloat
w4 := 0.0812743883_6157441197_2::DoubleFloat
half := 0.5::DoubleFloat

gauss9(f : FT -> FT, mid : FT, h : FT) : FT ==
h1 := h*x1
h2 := h*x2
res := w0*f(mid) + w1*(f(mid + h1) + f(mid - h1))
+ w2*(f(mid + h2) + f(mid - h2))
h3 := h*x3
h4 := h*x4
res := res + w3*(f(mid + h3) + f(mid - h3))
+ w4*(f(mid + h4) + f(mid - h4))
h*res

adaptive1(f : FT -> FT, mid : FT, h : FT, eps : FT,
n : Integer, g : (FT -> FT, FT, FT) -> FT, res0 : FT) : FT ==
n < 0 => error "too many recursion levels"
h1 := half*h
mid1 := mid - h1
mid2 := mid + h1
res1 := g(f, mid1, h1)
res2 := g(f, mid2, h1)
nres := res1 + res2
abs(res0 - nres) < eps => nres
eps1 := half*eps
adaptive1(f, mid1, h1, eps1, n - 1, g, res1)
+ adaptive1(f, mid2, h1, eps1, n - 1, g, res2)

adaptive(f : FT -> FT, a : FT, b : FT, eps : FT,
n : Integer, g : (FT -> FT, FT, FT) -> FT) : FT ==
mid := half*(a + b)
h := half*(b - a)
res0 := g(f, mid, h)
adaptive1(f, mid, h, eps, n, g, res0)

--------------------<cut here>-------------------

--
Waldek Hebisch

Bill Page

unread,

Oct 24, 2015, 11:26:42 PM10/24/15

to fricas-devel

Alasdair,

Have you been successful running the example below using SBCL? In my
case I get an error during quickload:

[package osicat]..................................
[package gsll].;
; caught ERROR:
; READ error during COMPILE-FILE:
;
; end of file on #<SB-IMPL::STRING-INPUT-STREAM {100AF084C3}>
;
; (in form starting at line: 68, column: 0, position: 2729)

--

My SuSE Linux SLES 11/SP3 installation is now a bit old and before
getting even this far I first had to recompile libffi from source.
Maybe something else is wrong but I am wondering if this works with
SBCL/

Thanks.

Bill Page.

Alfredo Portes

unread,

Oct 25, 2015, 2:50:18 AM10/25/15

to fricas-devel

It would be a nice addition to Fricas to include the FFI system that Gaby implemented in OpenAxiom.

https://parasol.tamu.edu/~gdr/OpenAxiom/ffi-draft.pdf

Alasdair McAndrew

unread,

Oct 25, 2015, 3:56:58 AM10/25/15

to fricas...@googlegroups.com

I have no particular attachment to Gauss-Kronrod integration - indeed Clenshaw-Curtis integration may well be easier to implement. My interest (as well as having a robust numeric integration routine in FriCAS) is seeing if an external numeric library can be used through FriCAS. And there are lots of these: QUADPACK, FFTW, ODE, ODEPACK to name but four Netlib ones. I don't see why any of us should have to reinvent the wheel when there are perfectly good open-source numeric routines already available. All I want is a FriCAS front end to them - and with arbitrary precision.

For example, the nodes and weights above could be obtained with:

digits(40)

outputGeneral(40)
x1 := solve(legendreP(9,x),1.0e-40)

xs := [rhs(s) for s in x1]
ws := [2/(1-s^2)/subst(D(legendreP(9,x),x),x=s)^2 for s in xs]

And then you have to compute the Kronrod nodes, which interpolate between the Gaussian nodes...

Waldek Hebisch

unread,

Oct 25, 2015, 8:06:39 AM10/25/15

to fricas...@googlegroups.com

Alasdair McAndrew wrote:
>
> I have no particular attachment to Gauss-Kronrod integration - indeed
> Clenshaw-Curtis integration may well be easier to implement. My interest
> (as well as having a robust numeric integration routine in FriCAS) is
> seeing if an external numeric library can be used through FriCAS. And
> there are lots of these: QUADPACK, FFTW, ODE, ODEPACK to name but four
> Netlib ones. I don't see why any of us should have to reinvent the wheel
> when there are perfectly good open-source numeric routines already
> available. All I want is a FriCAS front end to them - and with arbitrary
> precision.

Well, we can use routines from other projects if they are
available. However, note that popular libraries typically
are limited to double precision. To get arbitrary precision
in most cases we will have to roll our own. The routine
I wrote is limited to double precision, but once we work
out algorithms we should be able to produce version for
arbitarary precision (as you note below we can produce
nodes and weights for Gaussian quadrature of arbitrary
order). In fact, from my point of view main reason
for using foreign libraries is speed -- assuming that
authors put more effort into speed optimizations that
we want to. However quadratures (and to some degree
ODE solvers) are special -- their execution time
should be dominated by user provided routines. Assuming
that we use good algorithm we should be at least as
fast as external libraries. Namely with good algorithm
we will have the number of calls to user routines and
calling FriCAS routine (I assume that user routine is
written in FriCAS) from FriCAS is faster than calling
it from outside of FriCAS.

In case of elliptic functions I do not know of any
open source system having routines as good as ours.
OTOH Lapack and Blas are definitely worth reusing
(we still should roll our owen arbitrary precision
versions but for machine precision Lapeck and Blas
have definite speed advantage).

So part important part of work is to identify good
routines -- I hope on some first hand recommendations
for good routines.

>
> For example, the nodes and weights above could be obtained with:
>
> digits(40)
> outputGeneral(40)
> x1 := solve(legendreP(9,x),1.0e-40)
> xs := [rhs(s) for s in x1]
> ws := [2/(1-s^2)/subst(D(legendreP(9,x),x),x=s)^2 for s in xs]
>
> And then you have to compute the Kronrod nodes, which interpolate between
> the Gaussian nodes...

--
Waldek Hebisch

Kurt Pagani

unread,

Oct 25, 2015, 6:16:53 PM10/25/15

to fricas...@googlegroups.com

Hello Bill

I got the same error once; after

sudo apt-get install libgsl0-dev

it works fine.

Kurt

It might also be a path problem (not very probable):
http://stackoverflow.com/questions/28692610/install-gsll-on-sbcl-with-quicklisp

Kurt Pagani

unread,

Oct 25, 2015, 6:25:40 PM10/25/15

to fricas...@googlegroups.com

Hello Alasdair

Have you tried the following?
---

kfp@helix:~$ fricas -nox
viewman not present, disabling graphics
Checking for foreign routines
AXIOM="/usr/local/lib/fricas/target/x86_64-unknown-linux"
spad-lib="/usr/local/lib/fricas/target/x86_64-unknown-linux/lib/libspad.so"
foreign routines found
openServer result 0
FriCAS Computer Algebra System
Version: FriCAS 2014-12-18
Timestamp: Son Okt 25 19:16:29 CET 2015
-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave FriCAS and return to shell.
-----------------------------------------------------------------------------

(1) -> )fin
* (load "~/quicklisp/setup")

T
* (ql:quickload "gsll")
To load "gsll":
Load 1 ASDF system:
gsll
; Loading "gsll"
.............
("gsll")
* (gsll:integration-qng (lambda (x) (exp (- (* x x)))) 0.0 1.0)

0.7468241328124271
8.291413475940725e-15
21
* (|spad|)

(1) ->
Value = 0.7468241328124271
(1) -> )lisp (gsll:integration-qng (lambda (x) (exp (- (* x x)))) 0.0 1.0)

Value = 0.7468241328124271
(1) ->

---

If the above works then next you want to write a SPAD function for
"gsll:integration-qng". Why we get only one value returned is puzzling,
but I'll have a look into the code.

Best regards
Kurt

Waldek Hebisch

unread,

Oct 25, 2015, 7:12:25 PM10/25/15

to fricas...@googlegroups.com

Kurt Pagani wrote:
>
> Value = 0.7468241328124271
> (1) -> )lisp (gsll:integration-qng (lambda (x) (exp (- (* x x)))) 0.0 1.0)
>
> Value = 0.7468241328124271
> (1) ->
>
> ---
>
> If the above works then next you want to write a SPAD function for
> "gsll:integration-qng". Why we get only one value returned is puzzling,
> but I'll have a look into the code.

FriCAS does not support multiple values. We could modify ')lisp'
to print additional values, but '$Lisp' calls have no place
where to put extra values, so one can not access then from
Spad or interpreter.

--
Waldek Hebisch

Bill Page

unread,

Oct 25, 2015, 8:11:56 PM10/25/15

to fricas-devel

Kurt,

That did the trick! Thanks. It's kind of disappointing that QuickLisp
doesn't do a better job of informing about missing dependencies. I
wonder why one needs the -dev version?

Oh well.

Bill.

> --

> You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.

Bill Page

unread,

Oct 25, 2015, 9:51:03 PM10/25/15

to fricas-devel

Kurt,

When I try what you wrote I get the following:

wspage@suse:~> fricas -nox
Checking for foreign routines
AXIOM="/usr/local/lib/fricas/target/x86_64-suse-linux"
spad-lib="/usr/local/lib/fricas/target/x86_64-suse-linux/lib/libspad.so"

foreign routines found
openServer result 0
FriCAS Computer Algebra System
Version: FriCAS 2014-12-18

Timestamp: Sun Oct 25 21:18:13 EDT 2015

-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave FriCAS and return to shell.
-----------------------------------------------------------------------------

(1) -> )fin
* (load "~/quicklisp/setup")

While evaluating the form starting at line 129, column 0
of #P"~/quicklisp/setup.lisp":

debugger invoked on a SB-INT:EXTENSION-FAILURE in thread
#<THREAD "main thread" RUNNING {100450CC23}>:
Don't know how to REQUIRE SB-BSD-SOCKETS.
See also:
The SBCL Manual, Variable SB-EXT:*MODULE-PROVIDER-FUNCTIONS*
The SBCL Manual, Function REQUIRE

Type HELP for debugger help, or (SB-EXT:EXIT) to exit from SBCL.

restarts (invokable by number or by possibly-abbreviated name):
0: [TRY-RECOMPILING] Recompile impl and try loading it again
1: [RETRY ] Retry
loading FASL for #<CL-SOURCE-FILE "quicklisp" "impl">.
2: [ACCEPT ] Continue, treating
loading FASL for #<CL-SOURCE-FILE "quicklisp" "impl"> as
having been successful.
3: [RETRY ] Retry EVAL of current toplevel form.
4: [CONTINUE ] Ignore error and continue loading file
"/home/wspage/quicklisp/setup.lisp".
5: [ABORT ] Abort loading file "/home/wspage/quicklisp/setup.lisp".
6: Exit debugger, returning to top level.

(SB-IMPL::REQUIRE-ERROR "Don't know how to ~S ~A." REQUIRE SB-BSD-SOCKETS)
0]

---

I am not sure what is going on - maybe in the wrong default package?

The following recipe does work for me however:

wspage@suse:~> sbcl
This is SBCL 1.2.16, an implementation of ANSI Common Lisp.
More information about SBCL is available at <http://www.sbcl.org/>.

SBCL is free software, provided as is, with absolutely no warranty.
It is mostly in the public domain; some portions are provided under
BSD-style licenses. See the CREDITS and COPYING files in the
distribution for more information.

* (load "~/quicklisp/setup")

T

* (ql:quickload "gsll")
To load "gsll":
Load 1 ASDF system:
gsll
; Loading "gsll"

..............
("gsll")
* (load "load-fricas")

; file: /home/wspage/fricas-build/src/lisp/fricas-package.lisp
; in: DEFMACRO IN-PACKAGE
; (DEFMACRO FRICAS-LISP::IN-PACKAGE (PACKAGE &REST FRICAS-LISP::OPTIONS)
; `(IN-PACKAGE ,PACKAGE))
; --> SB-C::NAMED-DS-BIND SB-INT:BINDING*
; ==>
; (LET* ((#:G0
; (SB-C::CHECK-DS-LIST/&REST (CDR #:EXPR) 1 1
; '(# PACKAGE &REST FRICAS-LISP::OPTIONS)))
; (PACKAGE (POP #:G0))
; (FRICAS-LISP::OPTIONS #:G0))
; (BLOCK FRICAS-LISP::IN-PACKAGE `(IN-PACKAGE ,PACKAGE)))
;
; caught STYLE-WARNING:
; The variable OPTIONS is defined but never used.
;
; compilation unit finished
; caught 1 STYLE-WARNING condition
WARNING: Reference to deprecated function (SB-EXT:QUIT) from EXIT-WITH-STATUS
STYLE-WARNING: Undefined alien: "writeablep"
STYLE-WARNING: Undefined alien: "remove_directory"
STYLE-WARNING: Undefined alien: "open_server"
STYLE-WARNING: Undefined alien: "sock_get_int"
STYLE-WARNING: Undefined alien: "sock_send_int"
STYLE-WARNING: Undefined alien: "sock_get_float"
STYLE-WARNING: Undefined alien: "sock_send_float"
STYLE-WARNING: Undefined alien: "sock_send_string_len"
STYLE-WARNING: Undefined alien: "server_switch"
STYLE-WARNING: Undefined alien: "sock_send_signal"
STYLE-WARNING: Undefined alien: "sock_get_string_buf"
Re-reading compress.daase Re-reading interp.daase
FriCAS initialization: interpreter
FriCAS initialization: database
FriCAS initialization: constructors
FriCAS initialization: history

FriCAS Computer Algebra System
Version: FriCAS 2014-12-18

Timestamp: Sun Oct 25 21:18:13 EDT 2015

-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave FriCAS and return to shell.
-----------------------------------------------------------------------------

Re-reading compress.daase Re-reading interp.daase
Re-reading operation.daase
Re-reading browse.daase
Re-reading category.daase
Initial getdatabase
preloading /algebra/FR.fasl..skipped.
...

Checking for foreign routines
AXIOM=NIL
spad-lib="/lib/libspad.so"

FriCAS Computer Algebra System
Version: FriCAS 2014-12-18

Timestamp: Sun Oct 25 21:18:13 EDT 2015

-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave FriCAS and return to shell.
-----------------------------------------------------------------------------

T
* (in-package "BOOT")

#<PACKAGE "BOOT">
* (|spad|)

(1) -> )lisp (gsll:integration-qng (lambda (x) (exp (- (* x x)))) 0.0 1.0)

Value = 0.7468241328124271

---

I guess I can stick this stuff into '.sbclrc' and end up in FriCAS
that way, but would there be an easier and more portable way to
accomplish this if this should ever become a standard part of the
FriCAS distribution? Can the ql:quickload stuff be done during make
and become part of the FriCAS image?

Cheers,
Bill Page.

On 25 October 2015 at 18:25, Kurt Pagani <nil...@gmail.com> wrote:
> Hello Alasdair
>
>

Kurt Pagani

unread,

Oct 25, 2015, 9:55:21 PM10/25/15

to fricas...@googlegroups.com

> Don't know how to REQUIRE SB-BSD-SOCKETS

usually indicates that SBCL_HOME is not set.

Did you?

Bill Page

unread,

Oct 25, 2015, 10:23:02 PM10/25/15

to fricas-devel

Ah, OK cool.

So next: Write a small Lisp wrapper function to call the the foreign
routine and return it's values as a list or maybe an array?

---

wspage@suse:~> export SBCL_HOME=/usr/local/lib/sbcl

wspage@suse:~> fricas -nox
Checking for foreign routines
AXIOM="/usr/local/lib/fricas/target/x86_64-suse-linux"
spad-lib="/usr/local/lib/fricas/target/x86_64-suse-linux/lib/libspad.so"
foreign routines found
openServer result 0
FriCAS Computer Algebra System
Version: FriCAS 2014-12-18
Timestamp: Sun Oct 25 21:18:13 EDT 2015
-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave FriCAS and return to shell.
-----------------------------------------------------------------------------

(1) -> )fin
* (load "~/quicklisp/setup")

T
* (ql:quickload "gsll")
To load "gsll":
Load 1 ASDF system:
gsll
; Loading "gsll"
.............
("gsll")

* (|spad|)
(1) -> )lisp (gsll:integration-qng (lambda (x) (exp (- (* x x)))) 0.0 1.0)

Value = 0.7468241328124271

(1) ->

---

Kurt Pagani

unread,

Oct 25, 2015, 10:43:28 PM10/25/15

to fricas...@googlegroups.com

As you commanded (files attached) :)

What I do not understand yet: what's the meaning of the other two
values? Did you already look into the gsll sources?

*(gsll:integration-qng (lambda (x) (exp (- (* x x)))) 0.0 1.0)

0.7468241328124271
8.291413475940725e-15
21

kfp@helix:~/Development/GSL$ fricas -nox

viewman not present, disabling graphics

Checking for foreign routines
AXIOM="/usr/local/lib/fricas/target/x86_64-unknown-linux"
spad-lib="/usr/local/lib/fricas/target/x86_64-unknown-linux/lib/libspad.so"

foreign routines found
openServer result 0
FriCAS Computer Algebra System
Version: FriCAS 2014-12-18

Timestamp: Son Okt 25 19:16:29 CET 2015

-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave FriCAS and return to shell.
-----------------------------------------------------------------------------

(1) -> )co gsl
Compiling FriCAS source code from file
/home/kfp/Development/GSL/gsl.spad using old system compiler.
Value = T

To load "gsll":
Load 1 ASDF system:
gsll
; Loading "gsll"
.............

;
; compilation unit finished
; caught 1 STYLE-WARNING condition

; /home/kfp/Development/GSL/GSL.NRLIB/GSL.fasl written
; compilation finished in 0:00:00.054
------------------------------------------------------------------------
gsl is now explicitly exposed in frame frame1
gsl will be automatically loaded when needed from
/home/kfp/Development/GSL/GSL.NRLIB/GSL

(1) -> DF ==> DoubleFloat
f:=(x:DF):DF+->exp(-x^2)
g(x:DF):DF == exp(-x^2)
gslIntegrationQng(f,0.0::DF,1.0::DF)
Type:
Void
(2) ->
(2) theMap(*1;anonymousFunction;0;frame1;internal)
Type: (DoubleFloat ->
DoubleFloat)
(3) -> Function declaration g : DoubleFloat -> DoubleFloat has been added
to workspace.
Type:
Void
(4) ->
(4) 0.7468241328124271
Type:
DoubleFloat
(5) -> gslIntegrationQng(g,0.0::DF,1.0::DF)
Compiling function g with type DoubleFloat -> DoubleFloat

(5) 0.7468241328124271
Type:
DoubleFloat
(6) ->

gsl.lisp

gsl.spad

Bill Page

unread,

Oct 25, 2015, 10:54:25 PM10/25/15

to fricas-devel

You are FAST! That's amazing. Works exactly as advertised.

Re qng:

"The function returns the final approximation, result, an estimate of
the absolute error, abserr and the number of function evaluations
used, neval."

https://www.gnu.org/software/gsl/manual/html_node/QNG-non_002dadaptive-Gauss_002dKronrod-integration.html

This is great.

Bill.

Alasdair McAndrew

unread,

Oct 25, 2015, 11:10:05 PM10/25/15

to fricas...@googlegroups.com

Sweet! Works beautifully (once I set my SBCL_HOME variable); one line will do it:

(2) -> gslIntegrationQng(x+->exp(-x^2),0.0,1.0)

And I suppose adding other routines is just a matter of more (include ...) lines in gsl.input and corresponding lines in gsl.spad?

With regard to arbitrary precision, that might be another issue. However, I was recommended that it may be better to "roll our own" in which case the method of choice is probably tanh-sinh integration, which has all of the advantages of Gauss Kronrod (nested rules for error analysis) and is far easier to implement. See http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/dhb-tanh-sinh.pdf for some details.

--

Bill Page

unread,

Oct 26, 2015, 12:14:36 AM10/26/15

to fricas-devel

Alasdair and Kurt,

Here is one way to deal with multiple return values:

(1) -> )co gsl
Compiling FriCAS source code from file /home/wspage/gsl.spad using
old system compiler.

------------------------------------------------------------------------
gsl is now explicitly exposed in frame frame1
gsl will be automatically loaded when needed from

/home/wspage/GSL.NRLIB/GSL

(1) -> )r gsl-test
DF ==> DoubleFloat

Type: Void

f:=(x:DF):DF+->exp(-x^2)

(2) theMap(*1;anonymousFunction;0;frame1;internal)
Type: (DoubleFloat -> DoubleFloat)

gslIntegrationQng(f,0.0::DF,1.0::DF)

(3) [result= 0.7468241328124271,abserr= 8.291413475940725e-15,neval= 21]
Type: Record(result: DoubleFloat,abserr: DoubleFloat,neval: Integer)
(4) ->

---

In the revised files attached 'integration-qng-list' is a wrapper
routine that just returns the multiple values as a list of values.
Since lists are homogeneous in FriCAS we then have to sort out the
actual types of the result for SPAD in gslIntegrationQng and in this
example return it as a record.

---

Yes, adding other routines should be easy. The hard part of passing
FriCAS functions was done by Kurt (he just made it look easy! :)

I agree that arbitrary precision is indeed another issue and in my
opinion, in general not so interesting.

More tomorrow when I get a chance. I am wondering now about how best
to make GSL a permanent addition to FriCAS. Perhaps Waldek can advise.

Bill.

gsl.spad

gsl-test.input

gsl.lisp

Bill Page

unread,

Oct 26, 2015, 12:34:23 PM10/26/15

to fricas-devel

I just dropped this "proof-of-concept" level code into a repository at:

https://github.com/billpage/gsla

Title: GNU Scientific Library for Axiom (and FriCAS and OpenAxiom)

I hope there are no objections or let me know.

Please check it out and test it. Feel free to fork it, make changes, etc.

For now I have included Kurt, Alasdair and myself as "authors" though
really so far only Kurt has written any significant code.

OK?

Bill.

Kurt Pagani

unread,

Oct 26, 2015, 3:19:57 PM10/26/15

to fricas...@googlegroups.com

Hello Bill,

that's great. There is a lot in Lisp which is a "black spot" to me,
especially this "multiple value" stuff. Since this also seems to be
solved I think this could be a starting point.

Thank you and
best regards
Kurt

Kurt Pagani

unread,

Oct 26, 2015, 3:32:23 PM10/26/15

to fricas...@googlegroups.com

That's very "ok" :)
Unfortunately I have little time to contribute to this topic, but you
know much more about these things anyway. I also had a glance over the
FFI of OpenAxiom which Alfredo Portes posted recently here. Looks quite
interesting but here too, I have to avoid frittering ;)

I'll keep an eye on the repo :)

Many thanks
Kurt

Alasdair McAndrew

unread,

Oct 26, 2015, 8:44:32 PM10/26/15

to fricas...@googlegroups.com

Now that Karl has provided a fantastic working template for an integration routine from QUADPACK/GSL/GSLL it should be possible to extend those files to incorporate some of the other routines. For the next few days I'll not have much time, but I'll aim to fiddle later in the week. And I might even learn a little Lisp in the process!

Bill Page

unread,

Oct 26, 2015, 10:37:17 PM10/26/15

to fricas-devel

Here is a little design issue concerning names.

Right now we have the example of 'gslIntegrationQng' and I have
started to fill in some miscellaneous functions such as

gslLookup: String -> Symbol
++ \spad{\gslLookup} finds GSLL function by GSL name and
++ displays some documentation.
gslIsNaN?: DF -> Boolean
++ \spad{\gslIsNaN?} returns true if x is not-a-number.
gslIsInf: DF -> SingleInteger
++ \spad{\gslIsInf} returns +1 if x is positive infinity,
++ -1 if x is negative infinity and 0 otherwise. On
++ some platforms 1 is returned for negative infinity.
gslFinite?: DF -> Boolean
++ \spad{\gslFinite?} returns true if x is a real number,
++ and false if it is infinite or not-a-number.
gslFcmp: (DF,DF,DF) -> SingleInteger
++ \spad{\gslFcmp} determines whether x and y are approximately
++ equal to a relative accuracy epsilon.
gslIntegrationQng : (DF -> DF,DF,DF) -> Record(result:DF,
abserr:DF, neval:Integer)
++ \spad{\gslIntegrationQng} applies the Gauss-Kronrod 10-point,
++ 21-point, 43-point and 87-point integration rules in succession

but before going further it occurs to me to wonder if we really should
continue to use the prefix 'gsl' in all these function names? In
Axiom/FriCAS it is common to overload names and rely on package names
and types for disambiguation. It seems unnecessarily redundant to
write:

gslIntegrationQng$gsl

Also, the names that appear at the GSLL (lisp) level omit the gsl
prefix for a similar reason - Common Lisp packaging.

So I propose that we drop the 'gsl' part of the name and lowercase the
first letter, as is the Axiom custom.

lookup: String -> Symbol
++ \spad{\lookup} finds GSLL function by GSL name and
++ displays some documentation.
isNaN?: DF -> Boolean
++ \spad{\isNaN?} returns true if x is not-a-number.
isInf: DF -> SingleInteger
++ \spad{\isInf} returns +1 if x is positive infinity,
++ -1 if x is negative infinity and 0 otherwise. On
++ some platforms 1 is returned for negative infinity.
finite?: DF -> Boolean
++ \spad{\finite?} returns true if x is a real number,
++ and false if it is infinite or not-a-number.
fcmp: (DF,DF,DF) -> SingleInteger
++ \spad{\fcmp} determines whether x and y are approximately
++ equal to a relative accuracy epsilon.
integrationQng : (DF -> DF,DF,DF) -> Record(result:DF, abserr:DF,
neval:Integer)
++ \spad{\integrationQng} applies the Gauss-Kronrod 10-point,
++ 21-point, 43-point and 87-point integration rules in succession

What do you think?

On 26 October 2015 at 20:44, Alasdair McAndrew <amc...@gmail.com> wrote:
>
> Now that Karl has provided a fantastic working template for an integration
> routine from QUADPACK/GSL/GSLL it should be possible to extend those
> files to incorporate some of the other routines. For the next few days I'll
> not have much time, but I'll aim to fiddle later in the week. And I might
> even learn a little Lisp in the process!
>

My plan is to chip away at adding new routines as time permits. If
there is something in particular you would like to tackle, just let me
know.

Besides GSL are there some other libraries that might be interesting
and amenable to the QuickLisp approach?

Alasdair McAndrew

unread,

Oct 26, 2015, 11:20:44 PM10/26/15

to fricas...@googlegroups.com

We can just go the Maxima approach, and write quadQag, quadQng etc (Maxima has quad_qag, quad_qng) - this retains the naming convention from QUADPACK. Or we could have a domain QUAD which included operations qag, qng etc... like this:

v := quad(x+->exp(-x^2),0,1,1.0e-10)

qag(v)

As to other quicklisp libraries, BLAS, LAPACK? What about CL-BUCHBERGER - a common lisp implementation of Buchberger's algorithm for Groebner base computation? The place to search is quickdocs.org.

Note that QUADPACK is far from being the last word of numeric quadrature routines: I have been experimenting with the integral of sin(1/x) between x=0 and x=1, and pretty much every routine quits in disgust with an error. That's why I would also like to implement a scheme for tanh-sinh integration.

Alasdair McAndrew

unread,

Oct 27, 2015, 8:25:19 AM10/27/15

to fricas...@googlegroups.com

Sorry to be an idiot - but where's a good place to go to learn about spad/boot/lisp insofar as they relate to PanAxiom development? I've been looking over Kurt's gsl.lisp and gsl.spad and there are quite a few things which confuse me. I won't list them, as it would be a depressingly long list for two very short files...

Ralf Hemmecke

unread,

Oct 27, 2015, 8:50:23 AM10/27/15

to fricas...@googlegroups.com

On 10/27/2015 01:25 PM, Alasdair McAndrew wrote:
> Sorry to be an idiot

You are certainly not an idiot. The problem is that there is just not
every documentation put in such a way that one can easily find it. And
some documentation isn't even existing. In any case, I encourage you to
ask here on fricas-devel.

> - but where's a good place to go to learn about
> spad/boot/lisp insofar as they relate to PanAxiom development?

This one will certainly be helpful.

http://axiom-wiki.newsynthesis.org/images/axiom--test--1/src/boot/Makefile.pdf

I don't have a good link for LISP related to FriCAS.

But SPAD you can learn from the AXIOM book.

http://fricas.github.io/book.pdf

But to understand the principles, it's probably enough to read my
tutorial about SpadProgramming.

http://axiom-wiki.newsynthesis.org/ProgrammingSPAD

Furthermore (although it is not exactly SPAD), take a look into the
Aldor User Guide. Yes, there are differences:
http://axiom-wiki.newsynthesis.org/LanguageDifferences
but they are not so important. And overall the Aldor User Guide
concentrates on explaining the language while the AXIOM book focuses
more on how to use PanAxiom, but explains programming more as a side
remark than in full detail.

Ralf

Waldek Hebisch

unread,

Oct 27, 2015, 9:51:35 AM10/27/15

to fricas...@googlegroups.com

Alasdair McAndrew wrote:
>
> spad/boot/lisp insofar as they relate to PanAxiom development? I've been
> looking over Kurt's gsl.lisp and gsl.spad and there are quite a few things
> which confuse me. I won't list them, as it would be a depressingly long
> list for two very short files...

Why not? There is about 40 code lines in the two file, I hope
50-60 questions will cover it.

Concerning Lisp, the HyperSpec should anwer most your
questions:

http://www.lispworks.com/documentation/HyperSpec/Front/

--
Waldek Hebisch

Bill Page

unread,

Oct 27, 2015, 9:54:53 AM10/27/15

to fricas-devel

Alasdair,

There is stuff to read such as Ralf suggests but in my opinion just
hitting your head against something relevant hard enough for long
enough is much more effective :) Kurt's example code is short enough
not to cause much harm, so why not start just by asking one to two
questions rather than making a list? As I am sure you tell your
students, the first step is not letting it make you feel intimidated.
Just ask.

Concerning my question about naming conventions: My preference is to
avoid all redundant prefixes so Maxima isn't such a great model. I am
not sure I understand the intention of your examples

>> v := quad(x+->exp(-x^2),0,1,1.0e-10)
>> qag(v)

but yes, that is more or less what I would like. The term "domain" in
Axiom has a particular meaning: a data/mathematical type or class.
For the purposes of extending Axiom to include more operations on
existing domains, e.g. DoubleFloat, a "package" is more appropriate.
Because GSL is a rather large library, it might be appropriate to
define more than one package but right now we have only one named
'gsl'.

Actually, using a package/domain name like 'gsl' consisting of all
lower case letters is contrary to the general conventions in the
Axiom/FriCAS library. Unlike variable and functions names, these names
are supposed to start with a capital letter. Each package/domain also
has a short (< 6 characters) abbreviation, by convention written in
all upper case. So we should probably at least change the package
name from 'gsl' to 'Gsl'. In keeping with what I said above, maybe we
would want package names such as 'GSLintegrate' etc.

As you see in examples, in many contexts we can avoid explicitly using
the package/domain name since it is automatically deduced by the
interpreter.

Concerning your suggested packages: Note that there is already
extensive support in Axiom/FriCAS for Groebner basis compuations. Did
you have something specific in mind? BLAS and LAPACK might be
interesting. As I understand it there are already interfaces to these
packages in GSL.

One thing I tried last night (unsuccessfully) was to call
'lu-decomposition' the GSLL equivalent of 'gsl_linalg_LU_decomp'. I
have a problem converting the matrix passed by FriCAS to the required
matrix format. This is the kind of thing we need to resolve for using
more general routines.

BTW, There is or once was a policy promoting the use of plain text
email format on mailing lists. I notice yours arrive with a more
modern look. There might be some people around who prefer the older
format.

Bill.

> --
>
> --

> You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.

Alasdair McAndrew

unread,

Oct 27, 2015, 7:55:18 PM10/27/15

to fricas...@googlegroups.com

Orright... two questions about https://github.com/billpage/gsla/blob/master/gsl.spad, and in particular the final three lines:

gslIntegrationQng(f,a,b) ==
r:List DF:=INTEGRATION_-QNG_-LIST(mkLispFunction1(f@(DF->DF))$Lisp,a,b)$Lisp
[r(1),r(2),r(3) pretend Integer]

First the use of "pretend". As I understand it, "pretend" is a sort of syntactic fudge factor to allow Axiom to compile code in which the types may not be precisely fixed, or in which there may be a type mismatch. But why is it needed particularly here - in the definition of gslIntegrationQng? Surely an object is either an integer or it isn't, I would have thought, and we can define the final value "neval" to be an integer.

Second, the underscores in "INTEGRATION_-QNG_-LIST". Are these just to escape the hyphens, lest SPAD confuse them with minus signs? After all, the relevant function is defined in gsl.lisp as

defun integration-qng-list (f a b)
(multiple-value-list (integration-qng f a b)))

which incidentally would appear to indicate that one of the languages here is case-insensitive. Or is there another reason?

With regard to plain text format, I fully agree, and I'll do that if people want. On the other hand, it's nice to have some way of distinguishing code samples from text - the use of markdown is good for this, but I don't know if google groups supports it. (Well, it sort of does, but through a browser extension called "Markdown Here" which also supports LaTeX.)

I do find the mixture of SPAD and Lisp with their individual conventions and behaviours very confusing.

<Australian Accent>I tells ya, there's nuthin' like Axiom for makin' a bloke feel stupid!</Australian Accent>

-Alasdair

Alasdair McAndrew

unread,

Oct 27, 2015, 8:41:05 PM10/27/15

to fricas...@googlegroups.com

A quickie: ")read /home/amca/quicklisp/setup.lisp" works in FriCAS, but not in efricas. Why not? Also, the history mechanism in FriCAS in a console seems a bit odd: sometimes (well, often), a previous command, obtained using the arrow keys, will have extra characters, or chanrecters missing (I'm using Konsole: the KDE console, with my term=xterm.)

Bill Page

unread,

Oct 27, 2015, 9:10:20 PM10/27/15

to fricas-devel

On 27 October 2015 at 19:55, Alasdair McAndrew <amc...@gmail.com> wrote:
>
> Orright... two questions about https://github.com/billpage/gsla/blob/master/gsl.spad,
> and in particular the final three lines:
>
> gslIntegrationQng(f,a,b) ==
> r:List DF:=INTEGRATION_-QNG_-LIST(mkLispFunction1(f@(DF->DF))$Lisp,a,b)$Lisp
> [r(1),r(2),r(3) pretend Integer]
>
> First the use of "pretend". As I understand it, "pretend" is a sort of syntactic fudge
> factor to allow Axiom to compile code in which the types may not be precisely fixed,
> or in which there may be a type mismatch.

Well, sort of ... The first thing to remember is that underneath all
of Axiom/FriCAS there is Lisp, or maybe better said: Lisp provides a
kind of high level virtual machine or runtime environment, so
ultimately at the "bottom" every datum is represented by some Lisp
object. But unlike Axiom itself , Lisp objects are dynamically typed
- each value comes with it's own type - while in Axiom a value has a
type only in so far as it belongs to some domain. There is more that
can be said about this but enough for now.

The purpose of 'pretend' is to specify an Axiom domain for a given
value. This value may be considered by Axiom to belong to some other
domain but 'pretend' overrides this - kind of like a "cast" in other
languages or it might originate (as in our case) from the Lisp level
and so belong to no particular Axiom domain. (You should ask later on
about the use of the domain SExpression). In any case the use of
'pretend' is not type-safe. It is possible (even easy) to produce
unpredictable results and even Lisp errors if a value is treated as
belonging to a domain with an incompatible representation.

> But why is it needed particularly here - in the definition of gslIntegrationQng?
> Surely an object is either an integer or it isn't, I would have thought,

At the Lisp level an object has a type, so yes Lisp knows that the
third element of the list is an integer (or more specifically, it's
type is UNSIGNED-BYTE 64). but Axiom known nothing about that so we
are just telling it to treat it as if it belongs to the Integer
domain. Whether or not this is a good idea is highly dependent on the
nitty gritty details of the representation of Integer in Axom and the
implementation of UNSIGNED-BYTE 64 in Lisp (which may vary from one
implementation of Lisp, say SBCL, to another, say CLisp or ECL). It
works for us in SBCL but I could not guarantee that it will in other
Lisps. This is something to study more deeply and test in various
situations.

> and we can define the final value "neval" to be an integer.

But Axiom can't so that is exactly why the programmer needs to use 'pretend'.

>
> Second, the underscores in "INTEGRATION_-QNG_-LIST". Are these just
> to escape the hyphens, lest SPAD confuse them with minus
> signs?

Ah, an easy question with an easy answer: Yes.

> After all, the relevant function is defined in gsl.lisp as
>
> defun integration-qng-list (f a b)
> (multiple-value-list (integration-qng f a b)))
>
> which incidentally would appear to indicate that one of the languages here
> is case-insensitive. Or is there another reason?

Lisp is case-insensitive by default. Lisp uses a special notation if
symbols need to be in both upper and lower case. I think Kurt
mentioned this earlier. In Lisp

|abc|

with vertical bars, represents the symbol 'abc', while

abc

without vertical bars, represents the symbol 'ABC'.

>
> With regard to plain text format, I fully agree, and I'll do that if people
> want. On the other hand, it's nice to have some way of distinguishing
> code samples from text - the use of markdown is good for this, but I
> don't know if google groups supports it. (Well, it sort of does, but
> through a browser extension called "Markdown Here" which also
> supports LaTeX.)

I am flexible.

>
> I do find the mixture of SPAD and Lisp with their individual conventions
> and behaviours very confusing.
>

Granted. I did mention the BOOT language earlier, right? :) But in
Axiom terms this situation isn't very confusing (yet)!

> <Australian Accent>I tells ya, there's nuthin' like Axiom for makin' a
> bloke feel stupid!</Australian Accent>
>

Too bad there isn't a good way to include THAT in an email.

Cheers,
Bill.

Bill Page

unread,

Oct 27, 2015, 9:20:57 PM10/27/15

to fricas-devel

On 27 October 2015 at 20:41, Alasdair McAndrew <amc...@gmail.com> wrote:
>
> A quickie: ")read /home/amca/quicklisp/setup.lisp" works in FriCAS, but not
> in efricas. Why not?

Hmmm, not sure. I usually don't use efricas (emacs drives me crazy).
Maybe it is better to use something like

)lisp (load "quicklisp/setup.lisp")

> Also, the history mechanism in FriCAS in a console seems a bit odd: sometimes
> (well, often), a previous command, obtained using the arrow keys, will have extra
> characters, or chanrecters missing (I'm using Konsole: the KDE console, with my
> term=xterm.)
>

I suffer with the same sort of problem: backspace and inserting
corrections usually does not work. Axiom (and FriCAS by default) uses
an ancient program called 'clef' instead of the now much more common
'readline'. It is possible in principle to turn off clef when starting
FriCAS, e.g.

$ fricas -noclef

but then you have no history mechanism unless SBCL itself was compiled
with the option to include readline. I don't usually bother to try to
make this work. Maybe Waldek can say more about when and how this
works.

Cheers,
Bill.

Waldek Hebisch

unread,

Oct 27, 2015, 9:34:17 PM10/27/15

to fricas...@googlegroups.com

Alasdair McAndrew wrote:
>
> Orright... two questions about
> https://github.com/billpage/gsla/blob/master/gsl.spad, and in particular
> the final three lines:
>
> gslIntegrationQng(f,a,b) ==
> r:List DF:=INTEGRATION_-QNG_-LIST(mkLispFunction1(f@(DF->DF))$Lisp,a,b)$Lisp
> [r(1),r(2),r(3) pretend Integer]
>
> First the use of "pretend". As I understand it, "pretend" is a sort of
> syntactic fudge factor to allow Axiom to compile code in which the types
> may not be precisely fixed, or in which there may be a type mismatch. But
> why is it needed particularly here - in the definition of
> gslIntegrationQng? Surely an object is either an integer or it isn't, I
> would have thought, and we can define the final value "neval" to be an
> integer.

Expanding on Bill's answer: this code is an abuse (wrong). The
core thing is that FriCAS list are homogeneous: all elements have
the same type. In Lisp lists are allowed to have elements of
different Lisp types. The integration routine returns two
machine float and one machine integer. On Lisp side Bill packs
this into a list. But on Spad side this List in an illegal
object, so Bill lies and tells Spad compiler that the list
contains DoubleFloat-s. Then he has to use 'pretend' to
tell Spad compiler that the third element is in fact an
integer. Currently this works, but will break down if
Spad compiler becomes only slightly smarter: claiming that
something is both DoubleFloat and Integer is a contradiction
and may lead to wrong code or compile time error.

> Second, the underscores in "INTEGRATION_-QNG_-LIST". Are these just to
> escape the hyphens, lest SPAD confuse them with minus signs? After all, the
> relevant function is defined in gsl.lisp as
>
> defun integration-qng-list (f a b)
> (multiple-value-list (integration-qng f a b)))
>
> which incidentally would appear to indicate that one of the languages here
> is case-insensitive. Or is there another reason?

Underscores are just to escape hyphens. Note that hyphens
in "INTEGRATION-QNG-LIST" are a self-inflicted hardship.
On lisp side one can write:

(defun INTEGRATION_QNG_LIST (f a b) ....)

and on FriCAS side use the same "INTEGRATION_QNG_LIST". This is
because FriCAS treat underscores which are not needed as escapes
as literal underscores.

Concerning case-sensitivity, Lisp have schizophrenia here:
Lisp symbols are internally case sensitive, but by default
on input Lisp will uppercase them. My personal style is
to write symbols which are used from Spad in uppercase,
so that Spad version and Lisp version looks the same.

>
> With regard to plain text format, I fully agree, and I'll do that if people
> want. On the other hand, it's nice to have some way of distinguishing code
> samples from text - the use of markdown is good for this, but I don't know
> if google groups supports it. (Well, it sort of does, but through a
> browser extension called "Markdown Here" which also supports LaTeX.)

FYI, for me non-text part is junk that I delete without looking
at. Using plain text save me work on deleting it.

--
Waldek Hebisch

Bill Page

unread,

Oct 27, 2015, 10:11:31 PM10/27/15

to fricas-devel

On 27 October 2015 at 21:34, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
>
> Expanding on Bill's answer: this code is an abuse (wrong). The
> core thing is that FriCAS list are homogeneous: all elements have
> the same type. In Lisp lists are allowed to have elements of
> different Lisp types. The integration routine returns two
> machine float and one machine integer. On Lisp side Bill packs
> this into a list. But on Spad side this List in an illegal
> object, so Bill lies and tells Spad compiler that the list
> contains DoubleFloat-s. Then he has to use 'pretend' to
> tell Spad compiler that the third element is in fact an
> integer. Currently this works, but will break down if
> Spad compiler becomes only slightly smarter: claiming that
> something is both DoubleFloat and Integer is a contradiction
> and may lead to wrong code or compile time error.
>

Yes, guilty. :) But this was only the first solution I could think of
and so far this is just proof-of-concept code. There are other ways
in Lisp of dealing with functions that return multiple values but this
is only possible in FriCAS at a higher level, e.g. by use of Record.
Perhaps there is some easy way of using routines already available in
BOOT to construct a Record value directly? It would take me some time
to find it. Do you know Waldek? Or can you suggest a different
possible permanent solution?

Alasdair McAndrew

unread,

Oct 28, 2015, 1:40:35 AM10/28/15

to fricas...@googlegroups.com

Does SPAD support optional arguments to functions? The default absolute error (as defined in gsll: calculus/numerical-integration.lisp) is only 1.0d-5. We might wish to decrease this error. In which case the definition of gslIntegrationQng in gsl.spad needs to include some optional arguments which can then be passed to the function in gsl.lisp.

Ralf Hemmecke

unread,

Oct 28, 2015, 2:31:26 AM10/28/15

to fricas...@googlegroups.com

On 10/28/2015 02:34 AM, Waldek Hebisch wrote:
>> >
>> > With regard to plain text format, I fully agree, and I'll do that if people
>> > want. On the other hand, it's nice to have some way of distinguishing code
>> > samples from text - the use of markdown is good for this, but I don't know
>> > if google groups supports it. (Well, it sort of does, but through a
>> > browser extension called "Markdown Here" which also supports LaTeX.)
> FYI, for me non-text part is junk that I delete without looking
> at. Using plain text save me work on deleting it.

Similar for me. I usually don't look into
https://groups.google.com/forum/#!forum/fricas-devel. The reason is
simple, it uses a non-monospace font and thus formatted code looks ugly
and unreadable.

Often HTML mail comes with am additional plain text part.
Fortunately my mail program automatically hides the unreadable html part
from me, so it's only wasted memory and bandwidth.
If someone sends HTML-only, my mailer would show me an empty mail
(which, of course, I delete as junk).
HTML adds nothing to emails that I highly value.

Ralf

Ralf Hemmecke

unread,

Oct 28, 2015, 2:43:23 AM10/28/15

to fricas...@googlegroups.com

> Does SPAD support optional arguments to functions?

Aldor does, but to my knowledge SPAD does not. But you can overload
function names, see below.

Does that help?

Ralf

-- foo.spad
)abbrev package FOO Foo
Foo(T: Ring): with
foo: T -> T
foo: (T, T) -> T
== add
foo(x: T, y: T): T == 10*x+y
foo(x: T): T == foo(x, 1)
------------------------

(1) -> )co foo.spad

(1) -> foo(3)$Foo(Integer)

(1) 31
Type: Integer
(2) -> foo(1,7)$Foo(Integer)

(2) 17
Type: Integer
(3) -> foo(4.0)$Foo(Float)

(3) 41.0
Type: Float
(4) -> foo(1.0, 1.9)$Foo(Float)

(4) 11.9
Type: Float

Waldek Hebisch

unread,

Oct 28, 2015, 8:17:29 AM10/28/15

to fricas...@googlegroups.com

Bill Page wrote:
>
> Perhaps there is some easy way of using routines already available in
> BOOT to construct a Record value directly? It would take me some time
> to find it. Do you know Waldek? Or can you suggest a different
> possible permanent solution?

It seems there are no ready routine. But look at 'optMkRecord'
in 'g-opt.boot': it will generate Lisp code building a Record.

--
Waldek Hebisch

Alasdair McAndrew

unread,

Oct 28, 2015, 10:39:34 PM10/28/15

to fricas...@googlegroups.com

I tried to define

U==>Union(DoubleFloat, Integer)

and then have the return list r to be of type U:

r:List U:=...

Then I thought r could happily include both DF's and integers. But it doesn't work in the way I hoped it would. Drat.

Waldek Hebisch

unread,

Oct 28, 2015, 10:58:58 PM10/28/15

to fricas...@googlegroups.com

Alasdair McAndrew wrote:
>
> I tried to define
>
> U==>Union(DoubleFloat, Integer)
>
> and then have the return list r to be of type U:
>
> r:List U:=...
>
> Then I thought r could happily include both DF's and integers. But it
> doesn't work in the way I hoped it would. Drat.

Well, at low level elements of Union consist of tag and value.
Tag if a small integer identifying which variant is at force.
More precisely, Spad union is interpreted as Lisp pair,
consisting of tag and value:

(9) -> PRETTYPRINT(2::U)$Lisp
(0 . 2.0)

(9) ()
Type: SExpression

The dot above indicates that this is not a list but
a general pair (pairs in lists have rest of list as second
element). For some reason interpreter decided to treat 2
as a double, so we got first variant (with tag 0).

Having list of such pairs is quite different than list
which mixes integers and doubles:

(10) -> PRETTYPRINT([2::U])$Lisp
((0 . 2.0))

(10) ()
Type: SExpression
(11) -> PRETTYPRINT([2])$Lisp
(2)

(11) ()
Type: SExpression
(12) -> PRETTYPRINT([2::DoubleFloat])$Lisp
(2.0)

(12) ()
Type: SExpression

--
Waldek Hebisch

Alasdair McAndrew

unread,

Oct 28, 2015, 11:18:25 PM10/28/15

to fricas...@googlegroups.com

Not easy is it? We want a return value (of a function defined in a spad file) to consist of two DoubleFloat values and an integer. As you know, Bill Page had to use a list of DoubleFloats one of whose values was of type "pretend Integer", which is a kludge - although it works. We are told that spad doesn't support records - but is there some natural way of returning values of two different types?

Ralf Hemmecke

unread,

Oct 29, 2015, 5:17:45 AM10/29/15

to fricas...@googlegroups.com

>> U==>Union(DoubleFloat, Integer)

> (9) -> PRETTYPRINT(2::U)$Lisp
> (0 . 2.0)
>
> (9) ()
> Type: SExpression

When I look at this

(1) -> DF==>DoubleFloat;Z==>Integer
(2) -> R==>Record(result:DF, abserr:DF, neval:Integer)
(6) -> r: R := [2.0, 3.0, 17]

(6) [result= 2.0,abserr= 3.0,neval= 17]

Type: Record(result: DoubleFloat,abserr: DoubleFloat,neval: Integer)

(7) -> PRETTYPRINT(r)$Lisp
#(2.0 3.0 17)

it looks as if a record is internally stored as a lisp vector, no?

Concerning (from Alasdair's gsl.lisp)

; Return multiple values as a list
(defun integration-qng-list (f a b epsabs epsrel)
(multiple-value-list (integration-qng f a b epsabs epsrel)))

Maybe it would be easiest for Alasdair to turn the lisp value (that is
returned from the GSL function integration-png) into such an array.
I.e. instead of multiple-value-list he should call a function to produce
the array/record. (Sorry, no idea what that function is.)

In the .spad file instead of

r:List
DF:=INTEGRATION_-QNG_-LIST(mkLispFunction1(f@(DF->DF))$Lisp,a,b,epsabs,epsrel)$Lisp
[r(1),r(2),r(3) pretend Integer]

there would then simply by

R ==> Record(result:DF, abserr:DF, neval:Integer)
r: R :=
NTEGRATION_-QNG_-LIST(mkLispFunction1(f@(DF->DF))$Lisp,a,b,epsabs,epsrel)$Lisp

Ralf

Bill Page

unread,

Oct 29, 2015, 12:20:41 PM10/29/15

to fricas-devel

On 29 October 2015 at 05:17, Ralf Hemmecke <ra...@hemmecke.org> wrote:

> (7) -> PRETTYPRINT(r)$Lisp
> #(2.0 3.0 17)
> it looks as if a record is internally stored as a lisp vector, no?

No not quite. That is a literal vector in Lisp. But you are right
that this is the actual internal representation of something of type
Record. E.g.

(1) -> VECTOR(1.0::DF,2.0::DF,3)$Lisp pretend Record(a:DF,b:DF,c:Integer)

(1) [a= 1.0,b= 2.0,c= 3]
Type: Record(a: DoubleFloat,b: DoubleFloat,c: Integer)

>
> Concerning (from Alasdair's gsl.lisp)
>
> ; Return multiple values as a list
> (defun integration-qng-list (f a b epsabs epsrel)
> (multiple-value-list (integration-qng f a b epsabs epsrel)))
>
> Maybe it would be easiest for Alasdair to turn the lisp value (that is
> returned from the GSL function integration-png) into such an array.
> I.e. instead of multiple-value-list he should call a function to produce
> the array/record. (Sorry, no idea what that function is.)
>

Of course. Waldek did suggest where to look in the source code for
this, but it is certainly not the easiest option.

> In the .spad file instead of
>
> r:List
> DF:=INTEGRATION_-QNG_-LIST(mkLispFunction1(f@(DF->DF))$Lisp,a,b,epsabs,epsrel)$Lisp
> [r(1),r(2),r(3) pretend Integer]
>
> there would then simply by
>
> R ==> Record(result:DF, abserr:DF, neval:Integer)
> r: R :=
> NTEGRATION_-QNG_-LIST(mkLispFunction1(f@(DF->DF))$Lisp,a,b,epsabs,epsrel)$Lisp
>

Instead of speculating perhaps you should try it first? :)

>> Error detected within library code:
Bug: ridiculous record representation

But you had the right idea. This works:

; gsl.lisp
...
; Return multiple values as a vector (can be interpreted as Axiom Record)
(defun integration-qng-vector (f a b)
(apply #'vector (multiple-value-list (integration-qng f a b))))

-- gsl.spad:
...
gslIntegrationQng: (DF -> DF,DF,DF) -> Record(result:DF, abserr:DF,
neval:Integer)
...
gslIntegrationQng(f,a,b) ==
INTEGRATION_-QNG_-VECTOR(mkLispFunction1(f@(DF->DF))$Lisp,a,b)$Lisp

Bill Page

Bill Page

unread,

Oct 29, 2015, 12:25:28 PM10/29/15

to fricas-devel

Sorry. Previous email should say:

On 29 October 2015 at 12:20, Bill Page <bill...@newsynthesis.org> wrote:
> On 29 October 2015 at 05:17, Ralf Hemmecke <ra...@hemmecke.org> wrote:
>
>> (7) -> PRETTYPRINT(r)$Lisp
>> #(2.0 3.0 17)
>> it looks as if a record is internally stored as a lisp vector, no?
>

> *You are right* that this is the actual internal representation of something of
> type Record.
> ...

Waldek Hebisch

unread,

Oct 29, 2015, 12:59:07 PM10/29/15

to fricas...@googlegroups.com

Bill Page wrote:
>
> On 29 October 2015 at 05:17, Ralf Hemmecke <ra...@hemmecke.org> wrote:
>
> > (7) -> PRETTYPRINT(r)$Lisp
> > #(2.0 3.0 17)
> > it looks as if a record is internally stored as a lisp vector, no?
>
> No not quite. That is a literal vector in Lisp. But you are right
> that this is the actual internal representation of something of type
> Record. E.g.
>
> (1) -> VECTOR(1.0::DF,2.0::DF,3)$Lisp pretend Record(a:DF,b:DF,c:Integer)
>
> (1) [a= 1.0,b= 2.0,c= 3]
> Type: Record(a: DoubleFloat,b: DoubleFloat,c: Integer)

Note: this is for records of length > 2. Records of lennght 2 are
represented as Lisp pairs, the same as unions. In fact, without
knowing Spad type it is impossible to distinguish record from
union:

(11) -> U ==> Union(DoubleFloat, Integer)
Type: Void
(12) -> u := 2::U

(12) 2.0
Type: Union(DoubleFloat,...)
(13) -> PRETTYPRINT(u)$Lisp
(0 . 2.0)

(13) ()
Type: SExpression
(14) -> R ==> Record(ii : Integer, dd : DoubleFloat)
Type: Void
(15) -> r := [0, 2.0]$R

(15) [ii= 0,dd= 2.0]
Type: Record(ii: Integer,dd: DoubleFloat)
(16) -> PRETTYPRINT(r)$Lisp
(0 . 2.0)

(16) ()
Type: SExpression
(17) -> EQUAL(u, r)$Lisp

(17) T
Type: SExpression

> But you had the right idea. This works:
>
> ; gsl.lisp
> ...
> ; Return multiple values as a vector (can be interpreted as Axiom Record)
> (defun integration-qng-vector (f a b)
> (apply #'vector (multiple-value-list (integration-qng f a b))))
>

Yes. Let me just remark that the code above is fine for
experimantation. In long term we should not sprinkle
assumptions about representation of vectors all around.
So probably we should create something to create
records. Note that currently representation of record
only depends on its lenght (number of fields), but in
principle could depend on type. So one possibility
is to add real 'construct' function to Record and
get it if needed via lookup. Another it to
provide function that given type and values would
create record of given type.

Maybe extra remark why I am making fuss here: in the past
FriCAS used quite different representation of two dimensional
arrays. We could change representation because it was only
used via small number of functions and macros.

--
Waldek Hebisch

Ralf Hemmecke

unread,

Oct 29, 2015, 5:57:41 PM10/29/15

to fricas...@googlegroups.com

>> ; Return multiple values as a vector (can be interpreted as Axiom
>> Record) (defun integration-qng-vector (f a b) (apply #'vector
>> (multiple-value-list (integration-qng f a b))))

> Yes. Let me just remark that the code above is fine for
> experimantation. In long term we should not sprinkle assumptions
> about representation of vectors all around.

Right. I wouldn't want foo$Lisp in every second contribution package.

> So probably we should create something to create records.

Most welcome Is a specification and/or an API how one can get/pass data
from/to SPAD. Is the FFI in openaxiom such a thing? Can we copy it?

Ralf

Alasdair McAndrew

unread,

Oct 29, 2015, 6:57:00 PM10/29/15

to fricas...@googlegroups.com

A quick note: gsl.lisp and gsl.spad were originally written by Bill Page and made available on his site. So the fact that they work is entirely due to him. I have attempted to make some changes, but rather like barnacles on a ship's hull they don't do much except make the code more confusing. So you should refer to "Bill's gsl.spad" if it works, and "Alasdair's gsl.spad" if it doesn't!

Bill Page

unread,

Oct 29, 2015, 9:02:04 PM10/29/15

to fricas-devel

On 29 October 2015 at 18:56, Alasdair McAndrew <amc...@gmail.com> wrote:
>
> A quick note: gsl.lisp and gsl.spad were originally written by Bill Page
> and made available on his site. So the fact that they work is entirely due
> to him. I have attempted to make some changes, but rather like barnacles
> on a ship's hull they don't do much except make the code more confusing.
> So you should refer to "Bill's gsl.spad" if it works, and "Alasdair's gsl.spad"
> if it doesn't!
>

I always liked the sound of the git "blame" command ... :) In fact as
I tried to say clearly much earlier in this thread, the original
versions of 'gsl.lisp' and 'gsl.spad' were written by Kurt. Alasdair
and I are just making more or less successful changes to Kurt's
original model code.

Bill.

Alasdair McAndrew

unread,

Oct 29, 2015, 9:31:51 PM10/29/15

to fricas...@googlegroups.com

Yes, in fact the chain of programming goes Kurt, Bill, me.

On a side issue, I note that gsll works by providing a wrapper for gsl using the Foreign Function Interface (which is why we need the ffi libraries). This means that the SPAD function in gsl.spad calls a lisp function in gsl.lisp which calls the gsll function gsll:integration-qng which in turns calls the gsl function defined in /usr/local/src/gsl-1.16/integration/qag.c What this means is that if we wanted we could eliminate gsll completely and just use FFI to access gsl directly. (Of course I don't know how to do this!) But that's modern computing for you - there are chains of libraries and programs all depending on each other in different ways

Time for lunch.

Bill Page

unread,

Oct 29, 2015, 10:07:04 PM10/29/15

to fricas-devel

On 29 October 2015 at 21:31, Alasdair McAndrew <amc...@gmail.com> wrote:
>
> On a side issue, I note that gsll works by providing a wrapper for gsl using
> the Foreign Function Interface (which is why we need the ffi libraries).
> This means that the SPAD function in gsl.spad calls a lisp function in
> gsl.lisp which calls the gsll function gsll:integration-qng which in turns
> calls the gsl function defined in /usr/local/src/gsl-1.16/integration/qag.c

I see nothing wrong with this strategy since just calling C from Lisp
is the hardest part. GSLL solves many (but not all) of the problems
with type conversions.

> What this means is that if we wanted we could eliminate gsll completely
> and just use FFI to access gsl directly. (Of course I don't know how to do
> this!) But that's modern computing for you - there are chains of libraries
> and programs all depending on each other in different ways
>

Seems OK to me.

Waldek Hebisch

unread,

Oct 29, 2015, 10:28:29 PM10/29/15

to fricas...@googlegroups.com

Bill Page wrote:
>
> On 29 October 2015 at 21:31, Alasdair McAndrew <amc...@gmail.com> wrote:
> >
> > On a side issue, I note that gsll works by providing a wrapper for gsl using
> > the Foreign Function Interface (which is why we need the ffi libraries).
> > This means that the SPAD function in gsl.spad calls a lisp function in
> > gsl.lisp which calls the gsll function gsll:integration-qng which in turns
> > calls the gsl function defined in /usr/local/src/gsl-1.16/integration/qag.c
>
> I see nothing wrong with this strategy since just calling C from Lisp
> is the hardest part. GSLL solves many (but not all) of the problems
> with type conversions.

Not exactly. There are types in C and they map reasonably to
Spad types. Lisp encourages things which are awkward in
Spad. In both cases one needs to declare types at Spad
level and some cases needs wrappers. One needs to ensure
that libraries can be found, which for Lisp means that
we need to find bot Lisp library and C library.
--
Waldek Hebisch

Alasdair McAndrew

unread,

Oct 29, 2015, 11:41:35 PM10/29/15

to fricas...@googlegroups.com

Actually I don't know that gsll does fully solve the type problem. Gsll calls an FFI "groveller" method which does what it says: it "grovels" about and matches up types. So in a sense FFI does a great deal of the hard work and gsll is mainly a front end (although from my poking about the gsll directories it does do a lot of setting things up).

Kurt Pagani

unread,

Oct 29, 2015, 11:56:37 PM10/29/15

to fricas...@googlegroups.com

LOL, Alisdair was right: cui honorem, honorem
My humble contribution is common lore ... the back-breaking work lies
before you :)

Alasdair McAndrew

unread,

Oct 30, 2015, 12:10:16 AM10/30/15

to fricas...@googlegroups.com

Good to see you know your Vulgate. I like a little historical literacy in a news group.

Bill Page

unread,

Oct 30, 2015, 8:21:38 AM10/30/15

to fricas-devel

On 29 October 2015 at 22:28, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> Bill Page wrote:
>>
>> On 29 October 2015 at 21:31, Alasdair McAndrew <amc...@gmail.com> wrote:

>> > ...

>> > calls the gsl function defined in /usr/local/src/gsl-1.16/integration/qag.c
>>
>> I see nothing wrong with this strategy since just calling C from Lisp
>> is the hardest part. GSLL solves many (but not all) of the problems
>> with type conversions.
>
> Not exactly. There are types in C and they map reasonably to
> Spad types. Lisp encourages things which are awkward in
> Spad.

Surely less awkward than C.

> In both cases one needs to declare types at Spad
> level and some cases needs wrappers. One needs to ensure
> that libraries can be found, which for Lisp means that
> we need to find bot Lisp library and C library.

Well, Kurt showed how to call 'qng.c' will just 3 lines of code. I
would be happy to see a similar demonstration that involved Spad and
C.

Bill.

Alasdair McAndrew

unread,

Oct 30, 2015, 8:53:26 AM10/30/15

to fricas...@googlegroups.com

II am hoping for advice from folk higher up the food chain than me; I have an error message which is not helpful enough.

I am enlarging Kurt and Bill's gsl.* function to include the gsl/gsll function integration-qag. This can be invoked, for example as

)lisp (gsll:integration-qag (lambda (x) (exp (- (* x x)))) 0.0 1.0 :gauss21)

where the final :gauss21 is one of a number of parameters which define the integration rule to be used.

In my gsl.lisp file I have

(defvar intrules (list ':gauss15 ':gauss21 ':gauss31 ':gauss41 ':gauss51 ':gauss61))

(defun integration-qag-list (f a b n)
(multiple-value-list (integration-qag f a b (eval (nth n intrules)))))

and in my corresponding gsl.spad:

    gslIntegrationQag(f,a,b,n) ==
      r:List DF:=INTEGRATION_-QAG_-LIST(mkLispFunction1(f@(DF->DF))$Lisp,a,b,n)$Lisp
      [r(1),r(2),r(3) pretend Integer]

This all compiles OK. But in FriCAS:

(1) -> gslIntegrationQag(x+->exp(-x^2),0.0,1.0,2)

>> Error detected within library code:

index out of range

so clearly I've not set something up correctly - I suspect it's my cack-handed defvar list of parameters and my use of (nth n intrules). But I don't know why this isn't working, and what I should be doing instead.

Sigh. Maybe I need to drink more tea.

Waldek Hebisch

unread,

Oct 31, 2015, 12:58:54 AM10/31/15

to fricas...@googlegroups.com

Alasdair McAndrew wrote:
>
> I am enlarging Kurt and Bill's gsl.* function to include the gsl/gsll
> function integration-qag. This can be invoked, for example as
>
> )lisp (gsll:integration-qag (lambda (x) (exp (- (* x x)))) 0.0 1.0 :gauss21)
>
> where the final :gauss21 is one of a number of parameters which define the
> integration rule to be used.
>
> In my gsl.lisp file I have
>
> (defvar intrules (list ':gauss15 ':gauss21 ':gauss31 ':gauss41 ':gauss51
> ':gauss61))
>
> (defun integration-qag-list (f a b n)
> (multiple-value-list (integration-qag f a b (eval (nth n intrules)))))
>
> and in my corresponding gsl.spad:
>
> gslIntegrationQag(f,a,b,n) ==
> r:List DF:=INTEGRATION_-QAG_-LIST(mkLispFunction1(f@
> (DF->DF))$Lisp,a,b,n)$Lisp
> [r(1),r(2),r(3) pretend Integer]
>
> This all compiles OK. But in FriCAS:
>
> (1) -> gslIntegrationQag(x+->exp(-x^2),0.0,1.0,2)
>
> >> Error detected within library code:
> index out of range
>
> so clearly I've not set something up correctly - I suspect it's my
> cack-handed defvar list of parameters and my use of (nth n intrules). But
> I don't know why this isn't working, and what I should be doing instead.

Note that number 2 corresponds to :GAUSS31:

)lisp (defvar intrules (list ':gauss15 ':gauss21 ':gauss31 ':gauss41 ':ga)lisp (defvar intrules (list ':gauss15 ':gauss21 ':gauss31 ':gauss41 ':gauss51 ':gauss61))

Value = INTRULES
(1) -> )lisp (nth 2 intrules)

Value = :GAUSS31

AFAICS 'eval' in 'integration-qag-list' is not needed. But I do
not think this is reason for error. I would rather supect that
'integration-qag' returns less values than you expect. Anyway
in such situation I would do

)set break break

run the offending command

and then at Lisp debugger prompt (to avoid potentially too large
debugging printouts):

(setf sb-ext:*debug-print-variable-alist*
'((*print-array* . nil) (*print-length* . 10)))

Next:

backtrace

shows active calls and (if available) their arguments

--
Waldek Hebisch

Alasdair McAndrew

unread,

Oct 31, 2015, 5:41:24 AM10/31/15

to fricas...@googlegroups.com

Oops - yes you are quite right! The function integration-qag only returns two values, not three.

Thanks very much!

-Alasdair

Bill Page

unread,

Nov 2, 2015, 12:46:26 AM11/2/15

to fricas-devel

Alasdair, Kurt and Waldek,

I think I have found a reasonable way to avoid making too many
assumptions in the Lisp code concerning the representation of data.
The basic idea is the pass a function written in Spad to each Lisp
wrapper function. I refer to this function as 'per' after the macros
rep and per introduced by Aldor and used in some places in the Spad
code for type-safe conversions from to and from internal
representations.

Please see: https://github.com/billpage/gsla

I have inserted some of the relevant code below.

; gsl.lisp:
...
; The first parameter of each routine below, is a function 'per' written
; in Spad which must convert the internal representation of returned
; values to an appropriate external representation. The purpose of 'per'
; is to minimize the number of assumptions in the Lisp code concerning the
; representation of values.
; The internal representation is that obtained directly from GSLL and
; modified if necessary to be either atomic and compatible with Spad or a
; List of such returned values. When the output of GSL is multivalued the
; output of 'per' is usually of type Record.
; The rest of the parameters must be either atomic and compatible with Spad
; or a List of such input values. The input values are modified if necessary
; to be compatible with GSLL.

(defun gslintegrationqng (per f a b)
(apply (|mkLispFunction3| per) (multiple-value-list
(gsl:integration-qng (|mkLispFunction1| f) a b))))
...

-- gsl.spad
...
integrationQng: (DF -> DF,DF,DF) -> Record(result:DF, abserr:DF,
neval:Integer)
++ \spad{\integrationQng} applies the Gauss-Kronrod 10-point,
++ 21-point, 43-point and 87-point integration rules in succession
...
integrationQng(f,a,b) == GSLINTEGRATIONQNG(
-- convert returned values to Spad representation
(v1:DF,v2:DF,v3:Integer):Record(result:DF, abserr:DF,
neval:Integer)+->[v1,v2,v3],
f,a,b)$Lisp

Bill.

> Bill Page wrote:
>> This works:
>>
>> ; gsl.lisp
>> ...
>> ; Return multiple values as a vector (can be interpreted as Axiom Record)
>> (defun integration-qng-vector (f a b)
>> (apply #'vector (multiple-value-list (integration-qng f a b))))
>>
>

On 29 October 2015 at 12:59, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
>
> Note: this is for records of length > 2. Records of lennght 2 are
> represented as Lisp pairs, the same as unions. In fact, without
> knowing Spad type it is impossible to distinguish record from
> union:

> ...

> Yes. Let me just remark that the code above is fine for

> experimentation. In long term we should not sprinkle

Alasdair McAndrew

unread,

Nov 2, 2015, 6:51:37 PM11/2/15

to fricas...@googlegroups.com

Very nice! I have been trying to extend your (Bill's) original programs (with "pretend Integer") but these seem much nicer - although I don't fully understand how they work (nothing like starting from the bottom, is there?). You can also implement qagp - which requires a user-given list of singularities - with something like:

(defun gslintegration gagp (per f a b S)

(let* A (grid:make-foreign-array 'double-float :dimensions (length S) :initial-contents S))

(apply (|mkLispFunction2| per)
(multiple-value-list (gsl:integration-qagp (|mkLispFunction1| f) a b A))))

I haven't tested this yet, but I've tested its components. You can also include optional parameters (abserr, relerr etc). By the way, where is "per" defined?

-Alasdair

Bill Page

unread,

Nov 2, 2015, 8:43:32 PM11/2/15

to fricas-devel

On 2 November 2015 at 18:51, Alasdair McAndrew <amc...@gmail.com> wrote:
>
> Very nice! I have been trying to extend your (Bill's) original programs (with
> "pretend Integer") but these seem much nicer -

Although there is no explicit use of 'pretend' there is still some
implicit "pretending" going on behind the scenes. The main point
however is to try to do most of the work in Spad without pretend so
that things are as type-safe and representation independent as
possible.

> although I don't fully understand how they work (nothing like starting from
> the bottom, is there?).

I very much appreciate every time you ask a question...

> You can also implement qagp - which requires a user-given list of
> singularities - with something like:
>
> (defun gslintegration gagp (per f a b S)
> (let* A (grid:make-foreign-array 'double-float :dimensions (length S) :initial-contents S))
> (apply (|mkLispFunction2| per)
> (multiple-value-list (gsl:integration-qagp (|mkLispFunction1| f) a b A))))
>

Please is the slightly revised code for integrationQagp at
https://github.com/billpage/gsla

> I haven't tested this yet, but I've tested its components. You can also include optional
> parameters (abserr, relerr etc).

Can I leave that for you to do based on the example that I sent in
another email? :)

> By the way, where is "per" defined?
>

The function that is passed to the parameter 'per' is defined in the Spad code.

Cheers,
Bill Page.

Alasdair McAndrew

unread,

Nov 2, 2015, 11:47:34 PM11/2/15

to fricas...@googlegroups.com

Yep:

integrationQagp(x+->x^3*log(abs((x^2-1)*(x^2-2))),[0.0,1.0,sqrt(2.0),3.0])

returns the answer

[result= 52.74074838347258,abserr= 1.3130161420349395e-6]
Type: Record(result: DoubleFloat,abserr: DoubleFloat)

which is correct. Note that the absolute error is set in GSLL to be 1.0e-5, which is quite large (but on the other hand almost always ensures convergence), so unless you allow the user to enter an error value, the result will only be accurate to that value.

-Alasdair

Ralf Hemmecke

unread,

Nov 3, 2015, 1:34:27 AM11/3/15

to fricas...@googlegroups.com

On 11/03/2015 12:51 AM, Alasdair McAndrew wrote:
> By the way, where is "per" defined?

>> (defun gslintegrationqng (per f a b)
>> (apply (|mkLispFunction3| per) (multiple-value-list
>> (gsl:integration-qng (|mkLispFunction1| f) a b))))
>> ...
>>
>> -- gsl.spad
>> ...
>> integrationQng: (DF -> DF,DF,DF) -> Record(result:DF, abserr:DF,
>> neval:Integer)
>> ++ \spad{\integrationQng} applies the Gauss-Kronrod 10-point,
>> ++ 21-point, 43-point and 87-point integration rules in succession
>> ...
>> integrationQng(f,a,b) == GSLINTEGRATIONQNG(
>> -- convert returned values to Spad representation
>> (v1:DF,v2:DF,v3:Integer):Record(result:DF, abserr:DF,
>> neval:Integer)+->[v1,v2,v3],
>> f,a,b)$Lisp

As far as I can see, the definition of integrationQng calls
GSLINTEGRATIONQNG with 4 parameters. The first parameter is the per,
i.e. the (anonymous) function:

(v1:DF,v2:DF,v3:Integer): Record(result:DF, abserr:DF,neval:Integer) _
+-> [v1,v2,v3]

Ralf

Alasdair McAndrew

unread,

Nov 3, 2015, 8:03:27 PM11/3/15

to fricas...@googlegroups.com

A few thoughts:

How do we allow for optional parameters (epsabs, epsrel, limit etc)? The approach taken by GSLL, and therefore by me in my version of gsl.spad was to overload the integration function, and to give values for the extra parameters when they weren't given by the user. This all had to be done in order: epsabs, epsrel, limit. This means that you can't set the limit without setting epsabs and epsrel first. Is there a "named" way of entering optional parameters in spad: integrationQng(x+->sin(x^2),0.0,%pi/2,'limit=2000) say? This means that the user can set any optional parameters, in any order.
Improper integrals: GSLL has three separate functions: qagi (integrating from -inf to +inf), qagiu (integrating from a to inf, a being finite), and qagil (integrating from -inf to b). For each function it makes a substitution so that the integral is defined between 0 and 1 - you can see the substitutions in the gsl source file integration/qags.c. Thus the GSLL tests for qagi, qagiu, and qagil independently. The approach taken in Maxima is to have just one function, which they call quad_qagi, which can be used for any integral involving infinity, and which makes the appropriate substitution depending on which of the limits is finite.
Names of functions: gslIntegrationQng, integrationQng, quadQng...? (I like the last because it's short!)
Other GSL/GSLL functions... I suppose some numerical linear algebra (including eigensystems) would be nice, random number generation, and numerical ODE solving - GSL supports a host of different methods, mainly embedded Runge-Kutta methods of different orders. I think that FFTs and other transforms (wavelets) are better suited for purely numerical software, like GNU Octave. Also, the GSL implementation of FFT is not the most optimized: if we wanted one we should go for FFTW ("Fastest Fourier Transform in the West").

Oooh look - it's time to have a cup of tea!

-Alasdair

Ralf Hemmecke

unread,

Nov 4, 2015, 3:53:18 AM11/4/15

to fricas...@googlegroups.com

Well, if it were Aldor, you would have "keyword arguments".
See Section 6.3 of http://www.aldor.org/docs/aldorug.pdf .
I don't think that works in SPAD (yet).

However, what does the underlying GSL actually do? Do they allow
optional parameters? Or have they different function(name)s?

I actually would try to map the GSL function name to SPAD with the same
argument list.

> 2. Improper integrals: GSLL has three separate functions: qagi

> (integrating from -inf to +inf), qagiu (integrating from a to inf, a being
> finite), and qagil (integrating from -inf to b). For each function it
> makes a substitution so that the integral is defined between 0 and 1 - you
> can see the substitutions in the gsl source file integration/qags.c. Thus
> the GSLL tests for qagi, qagiu, and qagil independently. The approach
> taken in Maxima is to have just one function, which they call quad_qagi,
> which can be used for any integral involving infinity, and which makes the
> appropriate substitution depending on which of the limits is finite.

If you look at http://fricas.github.io/api/Infinity.html you find the
domain constructor OrderedCompletion. I bet you want OrderedCompletion
DoubleFloat (let's abbreviat that by OCDF).

So you would allow OCDF instead of just DF for the limits. Or maybe you
even want something like Segment OCDF, then it would be possible to write

foo(f, minusInfinity()..plusInfinity)
foo(f, 5.3..plusInfinity)
foo(f, 1.1..2.7)

for the integral bounds.

Inside that function you take what OrderedCompletion provides you to
single out finite or infinite bounds and call the respective gsll function.

finite? http://fricas.github.io/api/OrderedCompletion.html#index-1
infinite? http://fricas.github.io/api/OrderedCompletion.html#index-2
whatInfinity http://fricas.github.io/api/OrderedCompletion.html#index-8

> 3. Names of functions: gslIntegrationQng, integrationQng, quadQng...?

> (I like the last because it's short!)

Well, if quadQng is known and common in GSL, why not. Otherwise, you
should think of how a user would find your function in FriCAS. I guess,
it's never a bad idea to have 'integrate' or something similar in the
name, if the function deals with that. Don't worry about long names. The
user has the macro system to deal with too much typing.

Ralf

Bill Page

unread,

Nov 4, 2015, 7:58:05 AM11/4/15

to fricas-devel

On 3 November 2015 at 20:03, Alasdair McAndrew <amc...@gmail.com> wrote:
>
> A few thoughts:
>
> How do we allow for optional parameters (epsabs, epsrel, limit etc)? The approach taken by GSLL, and therefore by me in my version of gsl.spad was to overload the integration function, and to give values for the extra parameters when they weren't given by the user. This all had to be done in order: epsabs, epsrel, limit. This means that you can't set the limit without setting epsabs and epsrel first. Is there a "named" way of entering optional parameters in spad: integrationQng(x+->sin(x^2),0.0,%pi/2,'limit=2000) say? This means that the user can set any optional parameters, in any order.

The primary example of something like this is the options for the
'draw' command. You can see the source code in
fricas/src/algebra/drawopt.spad and how it is used in draw.spad. But
a simplified form of this is just

(1) -> vars:=OrderedVariableList [epsabs,epsrel,limit]

(1) OrderedVariableList([epsabs,epsrel,limit])
Type: Type
(2) -> opts:List Record(key:vars,val:Any)
Type: Void
(3) -> opts:=[[epsabs,1.0E-20],[limit,1000]]

(3) [[key= epsabs,val= 0.1 E -19],[key= limit,val= 1000]]
Type: List(Record(key: OrderedVariableList([epsabs,epsrel,limit]),val: Any))

(4) -> opts:=[[notakey,3]]

Cannot convert the value from type Variable(notakey) to
OrderedVariableList([epsabs,epsrel,limit]) .

Using this the syntax for options could be

integrationQng(x+->sin(x^2),0.0,%pi/2,[[limit,2000]])

The outter [ ] forms a list and the inner [option,value] forms a
record. Each function would have to be defined twice - once with
options and once without (to avoid having to write [[]] in each
function). To process the options one must scan the list for specific
keys.

(5) -> for opt in opts repeat output [lookup opt.key,opt.val]
[1,0.1 E -19]
[3,1000]
Type: Void

For more control and a slightly better syntax one can replace this
with a specific domain using this representation similar to DrawOpts.

Of course underneath at the Lisp/GSLL level you would have to pass all
the "optional" parameters and provide your own defaults.

> Improper integrals: GSLL has three separate functions: qagi (integrating from -inf to +inf), qagiu (integrating from a to inf, a being finite), and qagil (integrating from -inf to b). For each function it makes a substitution so that the integral is defined between 0 and 1 - you can see the substitutions in the gsl source file integration/qags.c. Thus the GSLL tests for qagi, qagiu, and qagil independently. The approach taken in Maxima is to have just one function, which they call quad_qagi, which can be used for any integral involving infinity, and which makes the appropriate substitution depending on which of the limits is finite.

I think maybe Ralf's proposal is a reasonable one. E.g.

(6) -> a:OrderedCompletion DoubleFloat
Type: Void
(7) -> a:=plusInfinity()

(7) + infinity
Type: OrderedCompletion(DoubleFloat)
(8) -> b:OrderedCompletion DoubleFloat
Type: Void
(9) -> b:=3.141592

(9) 3.141592
Type: OrderedCompletion(DoubleFloat)
(10) -> retract(b)@DoubleFloat

(10) 3.141592
Type: DoubleFloat
(11) -> b:=%pi/2

(11) 1.5707963267948966
Type: OrderedCompletion(DoubleFloat)

So then

integrationQng(f,0. 0::DF, plusInfinity())

> Names of functions: gslIntegrationQng, integrationQng, quadQng...? (I like the last because it's short!)

Actually as we implement more and more integrations functions it seems
reasonable to me that there be a package specifically for integration,
e.g. GSLintegration. Then the names could appear as qng$GSLintegration
or just 'qng' in the appropriate context.

I would also like to get rid of the requirement to pass a function and
allow just an Expression, e.g.

qng(sin(x^2), 0.0,%pi/2,[[limit,2000]])$GSLintegration

> Other GSL/GSLL functions... I suppose some numerical linear algebra (including eigensystems) would be nice, random number generation, and numerical ODE solving - GSL supports a host of different methods, mainly embedded Runge-Kutta methods of different orders. I think that FFTs and other transforms (wavelets) are better suited for purely numerical software, like GNU Octave. Also, the GSL implementation of FFT is not the most optimized: if we wanted one we should go for FFTW ("Fastest Fourier Transform in the West").
>

I have an application where I would like to use a good ODE solver. I
will work on that.

I was also looking into other interesting packages available in
QuickLisp and which might be amenable to the same interface approach
as GSLL. One thing that stood out was graphics such as 'vgplot' as an
interface to GnuPlot

https://github.com/volkers/vgplot

and also Cl-ppplot which is an interface to the even more capable

http://plplot.sourceforge.net/

Cheers,
Bill.

Mark Clements

unread,

Nov 12, 2015, 3:01:38 PM11/12/15

to fricas...@googlegroups.com

A further change could be to integrate for an Expression using a
SegmentBinding:

)co gsl
EF ==> Expression Float
DF ==> DoubleFloat
integrationQng(expr : EF, sb : SegmentBinding(EF)) :
Record(result: DF, abserr: DF, neval: Integer) ==
(var, a, b) := (variable sb, lo segment sb, hi segment sb)

integrationQng(makeFloatFunction(expr,var)$MakeFloatCompiledFunction(EF),a::DF,b::DF)$Gsl
integrationQng(exp(-y^2),y=0..1)

This is arguably more idiomatic.

--- Mark

Waldek Hebisch

unread,

Nov 13, 2015, 11:16:23 AM11/13/15

to fricas...@googlegroups.com

Ralf Hemmecke wrote:
>
> Well, if it were Aldor, you would have "keyword arguments".
> See Section 6.3 of http://www.aldor.org/docs/aldorug.pdf .
> I don't think that works in SPAD (yet).

AFAIK this is quite different that Lisp "keyword arguments":
in particular in Lisp everything is done at runtime. I do
not think this Aldor construct would help calling Lisp
routines (of course re-coding Lisp routines in Aldor is
easier due to this construct).

> However, what does the underlying GSL actually do? Do they allow
> optional parameters? Or have they different function(name)s?

AFAIK GSL has interface in C. In C in principle one can use
"varags" functions to have arbitarily messy untyped interface.
But normal practice is to have fixed argument list with specified
types. Normal practice in C is to use special values or extra
flag arguments. For example error of 0 may mean default value.
Or extra argument which says if requested error is absolute
or relative. AFAICS the whole multiple values and keyword
business comes entirely from Lisp wrappers. I would guess
that C code takes pointer to memory location where it can
store error estimate and number of steps. Usual convention
is that null pointer means that corresponding value is not
needed. In C one can also return a structure, but this is
less usual.

--
Waldek Hebisch

Reply all

Reply to author

Forward