Continued fraction/Arithmetic/G(matrix ng, continued fraction n1, continued fraction n2): Difference between revisions
Content added Content deleted
Line 3,254: | Line 3,254: | ||
! leading zero-bytes and is at least two bytes long. If b is one |
! leading zero-bytes and is at least two bytes long. If b is one |
||
! byte long and nonzero, use short division. |
! 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. |
|||
! |
! |
||
Line 3,266: | Line 3,260: | ||
integer :: j |
integer :: j |
||
integer :: qhat |
integer :: qhat |
||
logical :: carry |
|||
! |
! |
||
Line 3,311: | Line 3,305: | ||
call multiply_and_subtract (carry) |
call multiply_and_subtract (carry) |
||
q(j) = i2bk (qhat) |
q(j) = i2bk (qhat) |
||
if (carry |
if (carry) call add_back |
||
end do |
end do |
||
Line 3,432: | Line 3,426: | ||
subroutine multiply_and_subtract (carry) |
subroutine multiply_and_subtract (carry) |
||
logical, intent(out) :: carry |
|||
integer :: i |
integer :: i |
||
Line 3,453: | Line 3,447: | ||
u(j + i) = i2bk (diff) |
u(j + i) = i2bk (diff) |
||
end do |
end do |
||
carry = sub_carry |
carry = (sub_carry /= 0) |
||
end subroutine multiply_and_subtract |
end subroutine multiply_and_subtract |
||
Line 3,462: | Line 3,456: | ||
carry = 0 |
carry = 0 |
||
do i = 0, n - 1 |
do i = 0, n - 1 |
||
sum = bk2i (u(i)) + bk2i (v(i)) + carry |
sum = bk2i (u(j + i)) + bk2i (v(i)) + carry |
||
carry = ishft (sum, -8) |
carry = ishft (sum, -8) |
||
u(i) = i2bk (sum) |
u(j + i) = i2bk (sum) |
||
end do |
end do |
||
end subroutine add_back |
end subroutine add_back |
||
Line 3,535: | Line 3,529: | ||
! Shorten to the minimum number of bytes. |
! Shorten to the minimum number of bytes. |
||
i = size (a%bytes) |
i = size (a%bytes) |
||
if (1 < i) then |
|||
i |
do while (1 < i .and. a%bytes(i) == zero) |
||
i = i - 1 |
|||
end do |
|||
end if |
|||
if (i /= size (a%bytes)) then |
if (i /= size (a%bytes)) then |
||
allocate (fewer_bytes (i)) |
allocate (fewer_bytes (i)) |
||
Line 4,224: | Line 4,220: | ||
type(continued_fraction) :: method2 |
type(continued_fraction) :: method2 |
||
type(continued_fraction) :: method3 |
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) |
golden_ratio = constant_term_cf (1) |