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 23:18, 15 March 2023
, 1 year ago→{{header|Fortran}}
Line 3,254:
! leading zero-bytes and is at least two bytes long. If b is one
! byte long and nonzero, use short division.
!
Line 3,266 ⟶ 3,260:
integer :: j
integer :: qhat
!
Line 3,311 ⟶ 3,305:
call multiply_and_subtract (carry)
q(j) = i2bk (qhat)
if (carry
end do
Line 3,432 ⟶ 3,426:
subroutine multiply_and_subtract (carry)
integer :: i
Line 3,453 ⟶ 3,447:
u(j + i) = i2bk (diff)
end do
carry = (sub_carry /= 0)
end subroutine multiply_and_subtract
Line 3,462 ⟶ 3,456:
carry = 0
do i = 0, n - 1
sum = bk2i (u(j + i)) + bk2i (v(i)) + carry
carry = ishft (sum, -8)
u(j + i) = i2bk (sum)
end do
end subroutine add_back
Line 3,535 ⟶ 3,529:
! Shorten to the minimum number of bytes.
i = size (a%bytes)
do while (1 < i
end do
end if
if (i /= size (a%bytes)) then
allocate (fewer_bytes (i))
Line 4,224 ⟶ 4,220:
type(continued_fraction) :: method2
type(continued_fraction) :: method3
block
! Let us start with a test of the long division routine, with some
! numbers known to trigger a bug in the first version I
! posted. That version had a buggy "add_back" routine.
type(big_integer), allocatable :: a, b, q, r
a = string2big ("95292350834616415707142739736356097545484658215015733475&
&690528634954280668802285176284181482202905629004984123564335019024318905")
b = string2big ("63677949970178275389170357551071104191609722674550547056511830750")
call big_divrem (a, b, q, r)
if (big_sgn (((b * q) + r) - a) /= 0) error stop
a = string2big ("5286200801181288750950358142425694618335361315503743069535407838&
&1104411448878793976933118436177295215313131557463887718741957154")
b = string2big ("54401097470188014066128968444633185848791550678521")
call big_divrem (a, b, q, r)
if (big_sgn (((b * q) + r) - a) /= 0) error stop
a = string2big ("74352827755975214935544861176217106447734695144315262422&
&288346418457330023596489437154599028318030933202302606937951415862696330")
b = string2big ("291979433784649910583546698460221489986784915256036707914578&
&892106828527219639012712")
call big_divrem (a, b, q, r)
if (big_sgn (((b * q) + r) - a) /= 0) error stop
a = string2big ("7515839498480018152556264500298705705667515770181724145893&
&9544448656273749453586884931339958104923411455488844854494605760712247")
b = string2big ("8600698996698965932302079501896131441135807557744554902970070&
&402964318496325075877027770784963001")
call big_divrem (a, b, q, r)
if (big_sgn (((b * q) + r) - a) /= 0) error stop
a = string2big ("13370595927769020368832742717678609885835798503146654175875&
&149837801359758893206045930442389897206420586502996797614097489470778")
b = string2big ("871343613388")
call big_divrem (a, b, q, r)
if (big_sgn (((b * q) + r) - a) /= 0) error stop
end block
golden_ratio = constant_term_cf (1)
| |||