Interfacing C functions (with SBCL)

[Tobias Neumann]

I've noticed that FriCAS doesn't have a numerical implementation of polylog's yet.

So I will have to interface some C/Fortran code. For that I have constructed a small

example:

test.c:

void testfun(double arg, double* realpart, double* imagpart)

{

    *realpart = 1.0 * arg;

    *imagpart = 3.14 * arg;

    return;

}

gcc -shared -o libtest.so test.c

test.lsp:

(sb-alien:load-shared-object "./libtest.so")

(defun testfun (arg)

    (sb-alien:with-alien ((repart (double-float) 0.0d0) (impart (double-float) 0.0d0))

        (sb-alien:alien-funcall

          (sb-alien:extern-alien "testfun" (sb-alien:function sb-alien:void (double-float) (* double-float) (* double-float)))

          arg (sb-alien:addr repart) (sb-alien:addr impart)

         )

        (complex repart impart)

    )

)

And I can successfully use that function within sbcl:

sbcl --load "test.lsp"

(testfun 10.0d0)

#C(10.0d0 31.400000000000002d0)

Within FriCAS I ran ")lisp (load "test.lsp")" and it

seems to parse/load the file successfully, but I can't directly use that function, i.e.  (testfun 10.0d0) within ")fin", and neither from within FriCAS as TESTFUN(1$DoubleFloat)$Lisp:

debugger invoked on a TYPE-ERROR in thread

#<THREAD "main thread" RUNNING {1002410103}>:

  The value

    0.0d0

  is not of type

    (SB-ALIEN:ALIEN (* T))

  from the function type declaration.

As it turns out, somehow I can not give it an initial value here, and I just have to declare it with (repart (double-float)), which is fine, but I am wondering why.

The final step for me is to convert the SExpression into a FriCAS Complex(DoubleFloat). For that I looked into the SExpression domain, and it seems to be easiest to instead return (list repart impart) and then have something like: 

testfun(xx) == 

    res := TESTFUN(xx :: DoubleFloat)$Lisp

    complex(float car res, float car cdr res) 

Presumably I could convert this into Complex(DoubleFloat) within the Lisp code directly with some SPADCALL calls? But the current way seems to be perfectly fine, but I am wondering if this is the right/recommended way to do it. Overall I'm really happy how easy this turned out to be.

Thanks,

Tobias

[Qian Yun]

Hi,

I can not answer the first part of the question, about the initial value.

But for the second part, because in SPAD, the COMPLEX is
represented as
Rep := Record(real : R, imag : R)
which is essentially a Common Lisp cons pair, so you can do this directly:

testfun(xx) == TESTFUN(xx :: DoubleFloat)$Lisp pretend Complex(DoubleFloat)

AND use "(cons repart impart)" instead of "(complex repart impart)".

- Qian

> --
> 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
> <mailto:fricas-devel...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/fricas-devel/5eb8ad23-5c3e-4b0b-8810-9cccdc42d999n%40googlegroups.com
> <https://groups.google.com/d/msgid/fricas-devel/5eb8ad23-5c3e-4b0b-8810-9cccdc42d999n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Kurt Pagani]

The problem seems to be known for a long time. It works when "compiled" but not
when interpreted.

https://bugs.launchpad.net/sbcl/+bug/1731556
https://bugs.launchpad.net/sbcl/+bug/992362
https://bugs.launchpad.net/sbcl/+bug/734259


Note: I added a "(in-package :cl-user)" top in your test.lsp, otherwise you are
in package "BOOT" where symbol conflicts may occur, e.g. (coerce 0 'double-float):


FriCAS Computer Algebra System
Version: FriCAS 1.3.5
Timestamp: Mi, 17. Apr 2019 03:15:11
-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave FriCAS and return to shell.
-----------------------------------------------------------------------------

Function declaration sixel : TexFormat -> Void has been added to
workspace.
Type: Void
Function declaration quickLoad : String -> Void has been added to
workspace.
Type: Void
(3) -> )lisp (load "test.lsp")

Value = T
(3) -> )lisp (compile "test.lsp")

>> System error:
Invalid function name: "test.lsp"

(3) -> )lisp (compile testfun)

>> System error:
The variable TESTFUN is unbound.

(3) -> )lisp (compile 'testfun)

>> System error:
The function BOOT::TESTFUN is undefined.

(3) -> )lisp (compile 'cl-user::testfun)

Value = COMMON-LISP-USER::TESTFUN
(3) -> )lisp (testfun 10.0d0)

>> System error:
The function BOOT::TESTFUN is undefined.

(3) -> )lisp (CL-USER::testfun 10.0d0)

Value = #C(10.0 31.400000000000002)
(3) ->



This link may ALSO be helpful:
https://stackoverflow.com/questions/21353376/how-to-use-double-float

Although I find it realy ghastly interfacing C in fricas (but if, why not using
CFFI?)

> The value
> 0.0d0
> is not of type
> (SB-ALIEN:ALIEN (* T))

)lisp (defun |foo| (x) x)
foo(7)$Lisp -> 7
foo(7.2)$Lisp -> (265633114661417543270 . - 65)
foo(7.2::DoubleFloat)$Lisp -> 7.199999999999999

)lisp (coerce 0 'double-float)
)lisp (defun mycoerce (x) (cl-user::coerce x 'double-float))

[Waldek Hebisch]

On Mon, Apr 12, 2021 at 06:03:34PM -0700, Tobias Neumann wrote:
> I've noticed that FriCAS doesn't have a numerical implementation of
> polylog's yet.

I wonder what do you need. I have partial multiple precision
implementation, which does not cover full argument range but
seem to work well for some arguments.

Concerning calling to C/Fortran, maybe the following will
help. First, modification of your C file:

---------------<test.c cut here>-----------------

double
testfun1(double arg) {
return 3.14 * arg;
}

void
testfun2(double arg, double* realpart, double* imagpart) {

*realpart = 1.0 * arg;
*imagpart = 3.14 * arg;
}

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

Next, Lisp wrapper:

---------------<tst.lisp cut here>-----------------

(fricas-lisp::fricas-foreign-call |testfun1| "testfun1" fricas-lisp::double
(arg fricas-lisp::double))

(defun |testfun2| (arg)
(let* (
(tmp-bufr (sb-alien:make-alien sb-alien:double 1))
(tmp-bufi (sb-alien:make-alien sb-alien:double 1)))
(sb-alien:alien-funcall
(sb-alien:extern-alien "testfun2"
(sb-alien:function sb-alien:void
SB-ALIEN::double (* t) (* t)))
arg tmp-bufr tmp-bufi)
(cons (sb-alien:deref tmp-bufr 0)
(sb-alien:deref tmp-bufi 0))

)
)

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

C compilation:

gcc -Wall -O -shared -o libtest.so test.c

In FriCAS:

(1) -> )lisp (sb-alien::load-shared-object "./libtest.so")
Value = #P"./libtest.so"
(1) -> )lisp (load (compile-file "tst.lisp"))

; compiling file "/mnt/lv3/fricas/axp19/pp1/pp19/tst.lisp" (written 16 APR 2021 01:47:09 PM):

; wrote /mnt/lv3/fricas/axp19/pp1/pp19/tst.fasl
; compilation finished in 0:00:00.012
Value = T
(1) -> valr := testfun1(10.0::DoubleFloat)$Lisp pretend DoubleFloat

(1) 31.400000000000002
Type: DoubleFloat
(2) -> valc := testfun2(10.0::DoubleFloat)$Lisp pretend Complex(DoubleFloat)

(2) 10.0 + 31.400000000000002 %i
Type: Complex(DoubleFloat)

Comments:
- I am not sure what matters there, maybe fact that I compiled the code,
maybe direct use of low-level sbcl constructs
- when arguments are passed by value 'fricas-foreign-call' allows
relatively easy definition. ATM 'fricas-foreign-call' seem to
miss support for C 'long' but that could be easily added. Also,
needed symbols are not exported. For general use we probably
should rename it to 'fricas_foreign_call' and make names
lowercase
- pointer arguments are tricky. When doing call we may need to
allocate memory, but for performance we prefer passing addresses
of existing FriCAS data. We may need extra copies, for Complex
we need to cons wrapper node
- extra complication are caused by C structs

--
Waldek Hebisch