FortranAdvent Calendar 2024

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

Last updated at 2024-12-05Posted at 2024-12-05

やりたいこと

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

多倍長整数を扱いたい

皆さんは「多倍長整数を扱いたい!」と思ったことはあるでしょうか?
多倍長整数は言語によっては標準で対応しているものもあります.
(python, HaskellのInteger, JavaのBigInteger などなど...)

需要としては...

wikipedia によると公開鍵, オーバーフロー対策があるそうです.
競技プログラミングでは多倍長整数を扱えると有利になるときがあったりします.

実装方法

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

Fortran の実装は?

FM package

詳細はよく知りませんが, 多倍長整数や任意精度実数ができそうです.
GMP 使わなくても良いかもしれません(完).

実行環境

$ 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

とりあえず, iso_c_binding の練習として, Fortran から GMP を利用します.

ライブラリの使い方

wikipediaを参考にすると

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 は領域の解放をします.

実行結果

gcc test_gmp.c -lgmp && ./a.out
123456789101112131415161718
*
987654321191817161514131211
------------------------------
121932631236180233352728375173925092673910835836180498

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日目によると nm コマンドでバイナリの中身が見れるらしいです.

$ 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 の処理です.
マニュアルの 16 Internals に書いてあることを読むと,

type, bind(c) :: mpz_t
   integer(c_int) :: mp_size
   type(c_ptr) :: mp_d
   integer(c_int) :: mp_alloc
end type mpz_t

でいけそうです.

ただ, 内部のことは変わるかもしれないと書いてあります. 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.

test_gmp.f90 全体

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 側が変数が変更されることを知らない?