Continued fraction/Arithmetic/G(matrix ng, continued fraction n1, continued fraction n2): Difference between revisions

Content added Content deleted
Line 2,626: Line 2,626:
sqrt(2) * sqrt(2) => [2]
sqrt(2) * sqrt(2) => [2]
sqrt(2) / sqrt(2) => [1]
sqrt(2) / sqrt(2) => [1]
</pre>

=={{header|Fortran}}==
{{trans|Python}}
{{trans|ObjectIcon}}

This program includes a primitive module for multiple-precision integer arithmetic.

<syntaxhighlight lang="fortran">
!---------------------------------------------------------------------

module big_integers ! Big (but not very big) integers.

! NOTE: I assume that iachar and achar do not alter the most
! significant bit.

use, intrinsic :: iso_fortran_env, only: int16
implicit none
private

public :: big_integer
public :: integer2big
public :: string2big
public :: big2string
public :: big_sgn
public :: big_cmp, big_cmpabs
public :: big_neg, big_abs
public :: big_addabs, big_add
public :: big_subabs, big_sub
public :: big_mul ! One might also include a big_muladd.
public :: big_divrem ! One could also include big_div and big_rem.
public :: operator(+)
public :: operator(-)
public :: operator(*)

type :: big_integer
! The representation is sign-magnitude. The radix is 256, which
! is not speed-efficient, but which seemed relatively easy to
! work with if one were writing in standard Fortran (and assuming
! iachar and achar were "8-bit clean").
logical :: sign = .false. ! .false. => +sign, .true. => -sign.
character, allocatable :: bytes(:)
end type big_integer

character, parameter :: zero = achar (0)
character, parameter :: one = achar (1)

! An integer type capable of holding an unsigned 8-bit value.
integer, parameter :: bytekind = int16

interface operator(+)
module procedure big_add
end interface

interface operator(-)
module procedure big_neg
module procedure big_sub
end interface

interface operator(*)
module procedure big_mul
end interface

contains

elemental function logical2byte (bool) result (byte)
logical, intent(in) :: bool
character :: byte
if (bool) then
byte = one
else
byte = zero
end if
end function logical2byte

elemental function logical2i (bool) result (i)
logical, intent(in) :: bool
integer :: i
if (bool) then
i = 1
else
i = 0
end if
end function logical2i

elemental function byte2i (c) result (i)
character, intent(in) :: c
integer :: i
i = iachar (c)
end function byte2i

elemental function i2byte (i) result (c)
integer, intent(in) :: i
character :: c
c = achar (i)
end function i2byte

elemental function byte2bk (c) result (i)
character, intent(in) :: c
integer(bytekind) :: i
i = iachar (c, kind = bytekind)
end function byte2bk

elemental function bk2byte (i) result (c)
integer(bytekind), intent(in) :: i
character :: c
c = achar (i)
end function bk2byte

elemental function bk2i (i) result (j)
integer(bytekind), intent(in) :: i
integer :: j
j = int (i)
end function bk2i

elemental function i2bk (i) result (j)
integer, intent(in) :: i
integer(bytekind) :: j
j = int (iand (i, 255), kind = bytekind)
end function i2bk

! Left shift of the least significant 8 bits of a bytekind integer.
elemental function lshftbk (a, i) result (c)
integer(bytekind), intent(in) :: a
integer, intent(in) :: i
integer(bytekind) :: c
c = ishft (ibits (a, 0, 8 - i), i)
end function lshftbk

! Right shift of the least significant 8 bits of a bytekind integer.
elemental function rshftbk (a, i) result (c)
integer(bytekind), intent(in) :: a
integer, intent(in) :: i
integer(bytekind) :: c
c = ibits (a, i, 8 - i)
end function rshftbk

! Left shift an integer.
elemental function lshfti (a, i) result (c)
integer, intent(in) :: a
integer, intent(in) :: i
integer :: c
c = ishft (a, i)
end function lshfti

! Right shift an integer.
elemental function rshfti (a, i) result (c)
integer, intent(in) :: a
integer, intent(in) :: i
integer :: c
c = ishft (a, -i)
end function rshfti

function integer2big (i) result (a)
integer, intent(in) :: i
type(big_integer), allocatable :: a

!
! To write a more efficient implementation of this procedure is
! left as an exercise for the reader.
!

character(len = 100) :: buffer

write (buffer, '(I0)') i
a = string2big (trim (buffer))
end function integer2big

function string2big (s) result (a)
character(len = *), intent(in) :: s
type(big_integer), allocatable :: a

integer :: n, i, istart, iend
integer :: digit

if ((s(1:1) == '-') .or. s(1:1) == '+') then
istart = 2
else
istart = 1
end if

iend = len (s)

n = (iend - istart + 2) / 2

allocate (a)
allocate (a%bytes(n))

a%bytes = zero
do i = istart, iend
digit = ichar (s(i:i)) - ichar ('0')
if (digit < 0 .or. 9 < digit) error stop
a = short_multiplication (a, 10)
a = short_addition (a, digit)
end do
a%sign = (s(1:1) == '-')
call normalize (a)
end function string2big

function big2string (a) result (s)
type(big_integer), intent(in) :: a
character(len = :), allocatable :: s

type(big_integer), allocatable :: q
integer :: r
integer :: sgn

sgn = big_sgn (a)
if (sgn == 0) then
s = '0'
else
q = a
s = ''
do while (big_sgn (q) /= 0)
call short_division (q, 10, q, r)
s = achar (r + ichar ('0')) // s
end do
if (sgn < 0) s = '-' // s
end if
end function big2string

function big_sgn (a) result (sgn)
type(big_integer), intent(in) :: a
integer :: sgn

integer :: n, i

n = size (a%bytes)
i = 1
sgn = 1234
do while (sgn == 1234)
if (i == n + 1) then
sgn = 0
else if (a%bytes(i) /= zero) then
if (a%sign) then
sgn = -1
else
sgn = 1
end if
else
i = i + 1
end if
end do
end function big_sgn

function big_cmp (a, b) result (cmp)
type(big_integer(*)), intent(in) :: a, b
integer :: cmp

if (a%sign) then
if (b%sign) then
cmp = -big_cmpabs (a, b)
else
cmp = -1
end if
else
if (b%sign) then
cmp = 1
else
cmp = big_cmpabs (a, b)
end if
end if
end function big_cmp

function big_cmpabs (a, b) result (cmp)
type(big_integer(*)), intent(in) :: a, b
integer :: cmp

integer :: n, i
integer :: ia, ib

cmp = 1234
n = max (size (a%bytes), size (b%bytes))
i = n
do while (cmp == 1234)
if (i == 0) then
cmp = 0
else
ia = byteval (a, i)
ib = byteval (b, i)
if (ia < ib) then
cmp = -1
else if (ia > ib) then
cmp = 1
else
i = i - 1
end if
end if
end do
end function big_cmpabs

function big_neg (a) result (c)
type(big_integer), intent(in) :: a
type(big_integer), allocatable :: c
c = a
c%sign = .not. c%sign
end function big_neg

function big_abs (a) result (c)
type(big_integer), intent(in) :: a
type(big_integer), allocatable :: c
c = a
c%sign = .false.
end function big_abs

function big_add (a, b) result (c)
type(big_integer), intent(in) :: a
type(big_integer), intent(in) :: b
type(big_integer), allocatable :: c

logical :: sign

if (a%sign) then
if (b%sign) then ! a <= 0, b <= 0
c = big_addabs (a, b)
sign = .true.
else ! a <= 0, b >= 0
c = big_subabs (a, b)
sign = .not. c%sign
end if
else
if (b%sign) then ! a >= 0, b <= 0
c = big_subabs (a, b)
sign = c%sign
else ! a >= 0, b >= 0
c = big_addabs (a, b)
sign = .false.
end if
end if
c%sign = sign
end function big_add

function big_sub (a, b) result (c)
type(big_integer), intent(in) :: a
type(big_integer), intent(in) :: b
type(big_integer), allocatable :: c

logical :: sign

if (a%sign) then
if (b%sign) then ! a <= 0, b <= 0
c = big_subabs (a, b)
sign = .not. c%sign
else ! a <= 0, b >= 0
c = big_addabs (a, b)
sign = .true.
end if
else
if (b%sign) then ! a >= 0, b <= 0
c = big_addabs (a, b)
sign = .false.
else ! a >= 0, b >= 0
c = big_subabs (a, b)
sign = c%sign
end if
end if
c%sign = sign
end function big_sub

function big_addabs (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c

! Compute abs(a) + abs(b).

integer :: n, nc, i
logical :: carry
type(big_integer), allocatable :: tmp

n = max (size (a%bytes), size (b%bytes))
nc = n + 1

allocate(tmp)
allocate(tmp%bytes(nc))

call add_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call add_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry)
end do
tmp%bytes(nc) = logical2byte (carry)
call normalize (tmp)
c = tmp
end function big_addabs

function big_subabs (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c

! Compute abs(a) - abs(b). The result is signed.

integer :: n, i
logical :: carry
type(big_integer), allocatable :: tmp

n = max (size (a%bytes), size (b%bytes))
allocate(tmp)
allocate(tmp%bytes(n))

if (big_cmpabs (a, b) >= 0) then
tmp%sign = .false.
call sub_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call sub_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry)
end do
else
tmp%sign = .true.
call sub_bytes (get_byte (b, 1), get_byte (a, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call sub_bytes (get_byte (b, i), get_byte (a, i), carry, tmp%bytes(i), carry)
end do
end if
call normalize (tmp)
c = tmp
end function big_subabs

function big_mul (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c

!
! This is Knuth, Volume 2, Algorithm 4.3.1M.
!

integer :: na, nb, nc
integer :: i, j
integer :: ia, ib, ic
integer :: carry
type(big_integer), allocatable :: tmp

na = size (a%bytes)
nb = size (b%bytes)
nc = na + nb + 1

allocate (tmp)
allocate (tmp%bytes(nc))

tmp%bytes = zero
j = 1
do j = 1, nb
ib = byte2i (b%bytes(j))
if (ib /= 0) then
carry = 0
do i = 1, na
ia = byte2i (a%bytes(i))
ic = byte2i (tmp%bytes(i + j - 1))
ic = (ia * ib) + ic + carry
tmp%bytes(i + j - 1) = i2byte (iand (ic, 255))
carry = ishft (ic, -8)
end do
tmp%bytes(na + j) = i2byte (carry)
end if
end do
tmp%sign = (a%sign .neqv. b%sign)
call normalize (tmp)
c = tmp
end function big_mul

subroutine big_divrem (a, b, q, r)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable, intent(inout) :: q, r

!
! Division with a remainder that is never negative. Equivalently,
! this is floor division if the divisor is positive, and ceiling
! division if the divisor is negative.
!
! See Raymond T. Boute, "The Euclidean definition of the functions
! div and mod", ACM Transactions on Programming Languages and
! Systems, Volume 14, Issue 2, pp. 127-144.
! https://doi.org/10.1145/128861.128862
!

call nonnegative_division (a, b, .true., .true., q, r)
if (a%sign) then
if (big_sgn (r) /= 0) then
q = short_addition (q, 1)
r = big_sub (big_abs (b), r)
end if
q%sign = .not. b%sign
else
q%sign = b%sign
end if
end subroutine big_divrem

function short_addition (a, b) result (c)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), allocatable :: c

! Compute abs(a) + b.

integer :: na, nc, i
logical :: carry
type(big_integer), allocatable :: tmp

na = size (a%bytes)
nc = na + 1

allocate(tmp)
allocate(tmp%bytes(nc))

call add_bytes (a%bytes(1), i2byte (b), .false., tmp%bytes(1), carry)
do i = 2, na
call add_bytes (a%bytes(i), zero, carry, tmp%bytes(i), carry)
end do
tmp%bytes(nc) = logical2byte (carry)
call normalize (tmp)
c = tmp
end function short_addition

function short_multiplication (a, b) result (c)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), allocatable :: c

integer :: i, na, nc
integer :: ia, ic
integer :: carry
type(big_integer), allocatable :: tmp

na = size (a%bytes)
nc = na + 1

allocate (tmp)
allocate (tmp%bytes(nc))

tmp%sign = a%sign
carry = 0
do i = 1, na
ia = byte2i (a%bytes(i))
ic = (ia * b) + carry
tmp%bytes(i) = i2byte (iand (ic, 255))
carry = ishft (ic, -8)
end do
tmp%bytes(nc) = i2byte (carry)
call normalize (tmp)
c = tmp
end function short_multiplication

! Division without regard to signs.
subroutine nonnegative_division (a, b, want_q, want_r, q, r)
type(big_integer), intent(in) :: a, b
logical, intent(in) :: want_q, want_r
type(big_integer), intent(inout), allocatable :: q, r

integer :: na, nb
integer :: remainder

na = size (a%bytes)
nb = size (b%bytes)

! It is an error if b has "significant" zero-bytes or is equal to
! zero.
if (b%bytes(nb) == zero) error stop

if (nb == 1) then
if (want_q) then
call short_division (a, byte2i (b%bytes(1)), q, remainder)
else
block
type(big_integer), allocatable :: bit_bucket
call short_division (a, byte2i (b%bytes(1)), bit_bucket, remainder)
end block
end if
if (want_r) then
if (allocated (r)) deallocate (r)
allocate (r)
allocate (r%bytes(1))
r%bytes(1) = i2byte (remainder)
end if
else
if (na >= nb) then
call long_division (a, b, want_q, want_r, q, r)
else
if (want_q) q = string2big ("0")
if (want_r) r = a
end if
end if
end subroutine nonnegative_division

subroutine short_division (a, b, q, r)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), intent(inout), allocatable :: q
integer, intent(inout) :: r

!
! This is Knuth, Volume 2, Exercise 4.3.1.16.
!
! The divisor is assumed to be positive.
!

integer :: n, i
integer :: ia, ib, iq
type(big_integer), allocatable :: tmp

ib = b
n = size (a%bytes)

allocate (tmp)
allocate (tmp%bytes(n))

r = 0
do i = n, 1, -1
ia = (256 * r) + byte2i (a%bytes(i))
iq = ia / ib
r = mod (ia, ib)
tmp%bytes(i) = i2byte (iq)
end do
tmp%sign = a%sign
call normalize (tmp)
q = tmp
end subroutine short_division

subroutine long_division (a, b, want_quotient, want_remainder, quotient, remainder)
type(big_integer), intent(in) :: a, b
logical, intent(in) :: want_quotient, want_remainder
type(big_integer), intent(inout), allocatable :: quotient
type(big_integer), intent(inout), allocatable :: remainder

!
! This is Knuth, Volume 2, Algorithm 4.3.1D.
!
! We do not deal here with the signs of the inputs and outputs.
!
! It is assumed size(a%bytes) >= size(b%bytes), and that b has no
! leading zero-bytes and is at least two bytes long. If b is one
! byte long and nonzero, use short division.
!

!
! FIXME: THIS ROUTINE IS COMPLICATED AND HAS BRANCHES THAT ARE
! SELDOM RUN. IT HAS NOT BEEN EXTENSIVELY TESTED. Unit
! tests and bugfixes for it would be appreciated.
!

integer :: na, nb, m, n
integer :: num_lz, num_nonlz
integer :: j
integer :: qhat
integer :: carry

!
! We will NOT be working with VERY large numbers, and so it will
! be safe to put temporary storage on the stack. (Note: your
! Fortran might put this storage in a heap instead of the stack.)
!
! v = b, normalized to put its most significant 1-bit all the
! way left.
!
! u = a, shifted left by the same amount as b.
!
! q = the quotient.
!
! The remainder, although shifted left, will end up in u.
!
integer(bytekind) :: u(0:size (a%bytes) + size (b%bytes))
integer(bytekind) :: v(0:size (b%bytes) - 1)
integer(bytekind) :: q(0:size (a%bytes) - size (b%bytes))

na = size (a%bytes)
nb = size (b%bytes)

n = nb
m = na - nb

! In the most significant byte of the divisor, find the number of
! leading zero bits, and the number of bits after that.
block
integer(bytekind) :: tmp
tmp = byte2bk (b%bytes(n))
num_nonlz = bit_size (tmp) - leadz (tmp)
num_lz = 8 - num_nonlz
end block

call normalize_v (b%bytes) ! Make the most significant bit of v be one.
call normalize_u (a%bytes) ! Shifted by the same amount as v.

! Assure ourselves that the most significant bit of v is a one.
if (.not. btest (v(n - 1), 7)) error stop

do j = m, 0, -1
call calculate_qhat (qhat)
call multiply_and_subtract (carry)
q(j) = i2bk (qhat)
if (carry /= 0) call add_back
end do

if (want_quotient) then
if (allocated (quotient)) deallocate (quotient)
allocate (quotient)
allocate (quotient%bytes(m + 1))
quotient%bytes = bk2byte (q)
call normalize (quotient)
end if

if (want_remainder) then
if (allocated (remainder)) deallocate (remainder)
allocate (remainder)
allocate (remainder%bytes(n))
call unnormalize_u (remainder%bytes)
call normalize (remainder)
end if

contains

subroutine normalize_v (b_bytes)
character, intent(in) :: b_bytes(n)

!
! Normalize v so its most significant bit is a one. Any
! normalization factor that achieves this goal will suffice; we
! choose 2**num_lz. (Knuth uses (2**32) div (y[n-1] + 1).)
!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!

integer :: i
integer(bytekind) :: btmp(0:n - 1)

btmp = byte2bk (b_bytes)

v(0) = lshftbk (btmp(0), num_lz)
do i = 1, n - 1
v(i) = ior (lshftbk (btmp(i), num_lz), &
& rshftbk (btmp(i - 1), num_nonlz))
end do
end subroutine normalize_v

subroutine normalize_u (a_bytes)
character, intent(in) :: a_bytes(m + n)

!
! Shift a leftwards to get u. Shift by as much as b was shifted
! to get v.
!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!

integer :: i
integer(bytekind) :: atmp(0:m + n - 1)

atmp = byte2bk (a_bytes)

u(0) = lshftbk (atmp(0), num_lz)
do i = 1, m + n - 1
u(i) = ior (lshftbk (atmp(i), num_lz), &
& rshftbk (atmp(i - 1), num_nonlz))
end do
u(m + n) = rshftbk (atmp(m + n - 1), num_nonlz)
end subroutine normalize_u

subroutine unnormalize_u (r_bytes)
character, intent(out) :: r_bytes(n)

!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!

integer :: i
integer(bytekind) :: rtmp(0:n - 1)

do i = 0, n - 1
rtmp(i) = ior (rshftbk (u(i), num_lz), &
& lshftbk (u(i + 1), num_nonlz))
end do
rtmp(n - 1) = rshftbk (u(n - 1), num_lz)

r_bytes = bk2byte (rtmp)
end subroutine unnormalize_u

subroutine calculate_qhat (qhat)
integer, intent(out) :: qhat

integer :: itmp, rhat
logical :: adjust

itmp = ior (lshfti (bk2i (u(j + n)), 8), &
& bk2i (u(j + n - 1)))
qhat = itmp / bk2i (v(n - 1))
rhat = mod (itmp, bk2i (v(n - 1)))
adjust = .true.
do while (adjust)
if (rshfti (qhat, 8) /= 0) then
continue
else if (qhat * bk2i (v(n - 2)) &
& > ior (lshfti (rhat, 8), &
& bk2i (u(j + n - 2)))) then
continue
else
adjust = .false.
end if
if (adjust) then
qhat = qhat - 1
rhat = rhat + bk2i (v(n - 1))
if (rshfti (rhat, 8) == 0) then
adjust = .false.
end if
end if
end do
end subroutine calculate_qhat

subroutine multiply_and_subtract (carry)
integer, intent(out) :: carry

integer :: i
integer :: qhat_v
integer :: mul_carry, sub_carry
integer :: diff

mul_carry = 0
sub_carry = 0
do i = 0, n
! Multiplication.
qhat_v = mul_carry
if (i /= n) qhat_v = qhat_v + (qhat * bk2i (v(i)))
mul_carry = rshfti (qhat_v, 8)
qhat_v = iand (qhat_v, 255)

! Subtraction.
diff = bk2i (u(j + i)) - qhat_v + sub_carry
sub_carry = -(logical2i (diff < 0)) ! Carry 0 or -1.
u(j + i) = i2bk (diff)
end do
carry = sub_carry
end subroutine multiply_and_subtract

subroutine add_back
integer :: i, carry, sum

q(j) = q(j) - 1_bytekind
carry = 0
do i = 0, n - 1
sum = bk2i (u(i)) + bk2i (v(i)) + carry
carry = ishft (sum, -8)
u(i) = i2bk (sum)
end do
end subroutine add_back

end subroutine long_division

subroutine add_bytes (a, b, carry_in, c, carry_out)
character, intent(in) :: a, b
logical, value :: carry_in
character, intent(inout) :: c
logical, intent(inout) :: carry_out

integer :: ia, ib, ic

ia = byte2i (a)
if (carry_in) ia = ia + 1
ib = byte2i (b)
ic = ia + ib
c = i2byte (iand (ic, 255))
carry_out = (ic >= 256)
end subroutine add_bytes

subroutine sub_bytes (a, b, carry_in, c, carry_out)
character, intent(in) :: a, b
logical, value :: carry_in
character, intent(inout) :: c
logical, intent(inout) :: carry_out

integer :: ia, ib, ic

ia = byte2i (a)
ib = byte2i (b)
if (carry_in) ib = ib + 1
ic = ia - ib
carry_out = (ic < 0)
if (carry_out) ic = ic + 256
c = i2byte (iand (ic, 255))
end subroutine sub_bytes

function get_byte (a, i) result (byte)
type(big_integer), intent(in) :: a
integer, intent(in) :: i
character :: byte

if (size (a%bytes) < i) then
byte = zero
else
byte = a%bytes(i)
end if
end function get_byte

function byteval (a, i) result (v)
type(big_integer), intent(in) :: a
integer, intent(in) :: i
integer :: v

if (size (a%bytes) < i) then
v = 0
else
v = byte2i (a%bytes(i))
end if
end function byteval

subroutine normalize (a)
type(big_integer), intent(inout) :: a

integer :: i
character, allocatable :: fewer_bytes(:)

! Shorten to the minimum number of bytes.
i = size (a%bytes)
do while (1 < i .and. a%bytes(i) == zero)
i = i - 1
end do
if (i /= size (a%bytes)) then
allocate (fewer_bytes (i))
fewer_bytes = a%bytes(1:i)
call move_alloc (fewer_bytes, a%bytes)
end if

! If the magnitude is zero, then clear the sign bit.
if (size (a%bytes) == 1) then
if (a%bytes(1) == zero) then
a%sign = .false.
end if
end if
end subroutine normalize

end module big_integers

!---------------------------------------------------------------------

module continued_fractions

use, non_intrinsic :: big_integers
implicit none
private

public :: continued_fraction

public :: term_generator
public :: term_generator_procedure

public :: make_continued_fraction

public :: i2cf
public :: make_integer_continued_fraction
public :: make_integer_continued_fraction_from_integer

public :: constant_term_cf
public :: make_constant_term_continued_fraction
public :: make_constant_term_continued_fraction_from_integer

public :: apply_ng8
public :: apply_ng8_big_integers
public :: apply_ng8_integers
public :: ng8_coefficient_threshold
public :: ng8_term_threshold

public :: add_continued_fractions
public :: subtract_continued_fractions
public :: multiply_continued_fractions
public :: divide_continued_fractions

public :: cf2string
public :: continued_fraction_to_string_given_max_terms
public :: continued_fraction_to_string_with_default_max_terms
public :: default_continued_fraction_max_terms

type :: continued_fraction

class(continued_fraction_record), pointer, private :: p => null ()

contains

procedure, pass :: get_term => get_continued_fraction_term
procedure, pass :: term_exists => continued_fraction_term_exists
procedure, pass :: term => continued_fraction_term

procedure, pass :: to_string => continued_fraction_to_string_with_default_max_terms

procedure, pass :: add => add_continued_fractions
generic :: operator(+) => add

procedure, pass :: subtract => subtract_continued_fractions
generic :: operator(-) => subtract

procedure, pass :: multiply => multiply_continued_fractions
generic :: operator(*) => multiply

procedure, pass :: divide => divide_continued_fractions
generic :: operator(/) => divide

procedure, pass, private :: continued_fraction_make_new_ref
generic :: assignment(=) => continued_fraction_make_new_ref

final :: continued_fraction_final

end type continued_fraction

type :: continued_fraction_record
logical, private :: terminated = .false. ! No more terms?
integer, private :: m = 0 ! No. of terms memoized.
type(big_integer), private, allocatable :: memo(:) ! Memoized terms.
class(term_generator), pointer :: gen ! Where terms come from.
integer :: refcount = 0
contains
procedure, pass :: get_term => get_continued_fraction_record_term
procedure, pass :: term_exists => continued_fraction_record_term_exists
procedure, pass :: term => continued_fraction_record_term
final :: continued_fraction_record_final
end type continued_fraction_record

type, abstract :: term_generator
contains
procedure(term_generator_procedure), pass, deferred :: generate
end type term_generator

interface
subroutine term_generator_procedure (gen, term_exists, term)
import term_generator
import big_integer
class(term_generator), intent(inout) :: gen
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term
end subroutine term_generator_procedure
end interface

type, extends (term_generator) :: integer_term_generator
type(big_integer), allocatable :: term
logical :: no_more_terms = .false.
contains
procedure, pass :: generate => integer_term_generator_generate
end type integer_term_generator

type, extends (term_generator) :: constant_term_generator
type(big_integer), allocatable :: term
contains
procedure, pass :: generate => constant_term_generator_generate
end type constant_term_generator

type, extends (term_generator) :: ng8_term_generator
type(big_integer), allocatable :: a12, a1, a2, a
type(big_integer), allocatable :: b12, b1, b2, b
type(continued_fraction) :: x, y
integer :: ix = 0
integer :: iy = 0
logical :: x_overflow = .false.
logical :: y_overflow = .false.
contains
procedure, pass :: generate => ng8_term_generator_generate
end type ng8_term_generator

interface i2cf
module procedure make_integer_continued_fraction
module procedure make_integer_continued_fraction_from_integer
end interface i2cf

interface constant_term_cf
module procedure make_constant_term_continued_fraction
module procedure make_constant_term_continued_fraction_from_integer
end interface constant_term_cf

interface apply_ng8
module procedure apply_ng8_big_integers
module procedure apply_ng8_integers
end interface apply_ng8

interface cf2string
module procedure continued_fraction_to_string_given_max_terms
module procedure continued_fraction_to_string_with_default_max_terms
end interface cf2string

integer :: default_continued_fraction_max_terms = 20

type(big_integer), allocatable :: ng8_coefficient_threshold
type(big_integer), allocatable :: ng8_term_threshold

contains

subroutine continued_fraction_make_new_ref (dst, src)
class(continued_fraction), intent(inout) :: dst
class(continued_fraction), intent(in) :: src

if (associated (dst%p)) deallocate (dst%p)
dst%p => src%p
dst%p%refcount = dst%p%refcount + 1
end subroutine continued_fraction_make_new_ref

subroutine continued_fraction_final (cf)
type(continued_fraction), intent(inout) :: cf
cf%p%refcount = cf%p%refcount - 1
if (cf%p%refcount == 0) deallocate (cf%p)
end subroutine continued_fraction_final

function make_continued_fraction (gen) result (cf)
class(term_generator), pointer, intent(in) :: gen
type(continued_fraction) :: cf

allocate (cf%p)
allocate (cf%p%memo(0:31)) ! The starting size is arbitrary.
cf%p%gen => gen
cf%p%refcount = cf%p%refcount + 1
end function make_continued_fraction

subroutine continued_fraction_record_final (cfrec)
type(continued_fraction_record), intent(inout) :: cfrec
deallocate (cfrec%gen)
end subroutine continued_fraction_record_final

function make_integer_continued_fraction (bigint) result (cf)
type(big_integer), intent(in) :: bigint
type(continued_fraction) :: cf

class(integer_term_generator), pointer :: gen

allocate (gen)
gen%term = bigint
cf = make_continued_fraction (gen)
end function make_integer_continued_fraction

function make_integer_continued_fraction_from_integer (i) result (cf)
integer, intent(in) :: i
type(continued_fraction) :: cf
cf = make_integer_continued_fraction (integer2big (i))
end function make_integer_continued_fraction_from_integer

subroutine integer_term_generator_generate (gen, term_exists, term)
class(integer_term_generator), intent(inout) :: gen
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term

term_exists = (.not. gen%no_more_terms)
if (term_exists) term = gen%term
gen%no_more_terms = .true.
end subroutine integer_term_generator_generate

function make_constant_term_continued_fraction (bigint) result (cf)
type(big_integer), intent(in) :: bigint
type(continued_fraction) :: cf

class(constant_term_generator), pointer :: gen

allocate (gen)
gen%term = bigint
cf = make_continued_fraction (gen)
end function make_constant_term_continued_fraction

function make_constant_term_continued_fraction_from_integer (i) result (cf)
integer, intent(in) :: i
type(continued_fraction) :: cf
cf = make_constant_term_continued_fraction (integer2big (i))
end function make_constant_term_continued_fraction_from_integer

subroutine constant_term_generator_generate (gen, term_exists, term)
class(constant_term_generator), intent(inout) :: gen
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term

term_exists = .true.
if (term_exists) term = gen%term
end subroutine constant_term_generator_generate

function apply_ng8_big_integers (a12, a1, a2, a, &
& b12, b1, b2, b, x, y) result (cf)
type(big_integer), intent(in) :: a12, a1, a2, a
type(big_integer), intent(in) :: b12, b1, b2, b
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf

class(ng8_term_generator), pointer :: gen

allocate (gen)
gen%a12 = a12; gen%a1 = a1; gen%a2 = a2; gen%a = a
gen%b12 = b12; gen%b1 = b1; gen%b2 = b2; gen%b = b
gen%x = x
gen%y = y
cf = make_continued_fraction (gen)
end function apply_ng8_big_integers

function apply_ng8_integers (a12, a1, a2, a, &
& b12, b1, b2, b, x, y) result (cf)
integer, intent(in) :: a12, a1, a2, a
integer, intent(in) :: b12, b1, b2, b
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf

cf = apply_ng8_big_integers (integer2big (a12), &
& integer2big (a1), &
& integer2big (a2), &
& integer2big (a), &
& integer2big (b12), &
& integer2big (b1), &
& integer2big (b2), &
& integer2big (b), x, y)
end function apply_ng8_integers

subroutine ng8_term_generator_generate (gen, term_exists, term)
class(ng8_term_generator), intent(inout) :: gen
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term

logical :: done
logical :: b12z, b1z, b2z, bz
type(big_integer), allocatable :: q12, r12
type(big_integer), allocatable :: q1, r1
type(big_integer), allocatable :: q2, r2
type(big_integer), allocatable :: q, r
logical :: absorb_x, absorb_y, compare_fracs

done = .false.
do while (.not. done)
absorb_x = .false.
absorb_y = .false.
compare_fracs = .false.

b12z = (big_sgn (gen%b12) == 0)
b1z = (big_sgn (gen%b1) == 0)
b2z = (big_sgn (gen%b2) == 0)
bz = (big_sgn (gen%b) == 0)

if (b12z .and. b1z .and. b2z .and. bz) then
! There are no more terms.
term_exists = .false.
done = .true.
else if (b2z .and. bz) then
absorb_x = .true.
else if (b2z .or. bz) then
absorb_y = .true.
else if (b1z) then
absorb_x = .true.
else
call big_divrem (gen%a1, gen%b1, q1, r1)
call big_divrem (gen%a2, gen%b2, q2, r2)
call big_divrem (gen%a, gen%b, q, r)
if (.not. b12z) then
call big_divrem (gen%a12, gen%b12, q12, r12)
if (big_cmp (q, q1) /= 0) then
compare_fracs = .true.
else if (big_cmp (q, q2) /= 0) then
compare_fracs = .true.
else if (big_cmp (q, q12) /= 0) then
compare_fracs = .true.
else
call output_term
done = .true.
end if
end if
end if

if (compare_fracs) call compare_fractions (absorb_x, absorb_y)
if (absorb_x) call absorb_x_term
if (absorb_y) call absorb_y_term
end do

contains

subroutine output_term
gen%a12 = gen%b12; gen%a1 = gen%b1; gen%a2 = gen%b2; gen%a = gen%b
gen%b12 = r12; gen%b1 = r1; gen%b2 = r2; gen%b = r
term_exists = (.not. treat_as_infinite (q))
if (term_exists) term = q
end subroutine output_term

subroutine compare_fractions (absorb_x, absorb_y)
logical, intent(inout) :: absorb_x, absorb_y

! Rather than compare fractions, we will put the numerators over
! a common denominator of b1*b2*b, and then compare the new
! numerators.

type(big_integer), allocatable :: n1, n2, n

n1 = gen%a1 * gen%b2 * gen%b
n2 = gen%a2 * gen%b1 * gen%b
n = gen%a * gen%b1 * gen%b2
if (big_cmpabs (n1 - n, n2 - n) > 0) then
absorb_x = .true.
else
absorb_y = .true.
end if
end subroutine compare_fractions

subroutine absorb_x_term
logical :: term_exists
type(big_integer), allocatable :: term
type(big_integer), allocatable :: new_a12, new_a1, new_a2, new_a
type(big_integer), allocatable :: new_b12, new_b1, new_b2, new_b

if (gen%x_overflow) then
term_exists = .false.
else
term_exists = gen%x%term_exists(gen%ix)
end if
new_a2 = gen%a12
new_a = gen%a1
new_b2 = gen%b12
new_b = gen%b1
if (term_exists) then
term = gen%x%term(gen%ix)
new_a12 = gen%a2 + (gen%a12 * term)
new_a1 = gen%a + (gen%a1 * term)
new_b12 = gen%b2 + (gen%b12 * term)
new_b1 = gen%b + (gen%b1 * term)
if (any_too_big (new_a12, new_a1, new_a2, new_a, &
& new_b12, new_b1, new_b2, new_b)) then
gen%x_overflow = .true.
new_a12 = gen%a12
new_a1 = gen%a1
new_b12 = gen%b12
new_b1 = gen%b1
end if
else
new_a12 = gen%a12
new_a1 = gen%a1
new_b12 = gen%b12
new_b1 = gen%b1
end if
gen%a12 = new_a12; gen%a1 = new_a1; gen%a2 = new_a2; gen%a = new_a
gen%b12 = new_b12; gen%b1 = new_b1; gen%b2 = new_b2; gen%b = new_b
gen%ix = gen%ix + 1
end subroutine absorb_x_term

subroutine absorb_y_term
logical :: term_exists
type(big_integer), allocatable :: term
type(big_integer), allocatable :: new_a12, new_a1, new_a2, new_a
type(big_integer), allocatable :: new_b12, new_b1, new_b2, new_b

if (gen%y_overflow) then
term_exists = .false.
else
term_exists = gen%y%term_exists(gen%iy)
end if
new_a1 = gen%a12
new_a = gen%a2
new_b1 = gen%b12
new_b = gen%b2
if (term_exists) then
term = gen%y%term(gen%iy)
new_a12 = gen%a1 + (gen%a12 * term)
new_a2 = gen%a + (gen%a2 * term)
new_b12 = gen%b1 + (gen%b12 * term)
new_b2 = gen%b + (gen%b2 * term)
if (any_too_big (new_a12, new_a1, new_a2, new_a, &
& new_b12, new_b1, new_b2, new_b)) then
gen%y_overflow = .true.
new_a12 = gen%a12
new_a2 = gen%a2
new_b12 = gen%b12
new_b2 = gen%b2
end if
else
new_a12 = gen%a12
new_a2 = gen%a2
new_b12 = gen%b12
new_b2 = gen%b2
end if
gen%a12 = new_a12; gen%a1 = new_a1; gen%a2 = new_a2; gen%a = new_a
gen%b12 = new_b12; gen%b1 = new_b1; gen%b2 = new_b2; gen%b = new_b
gen%iy = gen%iy + 1
end subroutine absorb_y_term

function any_too_big (a, b, c, d, e, f, g, h) result (yes)
type(big_integer), intent(in) :: a, b, c, d, e, f, g, h
logical :: yes

if (too_big (a)) then
yes = .true.
else if (too_big (b)) then
yes = .true.
else if (too_big (c)) then
yes = .true.
else if (too_big (d)) then
yes = .true.
else if (too_big (e)) then
yes = .true.
else if (too_big (f)) then
yes = .true.
else if (too_big (g)) then
yes = .true.
else if (too_big (h)) then
yes = .true.
else
yes = .false.
end if
end function any_too_big

function too_big (coef) result (yes)
type(big_integer), intent(in) :: coef
logical :: yes

if (.not. allocated (ng8_coefficient_threshold)) then
! 2**512
ng8_coefficient_threshold = string2big ('1340780792994259709957&
&402499820584612747936582059239337772356144372176403007354&
&697680187429816690342769003185818648605085375388281194656&
&9946433649006084096')
end if

yes = (big_cmpabs (coef, ng8_coefficient_threshold) >= 0)
end function too_big

function treat_as_infinite (term) result (yes)
type(big_integer), intent(in) :: term
logical :: yes

if (.not. allocated (ng8_term_threshold)) then
! 2**64
ng8_term_threshold = string2big ('18446744073709551616')
end if

yes = (big_cmpabs (term, ng8_term_threshold) >= 0)
end function treat_as_infinite

end subroutine ng8_term_generator_generate

function add_continued_fractions (x, y) result (cf)
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf
cf = apply_ng8 (0, 1, 1, 0, 0, 0, 0, 1, x, y)
end function add_continued_fractions

function subtract_continued_fractions (x, y) result (cf)
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf
cf = apply_ng8 (0, 1, -1, 0, 0, 0, 0, 1, x, y)
end function subtract_continued_fractions

function multiply_continued_fractions (x, y) result (cf)
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf
cf = apply_ng8 (1, 0, 0, 0, 0, 0, 0, 1, x, y)
end function multiply_continued_fractions

function divide_continued_fractions (x, y) result (cf)
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf
cf = apply_ng8 (0, 1, 0, 0, 0, 0, 1, 0, x, y)
end function divide_continued_fractions

subroutine get_continued_fraction_term (cf, i, term_exists, term)
class(continued_fraction), intent(in) :: cf
integer, intent(in) :: i
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term

call get_continued_fraction_record_term (cf%p, i, term_exists, term)
end subroutine get_continued_fraction_term

subroutine get_continued_fraction_record_term (cfrec, i, term_exists, term)
class(continued_fraction_record), intent(inout) :: cfrec
integer, intent(in) :: i
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term

if (i < 0) error stop

call update (i + 1)
term_exists = (i < cfrec%m)
if (term_exists) term = cfrec%memo(i)

contains

subroutine update (needed)
integer :: needed
logical :: term_exists1
type(big_integer), allocatable :: term1

if (.not. cfrec%terminated .and. cfrec%m < needed) then
if (size (cfrec%memo) < needed) then
block ! Allocate more storage.
type(big_integer), allocatable :: memo1(:)
allocate (memo1(0 : (2 * needed) - 1))
memo1(0:cfrec%m - 1) = cfrec%memo(0:cfrec%m - 1)
call move_alloc (memo1, cfrec%memo)
end block
end if
do while (.not. cfrec%terminated .and. cfrec%m < needed)
call cfrec%gen%generate (term_exists1, term1)
if (term_exists1) then
cfrec%memo(cfrec%m) = term1
cfrec%m = cfrec%m + 1
else
cfrec%terminated = .true.
end if
end do
end if
end subroutine update

end subroutine get_continued_fraction_record_term

function continued_fraction_term_exists (cf, i) result (term_exists)
class(continued_fraction), intent(in) :: cf
integer, intent(in) :: i
logical :: term_exists
term_exists = continued_fraction_record_term_exists (cf%p, i)
end function continued_fraction_term_exists

function continued_fraction_record_term_exists (cfrec, i) result (term_exists)
class(continued_fraction_record), intent(inout) :: cfrec
integer, intent(in) :: i
logical :: term_exists

type(big_integer), allocatable :: term

call get_continued_fraction_record_term (cfrec, i, term_exists, term)
end function continued_fraction_record_term_exists

function continued_fraction_term (cf, i) result (term)
class(continued_fraction), intent(in) :: cf
integer, intent(in) :: i
type(big_integer), allocatable :: term
term = continued_fraction_record_term (cf%p, i)
end function continued_fraction_term

function continued_fraction_record_term (cfrec, i) result (term)
class(continued_fraction_record), intent(inout) :: cfrec
integer, intent(in) :: i
type(big_integer), allocatable :: term

logical :: term_exists

call get_continued_fraction_record_term (cfrec, i, term_exists, term)
if (.not. term_exists) error stop
end function continued_fraction_record_term

function continued_fraction_to_string_given_max_terms (cf, max_terms) result (s)
class(continued_fraction), intent(in) :: cf
integer, intent(in) :: max_terms
character(len = :), allocatable :: s
s = continued_fraction_record_to_string_given_max_terms (cf%p, max_terms)
end function continued_fraction_to_string_given_max_terms

function continued_fraction_record_to_string_given_max_terms (cfrec, max_terms) result (s)
class(continued_fraction_record), intent(inout) :: cfrec
integer, intent(in) :: max_terms
character(len = :), allocatable :: s

integer :: i
logical :: done

i = 0
s = '['
done = .false.
do while (.not. done)
if (.not. cfrec%term_exists(i)) then
s = s // "]"
done = .true.
else if (i == max_terms) then
s = s // ",...]"
done = .true.
else
select case (i)
case (0); continue
case (1); s = s // ";"
case default; s = s // ","
end select
s = s // big2string (cfrec%term(i))
i = i + 1
end if
end do
end function continued_fraction_record_to_string_given_max_terms

function continued_fraction_to_string_with_default_max_terms (cf) result (s)
class(continued_fraction), intent(in) :: cf
character(len = :), allocatable :: s
s = continued_fraction_record_to_string_with_default_max_terms (cf%p)
end function continued_fraction_to_string_with_default_max_terms

function continued_fraction_record_to_string_with_default_max_terms (cfrec) result (s)
class(continued_fraction_record), intent(inout) :: cfrec
character(len = :), allocatable :: s

integer :: max_terms

max_terms = max (default_continued_fraction_max_terms, 1)
s = continued_fraction_record_to_string_given_max_terms (cfrec, max_terms)
end function continued_fraction_record_to_string_with_default_max_terms

end module continued_fractions

!---------------------------------------------------------------------

program bivariate_continued_fraction_task

use, non_intrinsic :: big_integers
use, non_intrinsic :: continued_fractions
implicit none

type(continued_fraction) :: golden_ratio
type(continued_fraction) :: silver_ratio
type(continued_fraction) :: sqrt2
type(continued_fraction) :: one
type(continued_fraction) :: two
type(continued_fraction) :: three
type(continued_fraction) :: four
type(continued_fraction) :: method1
type(continued_fraction) :: method2
type(continued_fraction) :: method3

golden_ratio = constant_term_cf (1)
silver_ratio = constant_term_cf (2)
one = i2cf (1)
two = i2cf (2)
three = i2cf (3)
four = i2cf (4)
sqrt2 = silver_ratio - one

method1 = apply_ng8 (0, 1, 0, 0, 0, 0, 2, 0, silver_ratio, sqrt2)
method2 = apply_ng8 (1, 0, 0, 1, 0, 0, 0, 8, silver_ratio, silver_ratio)
method3 = (one / two) * (one + (one / sqrt2))

call show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2")
call show ("silver ratio", silver_ratio, "(1 + sqrt(2))")
call show ("sqrt(2)", sqrt2, "silver ratio minus 1")
call show ("13/11", i2cf (13) / i2cf (11), "")
call show ("22/7", i2cf (22) / i2cf (7), "")
call show ("1", one, "")
call show ("2", two, "")
call show ("3", three, "")
call show ("4", four, "")
call show ("(1 + 1/sqrt(2))/2", method1, "method 1")
call show ("(1 + 1/sqrt(2))/2", method2, "method 2")
call show ("(1 + 1/sqrt(2))/2", method3, "method 3")
call show ("sqrt(2) + sqrt(2)", sqrt2 + sqrt2, "")
call show ("sqrt(2) - sqrt(2)", sqrt2 - sqrt2, "")
call show ("sqrt(2) * sqrt(2)", sqrt2 * sqrt2, "")
call show ("sqrt(2) / sqrt(2)", sqrt2 / sqrt2, "")

contains

subroutine show (expression, cf, note)
character(len = *), intent(in) :: expression
class(continued_fraction), intent(in) :: cf
character(len = *), intent(in) :: note

write (*, '(A19, " => ", A, T73, A)') &
& expression, cf%to_string(), note
end subroutine show

end program bivariate_continued_fraction_task

!---------------------------------------------------------------------
</syntaxhighlight>

{{out}}
<pre>$ gfortran -O2 -g -fbounds-check -Wall -Wextra bivariate_continued_fraction_task.f90 && ./a.out
golden ratio => [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...] (1 + sqrt(5))/2
silver ratio => [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] (1 + sqrt(2))
sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] silver ratio minus 1
13/11 => [1;5,2]
22/7 => [3;7]
1 => [1]
2 => [2]
3 => [3]
4 => [4]
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 1
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 2
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 3
sqrt(2) + sqrt(2) => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt(2) - sqrt(2) => [0]
sqrt(2) * sqrt(2) => [2]
sqrt(2) / sqrt(2) => [1]
</pre>
</pre>


Retrieved from "https://rosettacode.org/wiki/Continued_fraction/Arithmetic/G(matrix_ng,_continued_fraction_n1,_continued_fraction_n2)"