module mod_natural_number
  implicit none

  type natural_number
    private
    integer(1), pointer :: nn(:) => null()
  contains
    final :: finalize
  end type natural_number

  interface finalize
    module procedure finalize
  end interface finalize

  interface to_string
    module procedure to_string
  end interface to_string

  interface show
    module procedure show
  end interface show

  interface size
    module procedure sizen
  end interface size

  interface lshift
    module procedure lshiftn
  end interface lshift

  interface rshift
    module procedure rshiftn
  end interface rshift

  interface natural_number
    module procedure new, new0, new1, new2
  end interface natural_number

  interface assignment(=)
    module procedure set, set0, set1, set2
  end interface assignment(=)

  interface cmp
    module procedure cmp, cmp0ns, cmp0sn, cmp1nl, cmp1ln, cmp2ni, cmp2in
  end interface cmp

  interface operator(+)
    module procedure add, add0ns, add0sn, add1nl, add1ln, add2ni, add2in
  end interface operator(+)

  interface operator(-)
    module procedure sub, sub0ns, sub0sn, sub1nl, sub1ln, sub2ni, sub2in
  end interface operator(-)

  interface operator(*)
    module procedure mul, mul0ns, mul0sn, mul1nl, mul1ln, mul2ni, mul2in
  end interface operator(*)

  interface operator(/)
    module procedure div, div0ns, div0sn, div1nl, div1ln, div2ni, div2in
  end interface operator(/)

  interface mod
    module procedure rem, rem0ns, rem0sn, rem1nl, rem1ln, rem2ni, rem2in
  end interface mod

  interface operator(**)
    module procedure pow, pow0ns, pow0sn, pow1nl, pow1ln, pow2ni, pow2in
  end interface operator(**)

  interface operator(>)
    module procedure gtt, gtt0ns, gtt0sn, gtt1nl, gtt1ln, gtt2ni, gtt2in
  end interface operator(>)

  interface operator(<)
    module procedure lst, lst0ns, lst0sn, lst1nl, lst1ln, lst2ni, lst2in
  end interface operator(<)

  interface operator(>=)
    module procedure goe, goe0ns, goe0sn, goe1nl, goe1ln, goe2ni, goe2in
  end interface operator(>=)

  interface operator(<=)
    module procedure loe, loe0ns, loe0sn, loe1nl, loe1ln, loe2ni, loe2in
  end interface operator(<=)

  interface operator(==)
    module procedure eql, eql0ns, eql0sn, eql1nl, eql1ln, eql2ni, eql2in
  end interface operator(==)

  interface operator(/=)
    module procedure neq, neq0ns, neq0sn, neq1nl, neq1ln, neq2ni, neq2in
  end interface operator(/=)

  interface gcd
    module procedure gcd, gcd0ns, gcd0sn, gcd1nl, gcd1ln, gcd2ni, gcd2in
  end interface gcd

  interface lcm
    module procedure lcm, lcm0ns, lcm0sn, lcm1nl, lcm1ln, lcm2ni, lcm2in
  end interface lcm

  interface and
    module procedure andn, and0ns, and0sn, and1nl, and1ln, and2ni, and2in
  end interface and

  interface or
    module procedure orn, or0ns, or0sn, or1nl, or1ln, or2ni, or2in
  end interface or

  interface xor
    module procedure xorn, xor0ns, xor0sn, xor1nl, xor1ln, xor2ni, xor2in
  end interface xor

  interface max
    module procedure maxn, max0ns, max0sn, max1nl, max1ln, max2ni, max2in
  end interface max

  interface min
    module procedure minn, min0ns, min0sn, min1nl, min1ln, min2ni, min2in
  end interface min

  interface dot_product
    module procedure dot_productn
  end interface dot_product

  interface mat_mul
    module procedure mat_mulnmv, mat_mulnvm, mat_mulnmm
  end interface mat_mul

contains

  subroutine finalize(n)
    type(natural_number), intent(inout) :: n
    if (associated(n%nn)) deallocate(n%nn)
  end

  function to_string(n) result(s)
    type(natural_number), intent(in) :: n
    character(sizen(n)) :: s
    integer :: i, u
    if (.not.associated(n%nn)) then
      s = "0"
      return
    end if
    u = sizen(n)
    do i = 1, u
      s(i:i) = achar(n%nn(u-i)+48)
    end do
  end

  subroutine show(n)
    type(natural_number), intent(in) :: n
    write(*,'(a)') to_string(n)
  end

  pure integer function sizen(n)
    type(natural_number), intent(in) :: n
    sizen = 1
    if (associated(n%nn)) sizen = ubound(n%nn,1)+1
  end

  logical function is_zero(n)
    type(natural_number), intent(in) :: n
    is_zero = .not.associated(n%nn)
  end

  subroutine clear(n)
    type(natural_number), intent(inout) :: n
    if (associated(n%nn)) deallocate(n%nn)
  end

  subroutine norm(n)
    type(natural_number), intent(inout) :: n
    integer(1), allocatable :: tmp(:)
    integer :: i
    if (.not.associated(n%nn)) return
    i = ubound(n%nn,1)
    if (n%nn(i) /= 0) return
    do while (i >= 0 .and. n%nn(i) == 0)
      i = i-1
    end do
    if (i < 0) then
      deallocate(n%nn)
      return
    end if
    allocate(tmp(0:i))
    tmp(0:i) = n%nn(0:i)
    deallocate(n%nn)
    allocate(n%nn(0:i))
    n%nn(0:i) = tmp(0:i)
    deallocate(tmp)
  end

  function lshiftn(x,k) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: k
    type(natural_number) :: n
    integer :: u
    if (.not.associated(x%nn)) return
    u = ubound(x%nn,1)
    allocate(n%nn(0:u+k))
    n%nn = 0
    n%nn(k:u+k) = x%nn(0:u)
  end

  function rshiftn(x,k) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: k
    type(natural_number) :: n
    integer :: u
    if (.not.associated(x%nn)) return
    u = ubound(x%nn,1)
    if (k > u) return
    allocate(n%nn(0:u-k))
    n%nn = 0
    n%nn(0:u-k) = x%nn(k:u)
  end

  function new() result(n)
    type(natural_number) :: n
    allocate(n%nn(0:0))
    n%nn(0) = 1
  end

  subroutine set(n,x)
    type(natural_number), intent(out) :: n
    type(natural_number), intent(in) :: x
    call clear(n)
    if (.not.associated(x%nn)) return
    allocate(n%nn(0:ubound(x%nn,1)))
    n%nn = x%nn
  end

  integer function cmp(x,y)
    type(natural_number), intent(in) :: x, y
    integer :: i, ux, uy
    cmp = 0
    ux = -1
    uy = -1
    if (associated(x%nn)) ux = ubound(x%nn,1)
    if (associated(y%nn)) uy = ubound(y%nn,1)
    cmp = ux-uy
    if (cmp /= 0) return
    do i = ux, 0, -1
      cmp = x%nn(i)-y%nn(i)
      if (cmp /= 0) return
    end do
  end

  function add(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n
    integer :: u, i
    if (.not.associated(y%nn)) then
      call set(n,x)
      return
    end if
    if (.not.associated(x%nn)) then
      call set(n,y)
      return
    end if
    u = max(ubound(x%nn,1),ubound(y%nn,1))
    allocate(n%nn(0:u+1))
    n%nn = 0
    do i = 0, ubound(x%nn,1)
      n%nn(i) = x%nn(i)
    end do
    do i = 0, ubound(y%nn,1)
      n%nn(i) = n%nn(i)+y%nn(i)
    end do
    do i = 0, u
      if (n%nn(i) > 9) then
        n%nn(i) = n%nn(i)-10
        n%nn(i+1) = n%nn(i+1)+1
      end if
    end do
    call norm(n)
  end

  function sub(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n
    integer :: u, i
    if (.not.associated(y%nn)) then
      call set(n,x)
      return
    end if
    if (.not.associated(x%nn) .or. ubound(x%nn,1) < ubound(y%nn,1)) then
      write(*,'(a)') "Error: Overflow (natural_number, sub)"
      stop
    end if
    u = ubound(x%nn,1)
    allocate(n%nn(0:u))
    n%nn = x%nn
    do i = 0, ubound(y%nn,1)
      n%nn(i) = n%nn(i)-y%nn(i)
    end do
    do i = 0, u-1
      if (n%nn(i) < 0) then
        n%nn(i) = n%nn(i)+10
        n%nn(i+1) = n%nn(i+1)-1
      end if
    end do
    if (n%nn(u) < 0) then
      write(*,'(a)') "Error: Overflow (natural_number, sub)"
      stop
    end if
    call norm(n)
  end

  function mul(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n
    n = karatsuba_mul(x,y)
  end

  function primitive_mul(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n
    integer :: ux, uy, u, i, j
    if (.not.(associated(x%nn).and.associated(y%nn))) return
    ux = ubound(x%nn,1)
    uy = ubound(y%nn,1)
    u = ux+uy+1
    allocate(n%nn(0:u))
    n%nn = 0
    do i = 0, ux
      do j = 0, uy
        n%nn(i+j) = n%nn(i+j)+x%nn(i)*y%nn(j)
        n%nn(i+j+1) = n%nn(i+j+1)+n%nn(i+j)/10
        n%nn(i+j) = mod(n%nn(i+j),10)
      end do
    end do
    call norm(n)
  end

  recursive function karatsuba_mul(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n, z0, z1, z2
    type(natural_number) :: x0, x1, y0, y1
    integer :: u, ux, uy
    if (.not.(associated(x%nn).and.associated(y%nn))) return
    ux = ubound(x%nn,1)
    uy = ubound(y%nn,1)
    u = max(ux,uy)
    if (u < 1024) then
      n = primitive_mul(x,y)
      return
    end if
    u = (u+1)/2-1
    allocate(x0%nn(0:min(ux,u)))
    x0%nn(0:min(ux,u)) = x%nn(0:min(ux,u))
    call norm(x0)
    if (ux >= u+1) then
      allocate(x1%nn(0:ux-u-1))
      x1%nn(0:ux-u-1) = x%nn(u+1:ux)
    end if
    allocate(y0%nn(0:min(uy,u)))
    y0%nn(0:min(uy,u)) = y%nn(0:min(uy,u))
    call norm(y0)
    if (uy >= u+1) then
      allocate(y1%nn(0:uy-u-1))
      y1%nn(0:uy-u-1) = y%nn(u+1:uy)
    end if
    z0 = karatsuba_mul(x0,y0)
    if (associated(x1%nn) .and. associated(y1%nn)) then
      z2 = karatsuba_mul(x1,y1)
    end if
    z1 = add(z0,z2)
    if (neq(x0,x1) .and. neq(y0,y1)) then
      if (gtt(x1,x0) .and. gtt(y1,y0)) then
        z1 = sub(z1,karatsuba_mul(sub(x1,x0),sub(y1,y0)))
      else if (gtt(x1,x0)) then
        z1 = add(z1,karatsuba_mul(sub(x1,x0),sub(y0,y1)))
      else if (gtt(y1,y0)) then
        z1 = add(z1,karatsuba_mul(sub(x0,x1),sub(y1,y0)))
      else
        z1 = sub(z1,karatsuba_mul(sub(x0,x1),sub(y0,y1)))
      end if
    end if
    z2 = lshiftn(z2,2*(u+1))
    z1 = lshiftn(z1,u+1)
    n = add(add(z0,z1),z2)
  end

  function div(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n, r
    call effective_div(x,y,n,r)
  end

  function rem(x,y) result(r)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n, r
    call effective_div(x,y,n,r)
  end

  subroutine naive_div(x,y,n,r)
    type(natural_number), intent(in) :: x, y
    type(natural_number), intent(out) :: n, r
    integer :: ux, uy, i, j
    if (.not.associated(y%nn)) then
      write(*,'(a)') "Error: Division by zero (natural_number, naive_div)"
      stop
    end if
    call clear(n)
    call set(r,x)
    if (.not.associated(x%nn)) return
    ux = ubound(x%nn,1)
    uy = ubound(y%nn,1)
    if (ux < uy) return
    allocate(n%nn(0:ux-uy))
    n%nn = 0
    do i = ux-uy, 0, -1
      do
        do j = 0, uy
          r%nn(i+j) = r%nn(i+j)-y%nn(j)
          if (j /= uy .and. r%nn(i+j) < 0) then
            r%nn(i+j) = r%nn(i+j)+10
            r%nn(i+j+1) = r%nn(i+j+1)-1
          end if
        end do
        if (r%nn(i+uy) < 0) then
          do j = 0, uy
            r%nn(i+j) = r%nn(i+j)+y%nn(j)
            if (j /= uy .and. r%nn(i+j) > 9) then
              r%nn(i+j) = r%nn(i+j)-10
              r%nn(i+j+1) = r%nn(i+j+1)+1
            end if
          end do
          if (i+uy > 0) then
            r%nn(i+uy-1) = r%nn(i+uy-1)+10*r%nn(i+uy)
            r%nn(i+uy) = 0
          end if
          exit
        end if
        n%nn(i) = n%nn(i)+1
      end do
    end do
    if (uy > 0 .and. r%nn(uy-1) > 9) then
      r%nn(uy) = r%nn(uy-1)/10
      r%nn(uy-1) = mod(r%nn(uy-1),10)
    end if
    call norm(n)
    call norm(r)
  end

  subroutine primitive_div(x,y,n,r)
    type(natural_number), intent(in) :: x, y
    type(natural_number), intent(out) :: n, r
    type(natural_number) :: ky
    integer :: ux, uy, i, j
    integer(1) :: k
    if (.not.associated(y%nn)) then
      write(*,'(a)') "Error: Division by zero (natural_number, primitive_div)"
      stop
    end if
    call clear(n)
    call set(r,x)
    if (.not.associated(x%nn)) return
    ux = ubound(x%nn,1)
    uy = ubound(y%nn,1)
    if (ux < uy) return
    allocate(n%nn(0:ux-uy))
    n%nn = 0
    allocate(ky%nn(0:uy))
    do i = ux-uy, 0, -1
      k = min(r%nn(i+uy)/y%nn(uy),9)
      ky%nn(0:uy) = k*y%nn(0:uy)
      do j = 0, uy-1
        ky%nn(j+1) = ky%nn(j+1)+ky%nn(j)/10
        ky%nn(j) = mod(ky%nn(j),10)
      end do
      r%nn(i:i+uy) = r%nn(i:i+uy)-ky%nn(0:uy)
      do j = i, i+uy-1
        if (r%nn(j) < 0) then
          r%nn(j) = r%nn(j)+10
          r%nn(j+1) = r%nn(j+1)-1
        end if
      end do
      do while (r%nn(i+uy) < 0)
        r%nn(i:i+uy) = r%nn(i:i+uy)+y%nn(0:uy)
        do j = i, i+uy-1
          if (r%nn(j) > 9) then
            r%nn(j) = r%nn(j)-10
            r%nn(j+1) = r%nn(j+1)+1
          end if
        end do
        k = k-1
      end do
      if (i+uy > 0) then
        r%nn(i+uy-1) = r%nn(i+uy-1)+10*r%nn(i+uy)
        r%nn(i+uy) = 0
      end if
      n%nn(i) = k
    end do
    if (uy > 0 .and. r%nn(uy-1) > 9) then
      r%nn(uy) = r%nn(uy-1)/10
      r%nn(uy-1) = mod(r%nn(uy-1),10)
    end if
    call norm(n)
    call norm(r)
  end

  subroutine effective_div(x,y,n,r)
    type(natural_number), intent(in) :: x, y
    type(natural_number), intent(out) :: n, r
    type(natural_number) :: ky
    integer :: ux, uy, i, j, m, tmp
    integer(1) :: k
    if (.not.associated(y%nn)) then
      write(*,'(a)') "Error: Division by zero (natural_number, effective_div)"
      stop
    end if
    call clear(n)
    call set(r,x)
    if (.not.associated(x%nn)) return
    ux = ubound(x%nn,1)
    uy = ubound(y%nn,1)
    if (ux < uy) return
    allocate(n%nn(0:ux-uy))
    n%nn = 0
    allocate(ky%nn(0:uy))
    do i = ux-uy, 0, -1
      if (uy == 0) then
        k = min(r%nn(i+uy)/y%nn(uy),9)
      else
        tmp = int(r%nn(i+uy),4)*10+int(r%nn(i+uy-1),4)
        tmp = tmp/int(10*y%nn(uy)+y%nn(uy-1),4)
        k = min(int(tmp,1),9)
      end if
      ky%nn(0:uy) = k*y%nn(0:uy)
      do j = 0, uy-1
        ky%nn(j+1) = ky%nn(j+1)+ky%nn(j)/10
        ky%nn(j) = mod(ky%nn(j),10)
      end do
      r%nn(i:i+uy) = r%nn(i:i+uy)-ky%nn(0:uy)
      do j = i, i+uy-1
        if (r%nn(j) < 0) then
          r%nn(j) = r%nn(j)+10
          r%nn(j+1) = r%nn(j+1)-1
        end if
      end do
      do while (r%nn(i+uy) < 0)
        r%nn(i:i+uy) = r%nn(i:i+uy)+y%nn(0:uy)
        do j = i, i+uy-1
          if (r%nn(j) > 9) then
            r%nn(j) = r%nn(j)-10
            r%nn(j+1) = r%nn(j+1)+1
          end if
        end do
        k = k-1
      end do
      if (i+uy > 0) then
        r%nn(i+uy-1) = r%nn(i+uy-1)+10*r%nn(i+uy)
        r%nn(i+uy) = 0
      end if
      n%nn(i) = k
    end do
    if (uy > 0 .and. r%nn(uy-1) > 9) then
      r%nn(uy) = r%nn(uy-1)/10
      r%nn(uy-1) = mod(r%nn(uy-1),10)
    end if
    call norm(n)
    call norm(r)
  end

  function half(x) result(n)
    type(natural_number), intent(in) :: x
    type(natural_number) :: n
    integer :: i
    if (.not.associated(x%nn)) return
    call set(n,x)
    do i = ubound(n%nn,1), 0, -1
      if (i /= 0) n%nn(i-1) = n%nn(i-1)+10*mod(n%nn(i),2)
      n%nn(i) = n%nn(i)/2
    end do
    call norm(n)
  end

  function pow(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n, p, m
    if (.not.associated(x%nn)) return
    n = new()
    call set(p,x)
    call set(m,y)
    do while (associated(m%nn))
      if (mod(m%nn(0),2) == 1) n = mul(n,p)
      p = mul(p,p)
      m = half(m)
    end do
    call norm(n)
  end

  logical function gtt(x,y)
    type(natural_number), intent(in) :: x, y
    gtt = cmp(x,y) > 0
  end

  logical function lst(x,y)
    type(natural_number), intent(in) :: x, y
    lst = cmp(x,y) < 0
  end

  logical function goe(x,y)
    type(natural_number), intent(in) :: x, y
    goe = cmp(x,y) >= 0
  end

  logical function loe(x,y)
    type(natural_number), intent(in) :: x, y
    loe = cmp(x,y) <= 0
  end

  logical function eql(x,y)
    type(natural_number), intent(in) :: x, y
    eql = cmp(x,y) == 0
  end

  logical function neq(x,y)
    type(natural_number), intent(in) :: x, y
    neq = cmp(x,y) /= 0
  end

  function gcd(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n, m, t
    if (cmp(x,y) > 0) then
      call set(n,x)
      call set(m,y)
    else
      call set(n,y)
      call set(m,x)
    end if
    do while (associated(m%nn))
      t = rem(n,m)
      call set(n,m)
      call set(m,t)
    end do
    call norm(n)
  end

  function lcm(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n
    n = div(mul(x,y),gcd(x,y))
  end

  function andn(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n, a, b, p, two
    call set(a,x)
    call set(b,y)
    p = new()
    two = add(p,p)
    do while (associated(a%nn) .and. associated(b%nn))
      if (mod(a%nn(0),2) == 1 .and. mod(b%nn(0),2) == 1) n = add(n,p)
      p = mul(p,two)
      a = half(a)
      b = half(b)
    end do
  end

  function orn(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n, a, b, p, two
    call set(a,x)
    call set(b,y)
    p = new()
    two = add(p,p)
    do while (associated(a%nn) .or. associated(b%nn))
      if ((associated(a%nn) .and. mod(a%nn(0),2) == 1) .or. &
      & (associated(b%nn) .and. mod(b%nn(0),2) == 1)) n = add(n,p)
      p = mul(p,two)
      a = half(a)
      b = half(b)
    end do
  end

  function xorn(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n, a, b, p, two
    call set(a,x)
    call set(b,y)
    p = new()
    two = add(p,p)
    do while (associated(a%nn) .or. associated(b%nn))
      if ((associated(a%nn) .and. mod(a%nn(0),2) == 1) .neqv. &
      & (associated(b%nn) .and. mod(b%nn(0),2) == 1)) n = add(n,p)
      p = mul(p,two)
      a = half(a)
      b = half(b)
    end do
  end

  function maxn(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n
    if (cmp(x,y) > 0) then
      call set(n,x)
    else
      call set(n,y)
    end if
  end

  function minn(x,y) result(n)
    type(natural_number), intent(in) :: x, y
    type(natural_number) :: n
    if (cmp(x,y) < 0) then
      call set(n,x)
    else
      call set(n,y)
    end if
  end

  function dot_productn(x,y) result(n)
    type(natural_number), intent(in) :: x(:), y(:)
    type(natural_number) :: n
    integer :: i
    if (size(x) /= size(y)) then
      write(*,'(a)') "Error: Invalid lengths of arrays (natural_number, dot_productn)"
      stop
    end if
    do i = 1, size(x)
      n = add(n,mul(x(i),y(i)))
    end do
  end

  function mat_mulnmv(x,y) result(n)
    type(natural_number), intent(in) :: x(:,:), y(:)
    type(natural_number) :: n(size(x,1))
    integer :: i
    if (size(x,2) /= size(y)) then
      write(*,'(a)') "Error: Invalid size (natural_number, mat_mulnmv)"
      stop
    end if
    do i = 1, size(x,1)
      n(i) = dot_productn(x(i,:),y)
    end do
  end

  function mat_mulnvm(x,y) result(n)
    type(natural_number), intent(in) :: x(:), y(:,:)
    type(natural_number) :: n(size(y,2))
    integer :: i
    if (size(x) /= size(y,1)) then
      write(*,'(a)') "Error: Invalid size (natural_number, mat_mulnvm)"
      stop
    end if
    do i = 1, size(y,2)
      n(i) = dot_productn(x,y(:,i))
    end do
  end

  function mat_mulnmm(x,y) result(n)
    type(natural_number), intent(in) :: x(:,:), y(:,:)
    type(natural_number) :: n(size(x,1),size(y,2))
    integer :: i
    if (size(x,2) /= size(y,1)) then
      write(*,'(a)') "Error: Invalid size (natural_number, mat_mulnmm)"
      stop
    end if
    do i = 1, size(y,2)
      n(:,i) = mat_mulnmv(x,y(:,i))
    end do
  end

  ! overload with character

  function new0(s) result(n)
    character(*), intent(in) :: s
    type(natural_number) :: n
    character(len(s)) :: t
    integer :: m, i
    t = trim(adjustl(s))
    m = len_trim(t)
    if (m == 0) return
    allocate(n%nn(0:m-1))
    do i = 1, m
      n%nn(m-i) = int(ichar(t(i:i))-48,1)
      if (n%nn(m-i) < 0 .or. 9 < n%nn(m-i)) then
        write(*,'(a)') "Error: Illegal number format (natural_number, new0)"
        stop
      end if
    end do
    call norm(n)
  end

  subroutine set0(n,s)
    type(natural_number), intent(out) :: n
    character(*), intent(in) :: s
    call set(n,new0(s))
  end

  integer function cmp0ns(x,s)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    cmp0ns = cmp(x,new0(s))
  end

  integer function cmp0sn(s,y)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    cmp0sn = cmp(new0(s),y)
  end

  function add0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = add(x,new0(s))
  end

  function add0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = add(new0(s),y)
  end

  function sub0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = sub(x,new0(s))
  end

  function sub0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = sub(new0(s),y)
  end

  function mul0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = mul(x,new0(s))
  end

  function mul0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = mul(new0(s),y)
  end

  function div0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = div(x,new0(s))
  end

  function div0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = div(new0(s),y)
  end

  function rem0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = rem(x,new0(s))
  end

  function rem0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = rem(new0(s),y)
  end

  function pow0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = pow(x,new0(s))
  end

  function pow0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = pow(new0(s),y)
  end

  logical function gtt0ns(x,s)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    gtt0ns = gtt(x,new0(s))
  end

  logical function gtt0sn(s,y)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    gtt0sn = gtt(new0(s),y)
  end

  logical function lst0ns(x,s)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    lst0ns = lst(x,new0(s))
  end

  logical function lst0sn(s,y)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    lst0sn = lst(new0(s),y)
  end

  logical function goe0ns(x,s)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    goe0ns = goe(x,new0(s))
  end

  logical function goe0sn(s,y)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    goe0sn = goe(new0(s),y)
  end

  logical function loe0ns(x,s)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    loe0ns = loe(x,new0(s))
  end

  logical function loe0sn(s,y)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    loe0sn = loe(new0(s),y)
  end

  logical function eql0ns(x,s)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    eql0ns = eql(x,new0(s))
  end

  logical function eql0sn(s,y)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    eql0sn = eql(new0(s),y)
  end

  logical function neq0ns(x,s)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    neq0ns = neq(x,new0(s))
  end

  logical function neq0sn(s,y)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    neq0sn = neq(new0(s),y)
  end

  function gcd0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = gcd(x,new0(s))
  end

  function gcd0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = gcd(new0(s),y)
  end

  function lcm0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = lcm(x,new0(s))
  end

  function lcm0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = lcm(new0(s),y)
  end

  function and0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = andn(x,new0(s))
  end

  function and0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = andn(new0(s),y)
  end

  function or0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = orn(x,new0(s))
  end

  function or0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = orn(new0(s),y)
  end

  function xor0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = xorn(x,new0(s))
  end

  function xor0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = xorn(new0(s),y)
  end

  function max0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = maxn(x,new0(s))
  end

  function max0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = maxn(new0(s),y)
  end

  function min0ns(x,s) result(n)
    type(natural_number), intent(in) :: x
    character(*), intent(in) :: s
    type(natural_number) :: n
    n = minn(x,new0(s))
  end

  function min0sn(s,y) result(n)
    character(*), intent(in) :: s
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = minn(new0(s),y)
  end

  ! overload with primitive integer(8)

  function new1(m) result(n)
    integer(8), intent(in) :: m
    type(natural_number) :: n
    character(20) :: s
    if (m < 0_8) then
      write(*,'(a)') "Error: Overflow (natural_number, new1)"
      stop
    end if
    if (m == 0_8) return
    write(s,'(i0)') m
    n = new0(s)
  end

  subroutine set1(n,m)
    type(natural_number), intent(out) :: n
    integer(8), intent(in) :: m
    call set(n,new1(m))
  end

  integer function cmp1nl(x,m)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    cmp1nl = cmp(x,new1(m))
  end

  integer function cmp1ln(m,y)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    cmp1ln = cmp(new1(m),y)
  end

  function add1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = add(x,new1(m))
  end

  function add1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = add(new1(m),y)
  end

  function sub1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = sub(x,new1(m))
  end

  function sub1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = sub(new1(m),y)
  end

  function mul1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = mul(x,new1(m))
  end

  function mul1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = mul(new1(m),y)
  end

  function div1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = div(x,new1(m))
  end

  function div1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = div(new1(m),y)
  end

  function rem1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = rem(x,new1(m))
  end

  function rem1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = rem(new1(m),y)
  end

  function pow1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = pow(x,new1(m))
  end

  function pow1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = pow(new1(m),y)
  end

  logical function gtt1nl(x,m)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    gtt1nl = gtt(x,new1(m))
  end

  logical function gtt1ln(m,y)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    gtt1ln = gtt(new1(m),y)
  end

  logical function lst1nl(x,m)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    lst1nl = lst(x,new1(m))
  end

  logical function lst1ln(m,y)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    lst1ln = lst(new1(m),y)
  end

  logical function goe1nl(x,m)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    goe1nl = goe(x,new1(m))
  end

  logical function goe1ln(m,y)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    goe1ln = goe(new1(m),y)
  end

  logical function loe1nl(x,m)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    loe1nl = loe(x,new1(m))
  end

  logical function loe1ln(m,y)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    loe1ln = loe(new1(m),y)
  end

  logical function eql1nl(x,m)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    eql1nl = eql(x,new1(m))
  end

  logical function eql1ln(m,y)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    eql1ln = eql(new1(m),y)
  end

  logical function neq1nl(x,m)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    neq1nl = neq(x,new1(m))
  end

  logical function neq1ln(m,y)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    neq1ln = neq(new1(m),y)
  end

  function gcd1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = gcd(x,new1(m))
  end

  function gcd1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = gcd(new1(m),y)
  end

  function lcm1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = lcm(x,new1(m))
  end

  function lcm1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = lcm(new1(m),y)
  end

  function and1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = andn(x,new1(m))
  end

  function and1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = andn(new1(m),y)
  end

  function or1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = orn(x,new1(m))
  end

  function or1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = orn(new1(m),y)
  end

  function xor1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = xorn(x,new1(m))
  end

  function xor1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = xorn(new1(m),y)
  end

  function max1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = maxn(x,new1(m))
  end

  function max1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = maxn(new1(m),y)
  end

  function min1nl(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer(8), intent(in) :: m
    type(natural_number) :: n
    n = minn(x,new1(m))
  end

  function min1ln(m,y) result(n)
    integer(8), intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = minn(new1(m),y)
  end

  ! overload with primitive integer

  function new2(m) result(n)
    integer, intent(in) :: m
    type(natural_number) :: n
    character(10) :: s
    if (m < 0) then
      write(*,'(a)') "Error: Overflow (natural_number, new2)"
      stop
    end if
    if (m == 0) return
    write(s,'(i0)') m
    n = new0(s)
  end

  subroutine set2(n,m)
    type(natural_number), intent(out) :: n
    integer, intent(in) :: m
    call set(n,new2(m))
  end

  integer function cmp2ni(x,m)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    cmp2ni = cmp(x,new2(m))
  end

  integer function cmp2in(m,y)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    cmp2in = cmp(new2(m),y)
  end

  function add2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = add(x,new2(m))
  end

  function add2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = add(new2(m),y)
  end

  function sub2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = sub(x,new2(m))
  end

  function sub2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = sub(new2(m),y)
  end

  function mul2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = mul(x,new2(m))
  end

  function mul2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = mul(new2(m),y)
  end

  function div2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = div(x,new2(m))
  end

  function div2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = div(new2(m),y)
  end

  function rem2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = rem(x,new2(m))
  end

  function rem2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = rem(new2(m),y)
  end

  function pow2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = pow(x,new2(m))
  end

  function pow2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = pow(new2(m),y)
  end

  logical function gtt2ni(x,m)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    gtt2ni = gtt(x,new2(m))
  end

  logical function gtt2in(m,y)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    gtt2in = gtt(new2(m),y)
  end

  logical function lst2ni(x,m)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    lst2ni = lst(x,new2(m))
  end

  logical function lst2in(m,y)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    lst2in = lst(new2(m),y)
  end

  logical function goe2ni(x,m)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    goe2ni = goe(x,new2(m))
  end

  logical function goe2in(m,y)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    goe2in = goe(new2(m),y)
  end

  logical function loe2ni(x,m)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    loe2ni = loe(x,new2(m))
  end

  logical function loe2in(m,y)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    loe2in = loe(new2(m),y)
  end

  logical function eql2ni(x,m)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    eql2ni = eql(x,new2(m))
  end

  logical function eql2in(m,y)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    eql2in = eql(new2(m),y)
  end

  logical function neq2ni(x,m)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    neq2ni = neq(x,new2(m))
  end

  logical function neq2in(m,y)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    neq2in = neq(new2(m),y)
  end

  function gcd2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = gcd(x,new2(m))
  end

  function gcd2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = gcd(new2(m),y)
  end

  function lcm2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = lcm(x,new2(m))
  end

  function lcm2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = lcm(new2(m),y)
  end

  function and2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = andn(x,new2(m))
  end

  function and2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = andn(new2(m),y)
  end

  function or2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = orn(x,new2(m))
  end

  function or2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = orn(new2(m),y)
  end

  function xor2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = xorn(x,new2(m))
  end

  function xor2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = xorn(new2(m),y)
  end

  function max2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = maxn(x,new2(m))
  end

  function max2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = maxn(new2(m),y)
  end

  function min2ni(x,m) result(n)
    type(natural_number), intent(in) :: x
    integer, intent(in) :: m
    type(natural_number) :: n
    n = minn(x,new2(m))
  end

  function min2in(m,y) result(n)
    integer, intent(in) :: m
    type(natural_number), intent(in) :: y
    type(natural_number) :: n
    n = minn(new2(m),y)
  end

end module mod_natural_number
program multiple_of_9
  use mod_natural_number
  implicit none
  character(200000) :: s
  type(natural_number) :: n
  read(*,'(a)') s
  n = natural_number(s)
  write(*,'(a)') trim(merge('Yes', 'No ', mod(n, 9) == 0))
end program multiple_of_9

Measured memory_bytes: 4348