Fortran から多倍長整数をGMP経由で扱ってみる

# やりたいこと
- Fortran ではできないこと, 面倒なことを C 言語を通してやりたい!

# 多倍長整数を扱いたい
皆さんは「多倍長整数を扱いたい!」と思ったことはあるでしょうか?
多倍長整数は言語によっては標準で対応しているものもあります.
(python, HaskellのInteger, JavaのBigInteger などなど...)
## 需要としては...
- wikipedia によると [公開鍵, オーバーフロー対策](https://ja.wikipedia.org/wiki/%E4%BB%BB%E6%84%8F%E7%B2%BE%E5%BA%A6%E6%BC%94%E7%AE%97#%E7%94%A8%E9%80%94) があるそうです.
- 競技プログラミングでは多倍長整数を扱えると有利になるときがあったりします.


## 実装方法
[ここ](https://qiita.com/square1001/items/1aa12e04934b6e749962)を参考にすると, 10進数で $N$ 桁の数を表現するには, サイズ $N$ の配列を用います.
足し算や引き算は各々の配列の各要素同士を足し引きして, その後に繰り上がりを処理すれば $O(N)$ です.
厄介なのが掛け算で, 単純に実装すると $O(N^2)$ かかります.
カラツバ法や高速フーリエ変換を実装するのは箸を持つより重いため, 迂回するのも選択肢としてはアリかもしれません.

## Fortran の実装は?
### [FM package](https://dmsmith.lmu.build/)
詳細はよく知りませんが, 多倍長整数や任意精度実数ができそうです.
GMP 使わなくても良いかもしれません(完).

### [他にもいろいろあります...](https://fortranwiki.org/fortran/show/Libraries)
- [MPFUN2020](https://www.davidhbailey.com/dhbsoftware/)
- [gmp-fortran-library](https://sourceforge.net/projects/gmp-fortran-library/)


# 実行環境
```bash
$ lsb_release -a
No LSB modules are available.
Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.6 LTS
Release:	20.04
Codename:	focal
$ gfortran --version
GNU Fortran (GCC) 9.2.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
$ gcc --version
gcc (GCC) 9.2.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
$ apt list libgmp10
一覧表示... 完了
libgmp10/focal-updates,focal-security,now 2:6.2.0+dfsg-4ubuntu0.1 amd64 [インストール済み、自動]
libgmp10/focal-updates,focal-security 2:6.2.0+dfsg-4ubuntu0.1 i386
```

# [GNU Multi-Precision Library](https://gmplib.org/)
とりあえず, `iso_c_binding` の練習として, Fortran から GMP を利用します.
## ライブラリの使い方
[wikipedia](https://ja.wikipedia.org/wiki/GNU_Multi-Precision_Library)を参考にすると

```test_gmp.c
#include <gmp.h>

int main(void)
{
  mpz_t x, y, res;

  mpz_init_set_str(x, "123456789101112131415161718", 10);
  mpz_init_set_str(y, "987654321191817161514131211", 10);
  mpz_init(res);

  mpz_mul(res, x, y);
  gmp_printf("%Zd\n*\n%Zd\n------------------------------\n%Zd\n", x, y, res);

  mpz_clear(x);
  mpz_clear(y);
  mpz_clear(res);

  return 0;
}
```

- `mpz_t` が多倍長整数を表す構造体の型です.
- `mpz_init_set_str` や `mpz_init` で文字列で表される多倍長整数や `0` を表現できます.
  + `mpz_init_set_str` は初期化と代入をする関数です.
  + `10` 進数にしています.
- `mpz_mul` は掛け算します.
- `gmp_printf` でフォーマット付きの出力ができます.
  + `%Z` は `mpz_t` を出力できます.
- `mpz_clear` は領域の解放をします.

実行結果
```bash
gcc test_gmp.c -lgmp && ./a.out
123456789101112131415161718
*
987654321191817161514131211
------------------------------
121932631236180233352728375173925092673910835836180498
```

- Haskell
```haskell
$ ghci
GHCi, version 8.6.5: http://www.haskell.org/ghc/  :? for help
Prelude> 123456789101112131415161718
123456789101112131415161718
Prelude> 123456789101112131415161718 * 987654321191817161514131211
121932631236180233352728375173925092673910835836180498
```

計算できています.

## Fortran から呼び出してみる
[Fortran アドベントカレンダー 2024 1日目](https://qiita.com/sakamoti/items/0fdc989ee530a5639f55) によると `nm` コマンドでバイナリの中身が見れるらしいです.
```bash
$ nm /usr/lib/x86_64-linux-gnu/libgmp.a
(略)
bin_ui.o:
                 U _GLOBAL_OFFSET_TABLE_
                 U __gmpn_rshift
                 U __gmpz_add
                 U __gmpz_add_ui
                 U __gmpz_addmul_ui
0000000000000350 T __gmpz_bin_ui
                 U __gmpz_clear
                 U __gmpz_cmp_ui
                 U __gmpz_divexact
                 U __gmpz_init
                 U __gmpz_init_set
                 U __gmpz_init_set_ui
                 U __gmpz_mul
(略)
printf.o:
                 U _GLOBAL_OFFSET_TABLE_
                 U __gmp_doprnt
                 U __gmp_fprintf_funs
0000000000000000 T __gmp_printf
                 U __stack_chk_fail
                 U stdout
(略)
```
どうやら, 関数の名前は `__gmp_z_init` みたいに `__` が付いているようです.
関数名が分かったので, 次は構造体の `mpz_t` の処理です.
[マニュアル](https://gmplib.org/gmp-man-6.3.0.pdf) の **16 Internals** に書いてあることを読むと,
```fortran
type, bind(c) :: mpz_t
   integer(c_int) :: mp_size
   type(c_ptr) :: mp_d
   integer(c_int) :: mp_alloc
end type mpz_t
```
でいけそうです.

:::note warn
ただ, 内部のことは変わるかもしれないと書いてあります. https://gmplib.org/gmp-man-6.3.0.pdf
> This chapter is provided only for informational purposes and the various internals described
here may change in future GMP releases. Applications expecting to be compatible with future
releases should use only the documented interfaces described in previous chapters.
:::

<details open><summary> test_gmp.f90 全体 </summary>

```test_gmp.f90
module gmp_m
  use, intrinsic :: iso_fortran_env
  use, intrinsic :: iso_c_binding
  implicit none
  type, bind(c) :: mpz_t
     integer(c_int) :: mp_size
     type(c_ptr) :: mp_d
     integer(c_int) :: mp_alloc
  end type mpz_t
  interface
     subroutine mpz_init_set_str(v, str, base) bind(c, name = "__gmpz_init_set_str")
       import mpz_t, c_char, c_int
       type(mpz_t), intent(inout) :: v
       character(c_char), intent(in) :: str
       integer(c_int), value :: base
     end subroutine mpz_init_set_str
  end interface
  interface
     subroutine mpz_init(v) bind(c, name = "__gmpz_init")
       import mpz_t
       type(mpz_t), intent(inout) :: v
     end subroutine mpz_init
  end interface
  interface
     subroutine mpz_mul(x, y, z) bind(c, name = "__gmpz_mul")
       import mpz_t
       type(mpz_t), intent(inout) :: x
       type(mpz_t), intent(in) :: y, z
     end subroutine mpz_mul
  end interface
  interface
     subroutine gmp_printf(fmt, x, y, z) bind(c, name = "__gmp_printf")
       import c_char, mpz_t
       character(c_char), intent(in) :: fmt
       type(mpz_t), intent(in) :: x, y, z
     end subroutine gmp_printf
  end interface
  interface
     subroutine mpz_clear(v) bind(c, name = "__gmpz_clear")
       import mpz_t
       type(mpz_t), intent(inout) :: v
     end subroutine mpz_clear
  end interface
end module gmp_m

program test_gmp
  use, intrinsic :: iso_fortran_env
  use, intrinsic :: iso_c_binding
  use gmp_m
  implicit none
  type(mpz_t) :: x, y, res

  call mpz_init_set_str(x, "123456789101112131415161718" // c_null_char, 10)
  call mpz_init_set_str(y, "987654321191817161514131211" // c_null_char, 10)
  call mpz_init(res)

  call mpz_mul(res, x, y)
  call gmp_printf("%Zd\n*\n%Zd\n------------------------------\n%Zd\n", x, y, res)

  call mpz_clear(x)
  call mpz_clear(y)
  call mpz_clear(res)
end program test_gmp
```

- `bind(c, name = "__gmpz_init")` のように, Cの名前(見えてる名前?)を指定します.
- `intent(inout)` にすべきか `value` にすべきかを決めます.
  + `type(mpz_t), intent(in)` でも通ってしまいます.
    - やはり安全ではない.
    - Fortran 側が変数が変更されることを知らない?


</details>

実行結果
```bash
 $ gfortran -fbackslash test_gmp.f90 -lgmp && ./a.out
123456789101112131415161718
*
987654321191817161514131211
------------------------------
121932631236180233352728375173925092673910835836180498
```

Fortran でもGMPを使うことができました.

# まとめ
今回は, GMP を Fortran で定義した `type(mpz_t)` を通じて扱うことができました.
他の演算も同様に呼び出すことが, おそらく, 可能でしょう.
ただ, コンパイラのチェックが効いていないようなので, C 言語並の安全性を手に入れてしまうのが懸念点ですね.
実用的には, 2024年にも更新されている `FM package` や `MPFUN2020` を使えば良いのだと思います.

# 参考

https://gmplib.org/

https://ja.wikipedia.org/wiki/%E4%BB%BB%E6%84%8F%E7%B2%BE%E5%BA%A6%E6%BC%94%E7%AE%97#%E7%94%A8%E9%80%94

https://ja.wikipedia.org/wiki/GNU_Multi-Precision_Library

https://dmsmith.lmu.build/

https://www.davidhbailey.com/dhbsoftware/

https://sourceforge.net/projects/gmp-fortran-library/

https://fortranwiki.org/fortran/show/Libraries

https://qiita.com/square1001/items/1aa12e04934b6e749962