Jump to content

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

Continued fraction/Arithmetic/G(matrix ng, continued fraction n1, continued fraction n2) (view source)

Revision as of 18:38, 14 March 2023

49,907 bytes added ,  1 year ago
→‎{{header|D}}
Line 2,626:
sqrt(2) * sqrt(2) => [2]
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>
 
Chemoelectric
1,448

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/338742"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information