From 4e0db4ae6a1f2e9592f33847889a112eb21f285a Mon Sep 17 00:00:00 2001
From: Blaise Bourdin <bourdin@mcmaster.ca>
Date: Mon, 19 Oct 2026 11:40:16 -0400
Subject: [PATCH 1/3] Require CMake 3.10

CMake 4 no longer accepts cmake_minimum_required versions older than 3.5.
---
 CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1bdc2ff..2dd6f52 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 2.8)
+cmake_minimum_required(VERSION 3.10)
 enable_language(C Fortran)
 project(symengine.f90)
 set (CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/mod")

From 3d88852dd7196a223a3b464842337bc28565e880 Mon Sep 17 00:00:00 2001
From: Blaise Bourdin <bourdin@mcmaster.ca>
Date: Mon, 19 Oct 2026 11:40:16 -0400
Subject: [PATCH 2/3] Add Lambdify

Lambdify turns expressions into numerical functions, like Lambdify in
symengine.py, using the lambda_real_double_visitor and llvm_double_visitor
C API:

    call f%init([ptr(x), ptr(y)], [ptr(x**2 + y**2)], cse=.true., backend='llvm')
    call f%eval(inps, outs)   ! inps(nin, npts) -> outs(nout, npts)
    call f%get_function(fn)   ! LLVM compiled function as a procedure pointer

The LLVM backend is compiled in when SymEngine provides it
(HAVE_SYMENGINE_LLVM). Requires the get_num_inputs/outputs, call_batch and
get_function functions of the SymEngine C wrapper.
---
 README.md                   |  27 ++++
 src/CMakeLists.txt          |   4 +
 src/lambdify.F90            | 276 ++++++++++++++++++++++++++++++++++++
 src/symengine.f90           |   2 +
 src/symengine_interface.f90 |  71 +++++++++-
 test/CMakeLists.txt         |   7 +
 test/test_lambdify.F90      | 131 +++++++++++++++++
 7 files changed, 517 insertions(+), 1 deletion(-)
 create mode 100644 src/lambdify.F90
 create mode 100644 test/test_lambdify.F90

diff --git a/README.md b/README.md
index e6812cb..4fd1696 100644
--- a/README.md
+++ b/README.md
@@ -15,3 +15,30 @@ Install SymEngine.f90:
     cmake -DCMAKE_PREFIX_PATH=/some/path .
     make
     make install
+
+## Lambdify
+
+`Lambdify` turns expressions into numerical functions, like `Lambdify` in
+symengine.py:
+
+```fortran
+use symengine
+type(Basic) :: x, y
+type(Lambdify) :: f
+real(c_double) :: out(1), inps(2, 100), outs(1, 100)
+procedure(lambdify_function), pointer :: fn
+
+x = Symbol('x')
+y = Symbol('y')
+call f%init([ptr(x), ptr(y)], [ptr(x**2 + y**2)])   ! or arrays of type(Basic)
+call f%eval([1.0d0, 2.0d0], out)                   ! one point
+call f%eval(inps, outs)                            ! inps(2, n) -> outs(1, n)
+```
+
+Optional arguments of `init`: `cse` (common subexpression elimination),
+`backend` (`'lambda'`, the default, or `'llvm'`) and `opt_level` (LLVM only).
+With `backend='llvm'`, `call f%get_function(fn)` returns the compiled function
+as a procedure pointer, called as `call fn(inps, outs)`. The LLVM backend
+requires SymEngine built with `WITH_LLVM=ON`; when building with fpm, also pass
+`--flag -DHAVE_SYMENGINE_LLVM`. A `Lambdify` object cannot be copied, and with
+the `'lambda'` backend each thread needs its own object.
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ccca0fc..a929563 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -9,9 +9,13 @@ set(SRC
     constants.f90
     dense_matrix.f90
     sets.f90
+    lambdify.F90
     symengine.f90
     )
 add_library(symengine_f90 ${SRC})
 target_link_libraries(symengine_f90 ${SYMENGINE_LIBRARIES})
+if (HAVE_SYMENGINE_LLVM)
+    target_compile_definitions(symengine_f90 PRIVATE HAVE_SYMENGINE_LLVM)
+endif()
 
 install(TARGETS symengine_f90 DESTINATION lib)
diff --git a/src/lambdify.F90 b/src/lambdify.F90
new file mode 100644
index 0000000..d095aa0
--- /dev/null
+++ b/src/lambdify.F90
@@ -0,0 +1,276 @@
+module symengine_lambdify
+
+! Turns symbolic expressions into numerical functions, like Lambdify in
+! symengine.py:
+!
+!     type(Lambdify) :: f
+!     call f%init([ptr(x), ptr(y)], [ptr(x**2 + y**2)])
+!     call f%eval([1.0d0, 2.0d0], outs)          ! one point
+!     call f%eval(inps, outs)                    ! inps(2, n), outs(1, n)
+!
+! Two backends are available:
+!  - 'lambda' (default): a tree of closures, always available.
+!  - 'llvm': machine code generated by the LLVM JIT, only if SymEngine was
+!    built with LLVM. get_function returns it as a procedure pointer.
+!
+! A Lambdify object owns a C++ object and cannot be copied. With the 'lambda'
+! backend, use one object per thread.
+
+use iso_c_binding, only: c_int, c_long, c_double, c_size_t, c_ptr, c_null_ptr, c_funptr, &
+                         c_associated, c_f_procpointer
+use symengine_interface
+use symengine_basic
+use exceptions
+
+implicit none
+
+integer, parameter :: BACKEND_NONE = 0, BACKEND_LAMBDA = 1, BACKEND_LLVM = 2
+
+type Lambdify
+    type(c_ptr), private :: ptr = c_null_ptr
+    integer, private :: backend = BACKEND_NONE
+contains
+    procedure :: lambdify_init_ptr, lambdify_init_basic
+    generic :: init => lambdify_init_ptr, lambdify_init_basic
+    procedure :: num_inputs => lambdify_num_inputs
+    procedure :: num_outputs => lambdify_num_outputs
+    procedure :: lambdify_eval_point, lambdify_eval_batch
+    generic :: eval => lambdify_eval_point, lambdify_eval_batch
+    procedure :: get_function => lambdify_get_function
+    procedure :: free => lambdify_free
+    procedure, private :: lambdify_assign
+    generic :: assignment(=) => lambdify_assign
+    final :: lambdify_final
+end type
+
+! Interface of the function returned by get_function
+abstract interface
+    subroutine lambdify_function(inps, outs) bind(c)
+        import :: c_double
+        real(c_double), dimension(*), intent(in) :: inps
+        real(c_double), dimension(*), intent(out) :: outs
+    end subroutine
+end interface
+
+private
+public :: Lambdify, lambdify_function
+
+contains
+
+subroutine lambdify_init_ptr(this, args, exprs, cse, backend, opt_level)
+    class(Lambdify), intent(inout) :: this
+    type(c_ptr), dimension(:), intent(in) :: args, exprs
+    logical, intent(in), optional :: cse
+    character(len=*), intent(in), optional :: backend
+    integer, intent(in), optional :: opt_level
+
+    type(c_ptr) :: vargs, vexprs
+    integer(c_int) :: perform_cse, level, exception
+    integer :: i
+
+    call this%free()
+
+    exception = 0
+    perform_cse = 0
+    if (present(cse)) then
+        if (cse) perform_cse = 1
+    end if
+    level = 3
+    if (present(opt_level)) level = int(opt_level, c_int)
+
+    this%backend = BACKEND_LAMBDA
+    if (present(backend)) then
+        select case (backend)
+        case ('lambda')
+            this%backend = BACKEND_LAMBDA
+        case ('llvm')
+#ifdef HAVE_SYMENGINE_LLVM
+            this%backend = BACKEND_LLVM
+#else
+            error stop "Lambdify: SymEngine was built without LLVM support"
+#endif
+        case default
+            error stop "Lambdify: backend must be 'lambda' or 'llvm'"
+        end select
+    end if
+
+    vargs = c_vecbasic_new()
+    vexprs = c_vecbasic_new()
+    do i = 1, size(args)
+        call handle_exception(c_vecbasic_push_back(vargs, args(i)))
+    end do
+    do i = 1, size(exprs)
+        call handle_exception(c_vecbasic_push_back(vexprs, exprs(i)))
+    end do
+
+    select case (this%backend)
+    case (BACKEND_LAMBDA)
+        this%ptr = c_lambda_real_double_visitor_new()
+        exception = c_lambda_real_double_visitor_init(this%ptr, vargs, vexprs, perform_cse)
+#ifdef HAVE_SYMENGINE_LLVM
+    case (BACKEND_LLVM)
+        this%ptr = c_llvm_double_visitor_new()
+        exception = c_llvm_double_visitor_init(this%ptr, vargs, vexprs, perform_cse, level)
+#endif
+    end select
+
+    call c_vecbasic_free(vargs)
+    call c_vecbasic_free(vexprs)
+    if (exception /= 0) then
+        call this%free()
+        call handle_exception(int(exception, c_long))
+    end if
+end subroutine
+
+subroutine lambdify_init_basic(this, args, exprs, cse, backend, opt_level)
+    class(Lambdify), intent(inout) :: this
+    class(Basic), dimension(:), intent(in) :: args, exprs
+    logical, intent(in), optional :: cse
+    character(len=*), intent(in), optional :: backend
+    integer, intent(in), optional :: opt_level
+
+    type(c_ptr), dimension(size(args)) :: pargs
+    type(c_ptr), dimension(size(exprs)) :: pexprs
+    integer :: i
+
+    do i = 1, size(args)
+        pargs(i) = args(i)%ptr
+    end do
+    do i = 1, size(exprs)
+        pexprs(i) = exprs(i)%ptr
+    end do
+    call this%lambdify_init_ptr(pargs, pexprs, cse, backend, opt_level)
+end subroutine
+
+function lambdify_num_inputs(this) result(res)
+    class(Lambdify), intent(in) :: this
+    integer :: res
+
+    res = 0
+    select case (this%backend)
+    case (BACKEND_LAMBDA)
+        res = int(c_lambda_real_double_visitor_get_num_inputs(this%ptr))
+#ifdef HAVE_SYMENGINE_LLVM
+    case (BACKEND_LLVM)
+        res = int(c_llvm_double_visitor_get_num_inputs(this%ptr))
+#endif
+    end select
+end function
+
+function lambdify_num_outputs(this) result(res)
+    class(Lambdify), intent(in) :: this
+    integer :: res
+
+    res = 0
+    select case (this%backend)
+    case (BACKEND_LAMBDA)
+        res = int(c_lambda_real_double_visitor_get_num_outputs(this%ptr))
+#ifdef HAVE_SYMENGINE_LLVM
+    case (BACKEND_LLVM)
+        res = int(c_llvm_double_visitor_get_num_outputs(this%ptr))
+#endif
+    end select
+end function
+
+! Evaluates at npts points. inps and outs are stored point after point.
+subroutine lambdify_call(this, inps, outs, npts)
+    class(Lambdify), intent(inout) :: this
+    real(c_double), dimension(*), intent(in) :: inps
+    real(c_double), dimension(*), intent(out) :: outs
+    integer, intent(in) :: npts
+
+    select case (this%backend)
+    case (BACKEND_LAMBDA)
+        call c_lambda_real_double_visitor_call_batch(this%ptr, outs, inps, int(npts, c_size_t))
+#ifdef HAVE_SYMENGINE_LLVM
+    case (BACKEND_LLVM)
+        call c_llvm_double_visitor_call_batch(this%ptr, outs, inps, int(npts, c_size_t))
+#endif
+    case default
+        error stop "Lambdify: eval called before init"
+    end select
+end subroutine
+
+! Evaluates at one point: inps(num_inputs), outs(num_outputs)
+subroutine lambdify_eval_point(this, inps, outs)
+    class(Lambdify), intent(inout) :: this
+    real(c_double), dimension(:), intent(in) :: inps
+    real(c_double), dimension(:), intent(out) :: outs
+
+    integer :: nin, nout
+
+    if (this%backend == BACKEND_NONE) error stop "Lambdify: eval called before init"
+    nin = this%num_inputs()
+    nout = this%num_outputs()
+    if (size(inps) /= nin .or. size(outs) /= nout) then
+        error stop "Lambdify: eval called with arrays of the wrong size"
+    end if
+    call lambdify_call(this, inps, outs, 1)
+end subroutine
+
+! Evaluates at several points: inps(num_inputs, npts), outs(num_outputs, npts)
+subroutine lambdify_eval_batch(this, inps, outs)
+    class(Lambdify), intent(inout) :: this
+    real(c_double), dimension(:, :), intent(in) :: inps
+    real(c_double), dimension(:, :), intent(out) :: outs
+
+    integer :: nin, nout
+
+    if (this%backend == BACKEND_NONE) error stop "Lambdify: eval called before init"
+    nin = this%num_inputs()
+    nout = this%num_outputs()
+    if (size(inps, 1) /= nin .or. size(outs, 1) /= nout .or. size(inps, 2) /= size(outs, 2)) then
+        error stop "Lambdify: eval called with arrays of the wrong shape"
+    end if
+    call lambdify_call(this, inps, outs, size(inps, 2))
+end subroutine
+
+! Returns the JIT compiled function (backend='llvm' only). It can be called
+! concurrently from several threads, and is valid as long as this object is
+! alive and not re-initialized.
+subroutine lambdify_get_function(this, fn)
+    class(Lambdify), intent(in) :: this
+    procedure(lambdify_function), pointer, intent(out) :: fn
+#ifdef HAVE_SYMENGINE_LLVM
+    type(c_funptr) :: cfn
+#endif
+
+    nullify(fn)
+    if (this%backend /= BACKEND_LLVM) then
+        error stop "Lambdify: get_function requires backend='llvm'"
+    end if
+#ifdef HAVE_SYMENGINE_LLVM
+    cfn = c_llvm_double_visitor_get_function(this%ptr)
+    call c_f_procpointer(cfn, fn)
+#endif
+end subroutine
+
+subroutine lambdify_free(this)
+    class(Lambdify), intent(inout) :: this
+
+    if (c_associated(this%ptr)) then
+        select case (this%backend)
+        case (BACKEND_LAMBDA)
+            call c_lambda_real_double_visitor_free(this%ptr)
+#ifdef HAVE_SYMENGINE_LLVM
+        case (BACKEND_LLVM)
+            call c_llvm_double_visitor_free(this%ptr)
+#endif
+        end select
+    end if
+    this%ptr = c_null_ptr
+    this%backend = BACKEND_NONE
+end subroutine
+
+subroutine lambdify_final(this)
+    type(Lambdify), intent(inout) :: this
+    call this%free()
+end subroutine
+
+subroutine lambdify_assign(a, b)
+    class(Lambdify), intent(inout) :: a
+    class(Lambdify), intent(in) :: b
+    error stop "Lambdify: objects cannot be copied, call init instead"
+end subroutine
+
+end module
diff --git a/src/symengine.f90 b/src/symengine.f90
index c1c8d26..bf8b5d5 100644
--- a/src/symengine.f90
+++ b/src/symengine.f90
@@ -11,6 +11,7 @@ module symengine
 use constants
 use dense_matrix
 use sets
+use symengine_lambdify
 
 implicit none
 
@@ -27,6 +28,7 @@ module symengine
 public :: DenseMatrix, transpose, ones, zeros, eye
 public :: emptyset, universalset, complexes, reals, rationals, integers, set_union, set_intersection
 public :: interval, finiteset
+public :: Lambdify, lambdify_function
 
 contains
 
diff --git a/src/symengine_interface.f90 b/src/symengine_interface.f90
index bddb9c2..968b63b 100644
--- a/src/symengine_interface.f90
+++ b/src/symengine_interface.f90
@@ -1,6 +1,6 @@
 module symengine_interface
 
-use iso_c_binding, only: c_size_t, c_ptr, c_long, c_int, c_char, c_double, c_f_pointer
+use iso_c_binding, only: c_size_t, c_ptr, c_long, c_int, c_char, c_double, c_f_pointer, c_funptr
 implicit none
 
 interface
@@ -415,6 +415,75 @@ integer(c_long) function c_basic_set_finiteset(s, container) bind(c, name='basic
         import :: c_long, c_ptr
         type(c_ptr), value :: s, container
     end function
+    ! Lambdify. The C functions return a symengine_exceptions_t, i.e. a C int.
+    type(c_ptr) function c_lambda_real_double_visitor_new() bind(c, name='lambda_real_double_visitor_new')
+        import :: c_ptr
+    end function
+    integer(c_int) function c_lambda_real_double_visitor_init(self, args, exprs, perform_cse) &
+            bind(c, name='lambda_real_double_visitor_init')
+        import :: c_int, c_ptr
+        type(c_ptr), value :: self, args, exprs
+        integer(c_int), value :: perform_cse
+    end function
+    integer(c_size_t) function c_lambda_real_double_visitor_get_num_inputs(self) &
+            bind(c, name='lambda_real_double_visitor_get_num_inputs')
+        import :: c_size_t, c_ptr
+        type(c_ptr), value :: self
+    end function
+    integer(c_size_t) function c_lambda_real_double_visitor_get_num_outputs(self) &
+            bind(c, name='lambda_real_double_visitor_get_num_outputs')
+        import :: c_size_t, c_ptr
+        type(c_ptr), value :: self
+    end function
+    subroutine c_lambda_real_double_visitor_call_batch(self, outs, inps, npts) &
+            bind(c, name='lambda_real_double_visitor_call_batch')
+        import :: c_ptr, c_double, c_size_t
+        type(c_ptr), value :: self
+        real(c_double), dimension(*) :: outs
+        real(c_double), dimension(*), intent(in) :: inps
+        integer(c_size_t), value :: npts
+    end subroutine
+    subroutine c_lambda_real_double_visitor_free(self) bind(c, name='lambda_real_double_visitor_free')
+        import :: c_ptr
+        type(c_ptr), value :: self
+    end subroutine
+    ! Only available if SymEngine was built with LLVM
+    type(c_ptr) function c_llvm_double_visitor_new() bind(c, name='llvm_double_visitor_new')
+        import :: c_ptr
+    end function
+    integer(c_int) function c_llvm_double_visitor_init(self, args, exprs, perform_cse, opt_level) &
+            bind(c, name='llvm_double_visitor_init')
+        import :: c_int, c_ptr
+        type(c_ptr), value :: self, args, exprs
+        integer(c_int), value :: perform_cse, opt_level
+    end function
+    integer(c_size_t) function c_llvm_double_visitor_get_num_inputs(self) &
+            bind(c, name='llvm_double_visitor_get_num_inputs')
+        import :: c_size_t, c_ptr
+        type(c_ptr), value :: self
+    end function
+    integer(c_size_t) function c_llvm_double_visitor_get_num_outputs(self) &
+            bind(c, name='llvm_double_visitor_get_num_outputs')
+        import :: c_size_t, c_ptr
+        type(c_ptr), value :: self
+    end function
+    subroutine c_llvm_double_visitor_call_batch(self, outs, inps, npts) &
+            bind(c, name='llvm_double_visitor_call_batch')
+        import :: c_ptr, c_double, c_size_t
+        type(c_ptr), value :: self
+        real(c_double), dimension(*) :: outs
+        real(c_double), dimension(*), intent(in) :: inps
+        integer(c_size_t), value :: npts
+    end subroutine
+    type(c_funptr) function c_llvm_double_visitor_get_function(self) &
+            bind(c, name='llvm_double_visitor_get_function')
+        import :: c_funptr, c_ptr
+        type(c_ptr), value :: self
+    end function
+    subroutine c_llvm_double_visitor_free(self) bind(c, name='llvm_double_visitor_free')
+        import :: c_ptr
+        type(c_ptr), value :: self
+    end subroutine
 end interface 
 
 contains
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 5dc82d3..df55e6b 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -9,3 +9,10 @@ endmacro(ADDTEST)
 project(tests)
 
 ADDTEST(test_basic)
+
+add_executable(test_lambdify test_lambdify.F90)
+target_link_libraries(test_lambdify symengine_f90)
+if (HAVE_SYMENGINE_LLVM)
+    target_compile_definitions(test_lambdify PRIVATE HAVE_SYMENGINE_LLVM)
+endif()
+add_test(test_lambdify ${PROJECT_BINARY_DIR}/test_lambdify)
diff --git a/test/test_lambdify.F90 b/test/test_lambdify.F90
new file mode 100644
index 0000000..a506218
--- /dev/null
+++ b/test/test_lambdify.F90
@@ -0,0 +1,131 @@
+program test_lambdify
+    use symengine
+    use iso_c_binding, only: c_double
+
+    implicit none
+
+    call test_point()
+    call test_batch()
+    call test_basic_array()
+#ifdef HAVE_SYMENGINE_LLVM
+    call test_llvm()
+#endif
+
+    print *, "Finishing"
+
+    contains
+
+    subroutine assert_close(a, b)
+        real(c_double), intent(in) :: a, b
+        if (abs(a - b) > 1e-12_c_double) then
+            print *, a, " /= ", b
+            error stop "Not equal"
+        end if
+    end subroutine
+
+    subroutine assert(a)
+        logical :: a
+        if (.not. a) then
+            error stop "Assertion failed"
+        end if
+    end subroutine
+
+    subroutine test_point()
+        type(Basic) :: x, y
+        type(Lambdify) :: f
+        real(c_double) :: outs(1)
+        real(c_double), parameter :: vals(3) = [-1.0_c_double, 0.0_c_double, 1.0_c_double]
+        integer :: c, i, j
+
+        x = Symbol('x')
+        y = Symbol('y')
+        do c = 0, 1
+            call f%init([ptr(x), ptr(y)], [ptr(parse('x^2 + y^2'))], cse=(c == 1))
+            call assert(f%num_inputs() == 2)
+            call assert(f%num_outputs() == 1)
+            do i = 1, 3
+                do j = 1, 3
+                    call f%eval([vals(i), vals(j)], outs)
+                    call assert_close(outs(1), vals(i)**2 + vals(j)**2)
+                end do
+            end do
+        end do
+        print *, "test_point: success"
+    end subroutine
+
+    subroutine test_batch()
+        type(Basic) :: x, y, z
+        type(Lambdify) :: f
+        real(c_double) :: inps(3, 2), outs(2, 2)
+
+        x = Symbol('x')
+        y = Symbol('y')
+        z = Symbol('z')
+        ! r = x + y*z + (y*z)**2, s = 2*x + y*z + (y*z)**2
+        call f%init([ptr(x), ptr(y), ptr(z)], &
+                    [ptr(x + y*z + (y*z)**2), ptr(2*x + y*z + (y*z)**2)], cse=.true.)
+        inps(:, 1) = [1.5_c_double, 2.0_c_double, 3.0_c_double]
+        inps(:, 2) = [1.0_c_double, 1.0_c_double, 1.0_c_double]
+        call f%eval(inps, outs)
+        call assert_close(outs(1, 1), 43.5_c_double)
+        call assert_close(outs(2, 1), 45.0_c_double)
+        call assert_close(outs(1, 2), 3.0_c_double)
+        call assert_close(outs(2, 2), 4.0_c_double)
+        print *, "test_batch: success"
+    end subroutine
+
+    subroutine test_basic_array()
+        type(Basic) :: args(2), exprs(2)
+        type(Lambdify) :: f
+        real(c_double) :: outs(2)
+
+        args(1) = Symbol('x')
+        args(2) = Symbol('y')
+        exprs(1) = sin(args(1)) * args(2)
+        exprs(2) = args(1) + args(2)**2
+        call f%init(args, exprs)
+        call f%eval([1.0_c_double, 2.0_c_double], outs)
+        call assert_close(outs(1), sin(1.0_c_double) * 2.0_c_double)
+        call assert_close(outs(2), 5.0_c_double)
+
+        ! Re-initializing frees the previous function
+        call f%init(args, exprs(2:2))
+        call assert(f%num_outputs() == 1)
+        call f%eval([1.0_c_double, 2.0_c_double], outs(1:1))
+        call assert_close(outs(1), 5.0_c_double)
+        print *, "test_basic_array: success"
+    end subroutine
+
+#ifdef HAVE_SYMENGINE_LLVM
+    subroutine test_llvm()
+        type(Basic) :: x, y
+        type(Lambdify) :: f
+        procedure(lambdify_function), pointer :: fn
+        real(c_double) :: inps(2, 9), outs(1, 9), out(1)
+        real(c_double), parameter :: vals(3) = [-1.0_c_double, 0.0_c_double, 1.0_c_double]
+        integer :: i, j
+
+        x = Symbol('x')
+        y = Symbol('y')
+        call f%init([ptr(x), ptr(y)], [ptr(x**2 + y**2)], backend='llvm')
+        do i = 1, 3
+            do j = 1, 3
+                inps(:, 3 * (i - 1) + j) = [vals(i), vals(j)]
+            end do
+        end do
+        call f%eval(inps, outs)
+        do j = 1, 9
+            call assert_close(outs(1, j), inps(1, j)**2 + inps(2, j)**2)
+        end do
+
+        ! The compiled function itself
+        call f%get_function(fn)
+        do j = 1, 9
+            call fn(inps(:, j), out)
+            call assert_close(out(1), inps(1, j)**2 + inps(2, j)**2)
+        end do
+        print *, "test_llvm: success"
+    end subroutine
+#endif
+
+end program

From 28ad9e59e817abcff7f47e7c2516fb04e55a57fc Mon Sep 17 00:00:00 2001
From: Blaise Bourdin <bourdin@mcmaster.ca>
Date: Mon, 19 Oct 2026 11:51:23 -0400
Subject: [PATCH 3/3] Document how to build with fpm

fpm does not link SymEngine by itself, and passes --link-flag before the
object files, which requires --whole-archive with the GNU linker.
---
 README.md | 29 ++++++++++++++++++++++++++---
 1 file changed, 26 insertions(+), 3 deletions(-)

diff --git a/README.md b/README.md
index 4fd1696..582ecb9 100644
--- a/README.md
+++ b/README.md
@@ -16,6 +16,29 @@ Install SymEngine.f90:
     make
     make install
 
+### Using fpm
+
+fpm does not know how SymEngine was built, so the libraries to link must be
+given with `--link-flag`: `-lsymengine`, followed, for a static SymEngine (the
+default), by the libraries it was built with (e.g. `-lgmp`, `-lLLVM`, `-lflint`,
+`-lmpfr`) and the C++ runtime. For example, for a SymEngine installed in
+`/some/path` and built with GMP and LLVM:
+
+On Linux (fpm passes `--link-flag` before the object files, so the GNU linker
+needs `--whole-archive` to keep SymEngine):
+
+    fpm test --flag "-DHAVE_SYMENGINE_LLVM" --link-flag "-L/some/path/lib \
+        -Wl,--whole-archive -lsymengine -Wl,--no-whole-archive \
+        -L$(llvm-config --libdir) -lgmp -lLLVM -lstdc++"
+
+On macOS, with GMP and LLVM from Homebrew:
+
+    fpm test --flag "-DHAVE_SYMENGINE_LLVM" --link-flag "-L/some/path/lib \
+        -L$(brew --prefix)/lib -L$(brew --prefix llvm)/lib -lsymengine -lgmp -lLLVM -lc++"
+
+`-DHAVE_SYMENGINE_LLVM` enables the LLVM backend of `Lambdify` and should only
+be passed if SymEngine was built with `WITH_LLVM=ON`.
+
 ## Lambdify
 
 `Lambdify` turns expressions into numerical functions, like `Lambdify` in
@@ -39,6 +62,6 @@ Optional arguments of `init`: `cse` (common subexpression elimination),
 `backend` (`'lambda'`, the default, or `'llvm'`) and `opt_level` (LLVM only).
 With `backend='llvm'`, `call f%get_function(fn)` returns the compiled function
 as a procedure pointer, called as `call fn(inps, outs)`. The LLVM backend
-requires SymEngine built with `WITH_LLVM=ON`; when building with fpm, also pass
-`--flag -DHAVE_SYMENGINE_LLVM`. A `Lambdify` object cannot be copied, and with
-the `'lambda'` backend each thread needs its own object.
+requires SymEngine built with `WITH_LLVM=ON` (with fpm, see above). A
+`Lambdify` object cannot be copied, and with the `'lambda'` backend each thread
+needs its own object.