ad_lib.f90.erb

<%
  LOOP_UNROLL_THRESHOLD = 10

  RKINDS = [
    32,
    64,
    128,
  ]
  IKINDS = [
    8,
    16,
    32,
    64,
  ]
  ORDERS = [
    1,
    2,
  ]
  NS = [
    1,
    2,
    5,
    7,
    8,
  ]
  RKINDS_ORDERS_NS = RKINDS.product(ORDERS, NS)
  FN1S = [
    :exp,
    :log,
    :abs,
    :sqrt,
    :sin,
  ]
  OPS = {
    add: :+,
    sub: :-,
    mul: :*,
    div: :/,
    pow: :**,
  }
  COMPS = {
    lt: :<,
    le: :<=,
    eq: :==,
    ne: :"/=",
    ge: :>=,
    gt: :>,
  }
%>
#include "fortran_lib.h"

! todo: optimize when necessary
#define exp_p1(x) (exp(x))
#define exp_p2(x) (exp(x))
#define log_p1(x) (1/(x))
#define log_p2(x) (-1/(x)**2)
#define abs_p1(x) (sign(real(1, kind=kind(x)), x))
#define sqrt_p1(x) (1/(2*sqrt(x)))
#define sqrt_p2(x) (-1/(4*x*sqrt(x)))
#define sin_p1(x) (cos(x))
#define sin_p2(x) (-sin(x))

! assuming i >= j
#define index_hess(i, j, n) (i + n*(j - 1) - (j*(j - 1))/2)

module ad_lib
   <% RKINDS.each{|k|%>
      use, intrinsic:: iso_fortran_env, only: real<%= k %>
   <% } %>
   <% IKINDS.each{|k|%>
      use, intrinsic:: iso_fortran_env, only: int<%= k %>
   <% } %>
   use, intrinsic:: iso_fortran_env, only: int32, int64
   USE_FORTRAN_LIB_H
   use, intrinsic:: iso_fortran_env, only: input_unit, output_unit, error_unit
   use, non_intrinsic:: comparable_lib, only: is_nan

   implicit none

   private
   ! accessors
   public:: kind, real, grad, hess
   public:: epsilon, tiny, huge
   public:: is_nan
   ! functions
   public:: dot_product
   public:: sum
   <% FN1S.each{|fn| %>
      public:: <%= fn %>
   <% } %>
   public:: min, max
   ! operators
   <% OPS.merge(COMPS).each{|_, op| %>
      public:: operator(<%= op %>)
   <% } %>
   public:: assignment(=)

   type, public:: DualNumber
   end type DualNumber

   <% RKINDS_ORDERS_NS.each{|kd, o, n| %>
      type, extends(DualNumber), public:: Dual<%= kd %>_<%= o %>_<%= n %>
         Real(kind=real<%= kd %>):: f = 0
         Real(kind=real<%= kd %>):: g(<%= n %>) = 0
         <% if o > 1 %>
            Real(kind=real<%= kd %>):: h(<%= (n*(n + 1))/2 %>) = 0
         <% end %>
      end type Dual<%= kd %>_<%= o %>_<%= n %>
   <% } %>

   ! Interface
   <% RKINDS_ORDERS_NS.each{|kd, o, n| %>
      interface kind
         module procedure kind_<%= kd %>_<%= o %>_<%= n %>
      end interface kind

      interface real
         module procedure real_<%= kd %>_<%= o %>_<%= n %>
      end interface real

      interface grad
         module procedure grad_<%= kd %>_<%= o %>_<%= n %>
      end interface grad

      <% [:epsilon, :tiny, :huge].each{|fn1| %>
         interface <%= fn1 %>
            module procedure <%= fn1 %>_<%= kd %>_<%= o %>_<%= n %>
         end interface <%= fn1 %>
      <% } %>

      interface is_nan
         module procedure is_nan_<%= kd %>_<%= o %>_<%= n %>
      end interface is_nan


      <% if o > 1 %>
         interface hess
            module procedure hess_<%= kd %>_<%= o %>_<%= n %>
         end interface hess
      <% end %>

      <% RKINDS.each{|kd2| %>
         interface dot_product
            module procedure dot_product_dd_<%= kd %>_<%= kd2 %>_<%= o %>_<%= n %>
            module procedure dot_product_dr_<%= kd %>_<%= kd2 %>_<%= o %>_<%= n %>
            module procedure dot_product_rd_<%= kd2 %>_<%= kd %>_<%= o %>_<%= n %>
         end interface dot_product
      <% } %>

      interface sum
         module procedure sum_<%= kd %>_<%= o %>_<%= n %>_1
      end interface sum

      <% FN1S.each{|fn1| %>
         interface <%= fn1 %>
            module procedure <%= fn1 %>_<%= kd %>_<%= o %>_<%= n %>
         end interface <%= fn1 %>
      <% } %>

      interface operator(-)
         module procedure neg_<%= kd %>_<%= o %>_<%= n %>
      end interface operator(-)

      interface max
         module procedure max_dd_<%= kd %>_<%= o %>_<%= n %>
         module procedure max_dr_<%= kd %>_<%= o %>_<%= n %>
         module procedure max_rd_<%= kd %>_<%= o %>_<%= n %>
      end interface max

      interface min
         module procedure min_dd_<%= kd %>_<%= o %>_<%= n %>
         module procedure min_dr_<%= kd %>_<%= o %>_<%= n %>
         module procedure min_rd_<%= kd %>_<%= o %>_<%= n %>
      end interface min

      <% [:d, :r, :i].each{|t| %>
         <% {d: RKINDS, r: RKINDS, i: IKINDS}.fetch(t).each{|kx| # x ∈ (d, r, i) %>
            <% OPS.merge(COMPS).each{|fn2, op| %>
               interface operator(<%= op %>)
                  module procedure <%= fn2 %>_d<%= t %>_<%= kd %>_<%= kx %>_<%= o %>_<%= n %>
                  <% unless t == :d %>
                     module procedure <%= fn2 %>_<%= t %>d_<%= kx %>_<%= kd %>_<%= o %>_<%= n %>
                  <% end %>
               end interface operator(<%= op %>)
            <% } %>

            <% unless t == :d %>
               interface assignment(=)
                  module procedure assign_d<%= t %>_<%= kd %>_<%= kx %>_<%= o %>_<%= n %>
               end interface assignment(=)
            <% end %>
         <% } %>
      <% } %>
   <% } %>

contains

   <% RKINDS_ORDERS_NS.each{|kd, o, n| %>
      elemental function kind_<%= kd %>_<%= o %>_<%= n %>(x) result(ret)
         type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
         Integer:: ret

         ret = <%= kd %>
      end function kind_<%= kd %>_<%= o %>_<%= n %>

      <% [:epsilon, :tiny, :huge].each{|fn1| %>
         elemental function <%= fn1 %>_<%= kd %>_<%= o %>_<%= n %>(x) result(ret)
            type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
            Real(kind=real<%= kd %>):: ret

            ret = <%= fn1 %>(ret)
         end function <%= fn1 %>_<%= kd %>_<%= o %>_<%= n %>
      <% } %>

      elemental function is_nan_<%= kd %>_<%= o %>_<%= n %>(x) result(ret)
         type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
         Logical:: ret

         ret = is_nan(x%f)
      end function is_nan_<%= kd %>_<%= o %>_<%= n %>

      elemental function real_<%= kd %>_<%= o %>_<%= n %>(x) result(ret)
         type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
         Real(kind=real<%= kd %>):: ret

         ret = x%f
      end function real_<%= kd %>_<%= o %>_<%= n %>

      pure function grad_<%= kd %>_<%= o %>_<%= n %>(x) result(ret)
         type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
         Real(kind=real<%= kd %>):: ret(<%= n %>)

         ret = x%g
      end function grad_<%= kd %>_<%= o %>_<%= n %>

      <% if o > 1 %>
         pure function hess_<%= kd %>_<%= o %>_<%= n %>(x) result(ret)
            type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
            Real(kind=real<%= kd %>):: ret(<%= n %>, <%= n %>)
            Integer:: i, j

            do, concurrent(j = 1:<%= n %>)
               ret(j, j) = x%h(index_hess(j, j, <%= n %>))
               do, concurrent(i = (j + 1):<%= n %>)
                  ret(i, j) = x%h(index_hess(i, j, <%= n %>))
                  ret(j, i) = ret(i, j)
               end do
            end do
         end function hess_<%= kd %>_<%= o %>_<%= n %>
      <% end %>

      <% RKINDS.each{|kd2| %>
         pure function dot_product_dd_<%= kd %>_<%= kd2 %>_<%= o %>_<%= n %>(x, y) result(z)
            type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x(:)
            type(Dual<%= kd2 %>_<%= o %>_<%= n %>), intent(in):: y(size(x))
            type(Dual<%= [kd, kd2].max %>_<%= o %>_<%= n %>):: z

            z = sum(x*y)
         end function dot_product_dd_<%= kd %>_<%= kd2 %>_<%= o %>_<%= n %>
      <% } %>

      <% RKINDS.each{|kr| %>
         <%
            t = :r
            t_decl = "Real(kind=real#{kr})"
            k_ret = [kd, kr].max
            d_r = "dot_product_d#{t}_#{kd}_#{kr}_#{o}_#{n}"
            r_d = "dot_product_#{t}d_#{kr}_#{kd}_#{o}_#{n}"
         %>
         pure function <%= d_r %>(x, y) result(z)
            type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x(:)
            <%= t_decl %>, intent(in):: y(size(x))
            type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

            z = sum(x*y)
         end function <%= d_r %>

         pure function <%= r_d %>(x, y) result(z)
            <%= t_decl %>, intent(in):: x(:)
            type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: y(size(x))
            type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

            z = sum(x*y)
         end function <%= r_d %>
      <% } %>

      pure function sum_<%= kd %>_<%= o %>_<%= n %>_1(x) result(y)
         Integer(kind=int64), parameter:: one = 1
         type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x(:)
         type(Dual<%= kd %>_<%= o %>_<%= n %>):: y

         y = sum_pairwise_<%= kd %>_<%= o %>_<%= n %>_1(x, one, size(x, kind=kind(one)))
      end function sum_<%= kd %>_<%= o %>_<%= n %>_1

      pure recursive function sum_pairwise_<%= kd %>_<%= o %>_<%= n %>_1(x, i1, i2) result(y)
         Integer(kind=int32), parameter:: blocksize = 2**5
         type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x(:)
         Integer(kind=int64), intent(in):: i1, i2
         type(Dual<%= kd %>_<%= o %>_<%= n %>):: y
         Integer(kind=kind(i1)):: i_mid

         if(i1 + blocksize > i2)then
            y = sum_seq_<%= kd %>_<%= o %>_<%= n %>_1(x, i1, i2)
         else
            i_mid = (i1 + i2)/2
            y &
               = sum_pairwise_<%= kd %>_<%= o %>_<%= n %>_1(x, i1, i_mid) &
               + sum_pairwise_<%= kd %>_<%= o %>_<%= n %>_1(x, i_mid + 1, i2)
         end if
      end function sum_pairwise_<%= kd %>_<%= o %>_<%= n %>_1

      pure function sum_seq_<%= kd %>_<%= o %>_<%= n %>_1(x, i1, i2) result(y)
         type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x(:)
         Integer(kind=int64), intent(in):: i1, i2
         type(Dual<%= kd %>_<%= o %>_<%= n %>):: y
         Integer(kind=kind(i1)):: i

         do i = i1, i2
            y = y + x(i)
         end do
      end function sum_seq_<%= kd %>_<%= o %>_<%= n %>_1

      <% FN1S.each{|fn1| %>
         elemental function <%= fn1 %>_<%= kd %>_<%= o %>_<%= n %>(x) result(y)
            type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
            type(Dual<%= kd %>_<%= o %>_<%= n %>):: y
            Real(kind=real<%= kd %>):: xf, fn_p1
            <% if o > 1 %>
               <% unless fn1 == :abs %>
                  <% if n > LOOP_UNROLL_THRESHOLD %>
                     Integer:: i, j
                  <% end %>
                  Real(kind=real<%= kd %>):: fn_p2, fn_p2_xj
               <% end %>
            <% end %>

            xf = x%f
            y%f = <%= fn1 %>(xf)
            fn_p1 = <%= fn1 %>_p1(xf)
            y%g = fn_p1*x%g
            <% if o > 1 %>
               y%h = fn_p1*x%h
               <% unless fn1 == :abs %>
                  fn_p2 = <%= fn1 %>_p2(xf)
                  <% if n > LOOP_UNROLL_THRESHOLD %>
                     do, concurrent(j = 1:<%= n %>)
                        fn_p2_xj = fn_p2*x%g(j)
                        do, concurrent(i = j:<%= n %>)
                           y%h(index_hess(i, j, <%= n %>)) = y%h(index_hess(i, j, <%= n %>)) + fn_p2_xj*x%g(i)
                        end do
                     end do
                  <% else %>
                     <% (1..n).each{|j| %>
                        fn_p2_xj = fn_p2*x%g(<%= j %>)
                        <% (j..n).each{|i| %>
                           y%h(index_hess(<%= i %>, <%= j %>, <%= n %>)) = y%h(index_hess(<%= i %>, <%= j %>, <%= n %>)) + fn_p2_xj*x%g(<%= i %>)
                        <% } %>
                     <% } %>
                  <% end %>
               <% end %>
            <% end %>
         end function <%= fn1 %>_<%= kd %>_<%= o %>_<%= n %>
      <% } %>

      elemental function neg_<%= kd %>_<%= o %>_<%= n %>(x) result(y)
         type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
         type(Dual<%= kd %>_<%= o %>_<%= n %>):: y

         y%f = -x%f
         y%g = -x%g
         <% if o > 1 %>
            y%h = -x%h
         <% end %>
      end function neg_<%= kd %>_<%= o %>_<%= n %>

      <% RKINDS.each{|kd2| %>
         <% suffix = "_dd_#{kd}_#{kd2}_#{o}_#{n}"%>

         <% [:add, :sub].each{|fn2| %>
            <% op = OPS[fn2] %>
            elemental function <%= fn2 %><%= suffix %>(x, y) result(z)
               type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
               type(Dual<%= kd2 %>_<%= o %>_<%= n %>), intent(in):: y
               type(Dual<%= [kd, kd2].max %>_<%= o %>_<%= n %>):: z

               z%f = x%f <%= op %> y%f
               z%g = x%g <%= op %> y%g
               <% if o > 1 %>
                  z%h = x%h <%= op %> y%h
               <% end %>
            end function <%= fn2 %><%= suffix %>
         <% } %>

         elemental function mul_dd_<%= kd %>_<%= kd2 %>_<%= o %>_<%= n %>(x, y) result(z)
            type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
            type(Dual<%= kd2 %>_<%= o %>_<%= n %>), intent(in):: y
            type(Dual<%= [kd, kd2].max %>_<%= o %>_<%= n %>):: z
            Real(kind=real<%= [kd, kd2].max %>):: xf, yf
            <% if o > 1 %>
               <% if n > LOOP_UNROLL_THRESHOLD %>
                  Integer:: i, j
               <% end %>
               Real(kind=real<%= [kd, kd2].max %>):: xj, yj
            <% end %>

            xf = x%f
            yf = y%f
            z%f = xf*yf

            z%g = x%g*yf + y%g*xf
            <% if o > 1 %>
               <% if n > LOOP_UNROLL_THRESHOLD %>
                  do, concurrent(j = 1:<%= n %>)
                     xj = x%g(j)
                     yj = y%g(j)
                     do, concurrent(i = j:<%= n %>)
                        z%h(index_hess(i, j, <%= n %>)) = x%h(index_hess(i, j, <%= n %>))*yf + y%h(index_hess(i, j, <%= n %>))*xf + x%g(i)*yj + y%g(i)*xj
                     end do
                  end do
               <% else %>
                  <% (1..n).each{|j| %>
                     xj = x%g(<%= j %>)
                     yj = y%g(<%= j %>)
                     <% (j..n).each{|i| %>
                        z%h(index_hess(<%= i %>, <%= j %>, <%= n %>)) = x%h(index_hess(<%= i %>, <%= j %>, <%= n %>))*yf + y%h(index_hess(<%= i %>, <%= j %>, <%= n %>))*xf + x%g(<%= i %>)*yj + y%g(<%= i %>)*xj
                     <% } %>
                  <% } %>
               <% end %>
            <% end %>
         end function mul_dd_<%= kd %>_<%= kd2 %>_<%= o %>_<%= n %>

         elemental function div_dd_<%= kd %>_<%= kd2 %>_<%= o %>_<%= n %>(x, y) result(z)
            type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
            type(Dual<%= kd2 %>_<%= o %>_<%= n %>), intent(in):: y
            type(Dual<%= [kd, kd2].max %>_<%= o %>_<%= n %>):: z
            Real(kind=real<%= [kd, kd2].max %>):: xf, yf, div_y, div_y2, x_div_y2, fn_py
            <% if o > 1 %>
               <% if n > LOOP_UNROLL_THRESHOLD %>
                  Integer:: i, j
               <% end %>
               Real(kind=real<%= [kd, kd2].max %>):: xj, yj, fn_pxy, fn_pyy, fn_pxx_xj_add_fn_pxy_yj, fn_pyx_xj_add_fn_pyy_yj
            <% end %>

            xf = x%f
            yf = y%f
            div_y = 1/yf
            z%f = xf*div_y

            ! fn_px = div_y
            div_y2 = div_y*div_y
            x_div_y2 = xf*div_y2
            fn_py = -x_div_y2
            z%g = x%g*div_y + y%g*fn_py
            <% if o > 1 %>
               ! fn_pxx = 0
               fn_pxy = -div_y2
               fn_pyy = 2*x_div_y2*div_y
               <% if n > LOOP_UNROLL_THRESHOLD %>
                  do, concurrent(j = 1:<%= n %>)
                     xj = x%g(j)
                     yj = y%g(j)
                     fn_pxx_xj_add_fn_pxy_yj = fn_pxy*yj ! + fn_pxx*xj
                     fn_pyx_xj_add_fn_pyy_yj = fn_pxy*xj + fn_pyy*yj
                     do, concurrent(i = j:<%= n %>)
                        z%h(index_hess(i, j, <%= n %>)) = x%h(index_hess(i, j, <%= n %>))*div_y + y%h(index_hess(i, j, <%= n %>))*fn_py + x%g(i)*fn_pxx_xj_add_fn_pxy_yj + y%g(i)*fn_pyx_xj_add_fn_pyy_yj
                     end do
                  end do
               <% else %>
                  <% (1..n).each{|j| %>
                     xj = x%g(<%= j %>)
                     yj = y%g(<%= j %>)
                     fn_pxx_xj_add_fn_pxy_yj = fn_pxy*yj ! + fn_pxx*xj
                     fn_pyx_xj_add_fn_pyy_yj = fn_pxy*xj + fn_pyy*yj
                     <% (j..n).each{|i| %>
                        z%h(index_hess(<%= i %>, <%= j %>, <%= n %>)) = x%h(index_hess(<%= i %>, <%= j %>, <%= n %>))*div_y + y%h(index_hess(<%= i %>, <%= j %>, <%= n %>))*fn_py + x%g(<%= i %>)*fn_pxx_xj_add_fn_pxy_yj + y%g(<%= i %>)*fn_pyx_xj_add_fn_pyy_yj
                     <% } %>
                  <% } %>
               <% end %>
            <% end %>
         end function div_dd_<%= kd %>_<%= kd2 %>_<%= o %>_<%= n %>

         elemental function pow_dd_<%= kd %>_<%= kd2 %>_<%= o %>_<%= n %>(x, y) result(z)
            type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
            type(Dual<%= kd2 %>_<%= o %>_<%= n %>), intent(in):: y
            type(Dual<%= [kd, kd2].max %>_<%= o %>_<%= n %>):: z
            Real(kind=real<%= [kd, kd2].max %>):: xf, yf, fn_px, fn_py, x_pow_y, x_pow_y_sub_1, log_x
            <% if o > 1 %>
               <% if n > LOOP_UNROLL_THRESHOLD %>
                  Integer:: i, j
               <% end %>
               Real(kind=real<%= [kd, kd2].max %>):: xj, yj, fn_pxx, fn_pxy, fn_pyy, fn_pxx_xj_add_fn_pxy_yj, fn_pyx_xj_add_fn_pyy_yj
            <% end %>

            xf = x%f
            yf = y%f
            x_pow_y = xf**yf
            z%f = x_pow_y
            ! this is not mathematically correct but might be useful enough in practice.
            ! todo: `0**0`
            if(xf == 0) return

            x_pow_y_sub_1 = xf**(yf - 1)
            fn_px = yf*x_pow_y_sub_1
            log_x = log(xf)
            fn_py = log_x*x_pow_y
            z%g = x%g*fn_px + y%g*fn_py
            <% if o > 1 %>
               fn_pxx = yf*(yf - 1)*xf**(yf - 2)
               fn_pxy = (1 + yf*log_x)*x_pow_y_sub_1
               fn_pyy = log_x*fn_py ! = log_x**2*x_pow_y
               <% if n > LOOP_UNROLL_THRESHOLD %>
                  do, concurrent(j = 1:<%= n %>)
                     xj = x%g(j)
                     yj = y%g(j)
                     fn_pxx_xj_add_fn_pxy_yj = fn_pxx*xj + fn_pxy*yj
                     fn_pyx_xj_add_fn_pyy_yj = fn_pxy*xj + fn_pyy*yj
                     do, concurrent(i = j:<%= n %>)
                        z%h(index_hess(i, j, <%= n %>)) = x%h(index_hess(i, j, <%= n %>))*fn_px + y%h(index_hess(i, j, <%= n %>))*fn_py + x%g(i)*fn_pxx_xj_add_fn_pxy_yj + y%g(i)*fn_pyx_xj_add_fn_pyy_yj
                     end do
                  end do
               <% else %>
                  <% (1..n).each{|j| %>
                     xj = x%g(<%= j %>)
                     yj = y%g(<%= j %>)
                     fn_pxx_xj_add_fn_pxy_yj = fn_pxx*xj + fn_pxy*yj
                     fn_pyx_xj_add_fn_pyy_yj = fn_pxy*xj + fn_pyy*yj
                     <% (j..n).each{|i| %>
                        z%h(index_hess(<%= i %>, <%= j %>, <%= n %>)) = x%h(index_hess(<%= i %>, <%= j %>, <%= n %>))*fn_px + y%h(index_hess(<%= i %>, <%= j %>, <%= n %>))*fn_py + x%g(<%= i %>)*fn_pxx_xj_add_fn_pxy_yj + y%g(<%= i %>)*fn_pyx_xj_add_fn_pyy_yj
                     <% } %>
                  <% } %>
               <% end %>
            <% end %>
         end function pow_dd_<%= kd %>_<%= kd2 %>_<%= o %>_<%= n %>
      <% } %>

      <% [:r, :i].each{|t| %>
         <% (t == :r ? RKINDS : IKINDS).each{|kri| %>
            <%
              t_decl = t == :r ? "Real(kind=real#{kri})" : "Integer(kind=int#{kri})"
              # this promotion rule emulates Fortran
              k_ret = t == :r ? [kd, kri].max : kd
              d_ri = "_d#{t}_#{kd}_#{kri}_#{o}_#{n}"
              ri_d = "_#{t}d_#{kri}_#{kd}_#{o}_#{n}"
            %>
            elemental function add<%= d_ri %>(x, y) result(z)
               type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
               <%= t_decl %>, intent(in):: y
               type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

               z%f = x%f + y
               z%g = x%g
               <% if o > 1 %>
                  z%h = x%h
               <% end %>
            end function add<%= d_ri %>

            elemental function add<%= ri_d %>(x, y) result(z)
               <%= t_decl %>, intent(in):: x
               type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: y
               type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

               z%f = x + y%f
               z%g = y%g
               <% if o > 1 %>
                  z%h = y%h
               <% end %>
            end function add<%= ri_d %>

            elemental function sub<%= d_ri %>(x, y) result(z)
               type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
               <%= t_decl %>, intent(in):: y
               type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

               z%f = x%f - y
               z%g = x%g
               <% if o > 1 %>
                  z%h = x%h
               <% end %>
            end function sub<%= d_ri %>

            elemental function sub<%= ri_d %>(x, y) result(z)
               <%= t_decl %>, intent(in):: x
               type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: y
               type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

               z%f = x - y%f
               z%g = -y%g
               <% if o > 1 %>
                  z%h = -y%h
               <% end %>
            end function sub<%= ri_d %>

            elemental function mul<%= d_ri %>(x, y) result(z)
               type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
               <%= t_decl %>, intent(in):: y
               type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

               z%f = x%f*y
               z%g = x%g*y
               <% if o > 1 %>
                  z%h = x%h*y
               <% end %>
            end function mul<%= d_ri %>

            elemental function mul<%= ri_d %>(x, y) result(z)
               <%= t_decl %>, intent(in):: x
               type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: y
               type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

               z%f = x*y%f
               z%g = x*y%g
               <% if o > 1 %>
                  z%h = x*y%h
               <% end %>
            end function mul<%= ri_d %>

            elemental function div<%= d_ri %>(x, y) result(z)
               type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
               <%= t_decl %>, intent(in):: y
               type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

               z%f = x%f/y
               z%g = x%g/y
               <% if o > 1 %>
                  z%h = x%h/y
               <% end %>
            end function div<%= d_ri %>

            elemental function div<%= ri_d %>(x, y) result(z)
               <%= t_decl %>, intent(in):: x
               type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: y
               type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

               ! todo: optimize
               z = Dual<%= k_ret %>_<%= o %>_<%= n %>(x)/y
            end function div<%= ri_d %>

            elemental function pow<%= d_ri %>(x, y) result(z)
               type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
               <%= t_decl %>, intent(in):: y
               type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

               ! todo: optimize
               z = x**Dual<%= k_ret %>_<%= o %>_<%= n %>(y)
            end function pow<%= d_ri %>

            elemental function pow<%= ri_d %>(x, y) result(z)
               <%= t_decl %>, intent(in):: x
               type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: y
               type(Dual<%= k_ret %>_<%= o %>_<%= n %>):: z

               ! todo: optimize
               z = Dual<%= k_ret %>_<%= o %>_<%= n %>(x)**y
            end function pow<%= ri_d %>
         <% } %>
      <% } %>

      <% [[:min, :>], [:max, :<]].each{|fn, comp| %>
         <% [[:d, :d], [:d, :r], [:r, :d]].each{|t1, t2| %>
            <%
               decls = {
                  d: "type(Dual#{kd}_#{o}_#{n})",
                  r: "Real(kind=real#{kd})",
               }
            %>
            elemental function <%= fn %>_<%= t1 %><%= t2 %>_<%= kd %>_<%= o %>_<%= n %>(x, y) result(ret)
               <%= decls.fetch(t1) %>, intent(in):: x
               <%= decls.fetch(t2) %>, intent(in):: y
               type(Dual<%= kd %>_<%= o %>_<%= n %>):: ret

               if(is_nan(x))then
                  ret = y
               else if(is_nan(y))then
                  ret = x
               else if(x <%= comp %> y)then
                  ret = y
               else
                  ret = x
               end if
            end function <%= fn %>_<%= t1 %><%= t2 %>_<%= kd %>_<%= o %>_<%= n %>
         <% } %>
      <% } %>

      <% [:d, :r, :i].each{|t| %>
         <% {d: RKINDS, r: RKINDS, i: IKINDS}.fetch(t).each{|kx| # x ∈ (d, r, i) %>
            <%
               t_decl = {
                  d: "type(Dual#{kx}_#{o}_#{n})",
                  r: "Real(kind=real#{kx})",
                  i: "Integer(kind=int#{kx})",
               }.fetch(t)
               y = if t == :d
                  "y%f"
               else
                  "y"
               end
            %>
            <% COMPS.each{|fn2, op| %>
               elemental function <%= fn2 %>_d<%= t %>_<%= kd %>_<%= kx %>_<%= o %>_<%= n %>(x, y) result(ret)
                  type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: x
                  <%= t_decl %>, intent(in):: y
                  Logical:: ret

                  ret = x%f <%= op %> <%= y %>
               end function <%= fn2 %>_d<%= t %>_<%= kd %>_<%= kx %>_<%= o %>_<%= n %>

               <% unless t == :d %>
                  elemental function <%= fn2 %>_<%= t %>d_<%= kx %>_<%= kd %>_<%= o %>_<%= n %>(x, y) result(ret)
                     <%= t_decl %>, intent(in):: x
                     type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(in):: y
                     Logical:: ret

                     ret = x <%= op %> y%f
                  end function <%= fn2 %>_<%= t %>d_<%= kx %>_<%= kd %>_<%= o %>_<%= n %>
               <% end %>
            <% } %>

            <% unless t == :d %>
               elemental subroutine assign_d<%= t %>_<%= kd %>_<%= kx %>_<%= o %>_<%= n %>(to, from)
                  type(Dual<%= kd %>_<%= o %>_<%= n %>), intent(out):: to
                  <%= t_decl %>, intent(in):: from

                  ! to%g and to%h are automatically initialized to 0
                  to%f = from
               end subroutine assign_d<%= t %>_<%= kd %>_<%= kx %>_<%= o %>_<%= n %>
            <% end %>
         <% } %>
      <% } %>
   <% } %>
end module ad_lib