96 lines
2.5 KiB
Fortran
96 lines
2.5 KiB
Fortran
module modular_arithmetic
|
|
implicit none
|
|
|
|
type :: modular
|
|
integer :: val
|
|
integer :: modulus
|
|
end type modular
|
|
|
|
interface operator(+)
|
|
module procedure modular_modular_add
|
|
module procedure modular_integer_add
|
|
end interface operator(+)
|
|
|
|
interface operator(**)
|
|
module procedure modular_integer_pow
|
|
end interface operator(**)
|
|
|
|
contains
|
|
|
|
function modular_modular_add (a, b) result (c)
|
|
type(modular), intent(in) :: a
|
|
type(modular), intent(in) :: b
|
|
type(modular) :: c
|
|
|
|
if (a%modulus /= b%modulus) error stop
|
|
c%val = modulo (a%val + b%val, a%modulus)
|
|
c%modulus = a%modulus
|
|
end function modular_modular_add
|
|
|
|
function modular_integer_add (a, i) result (c)
|
|
type(modular), intent(in) :: a
|
|
integer, intent(in) :: i
|
|
type(modular) :: c
|
|
|
|
c%val = modulo (a%val + i, a%modulus)
|
|
c%modulus = a%modulus
|
|
end function modular_integer_add
|
|
|
|
function modular_integer_pow (a, i) result (c)
|
|
type(modular), intent(in) :: a
|
|
integer, intent(in) :: i
|
|
type(modular) :: c
|
|
|
|
! One cannot simply use the integer ** operator and then compute
|
|
! the least residue, because the integers will overflow. Let us
|
|
! instead use the right-to-left binary method:
|
|
! https://en.wikipedia.org/w/index.php?title=Modular_exponentiation&oldid=1136216610#Right-to-left_binary_method
|
|
|
|
integer :: modulus
|
|
integer :: base
|
|
integer :: exponent
|
|
|
|
modulus = a%modulus
|
|
exponent = i
|
|
|
|
if (modulus < 1) error stop
|
|
if (exponent < 0) error stop
|
|
|
|
c%modulus = modulus
|
|
if (modulus == 1) then
|
|
c%val = 0
|
|
else
|
|
c%val = 1
|
|
base = modulo (a%val, modulus)
|
|
do while (exponent > 0)
|
|
if (modulo (exponent, 2) /= 0) then
|
|
c%val = modulo (c%val * base, modulus)
|
|
end if
|
|
exponent = exponent / 2
|
|
base = modulo (base * base, modulus)
|
|
end do
|
|
end if
|
|
end function modular_integer_pow
|
|
|
|
end module modular_arithmetic
|
|
|
|
! If one uses the extension .F90 instead of .f90, then gfortran will
|
|
! pass the program through the C preprocessor. Thus one can write f(x)
|
|
! without considering the type of the argument
|
|
#define f(x) ((x)**100 + (x) + 1)
|
|
|
|
program modular_arithmetic_task
|
|
use, intrinsic :: iso_fortran_env
|
|
use, non_intrinsic :: modular_arithmetic
|
|
implicit none
|
|
|
|
type(modular) :: x, y
|
|
|
|
x = modular(10, 13)
|
|
y = f(x)
|
|
|
|
write (*, '(" modulus 13: ", I0)') y%val
|
|
write (*, '("floating point: ", E55.50)') f(10.0_real64)
|
|
|
|
end program modular_arithmetic_task
|