! =============================================================================
! metallic.f90  --  Metallic Ratios via Lucas-like sequences
!
! The b-th Metallic Ratio is the positive root of  x^2 - bx - 1 = 0,
! evaluated by the closed form  (b + sqrt(b^2 + 4)) / 2.
!
! It can also be approximated by the ratio of successive terms of the
! Lucas-like recurrence  x(n) = b*x(n-1) + x(n-2),  started at x(1)=x(2)=1.
! As n grows the ratio x(n+1)/x(n) converges to the b-th Metallic Ratio.
!
! This program:
!   - Generates and displays the first 15 terms for b = 0 .. 9
!   - Iterates the recurrence until the ratio stabilises to 32 decimal places
!   - Reports the converged value and the iteration count
!   - Stretch goal: repeats for b=1 (Golden Ratio) to 256 decimal places
!
! Because standard floating-point types top out at ~18 significant decimal
! digits, all ratio arithmetic is done in a hand-rolled arbitrary-precision
! integer module (bigint_mod).  The ratio a/b is extracted one decimal digit
! at a time by the long-division routine bg_div_dec.
!
! Compile:  gfortran -O2 -o metallic metallic.f90
! =============================================================================

! =============================================================================
! MODULE bigint_mod
!
! Minimal arbitrary-precision non-negative integer library.
!
! Representation
! --------------
! A big integer is stored as an array of "limbs"  a(1..na)  in base BB=10^9.
!   a(1)  holds the LEAST-significant 9 decimal digits.
!   a(na) holds the MOST-significant (possibly shorter) group.
! Only limbs a(1)..a(na) are active; a(na+1)..a(ML) are ignored (but are
! kept zeroed by all routines to avoid stale data causing wrong results).
!
! Why base 10^9?
!   Each limb is an int64 (range up to ~9.2e18).  During the multiply-by-10
!   inner loop we temporarily hold  a(i)*10 + carry  which is at most
!   (10^9 - 1)*10 + (10^9 - 1) = ~10^10, well within int64.  Using a power
!   of ten as the base also makes printing decimal digits straightforward.
!
! Public interface
! ----------------
!   bg_set      -- initialise a big-int from a small (< BB) integer
!   bg_copy     -- copy one big-int to another
!   bg_add      -- c = a + b
!   bg_mul_s    -- c = a * k  (k is a small ordinary integer, 0 <= k <= ~100)
!   bg_sub      -- c = a - b  (a >= b must hold)
!   bg_cmp      -- compare two big-ints, return +1/0/-1
!   bg_div_dec  -- divide a/b and return the result as a decimal string
!                  to an arbitrary number of places
! =============================================================================
module bigint_mod
   use iso_fortran_env, only: int64
   implicit none

   ! ML  : maximum number of limbs per big-int.
   !       ML=200 supports numbers with up to 200*9 = 1800 decimal digits,
   !       which is vastly more than the ~300 digits needed here.
   integer, parameter :: ml = 200

   ! BB  : the limb base.  Every limb satisfies  0 <= limb < BB.
   integer(kind=int64), parameter :: bb = 1000000000_int64  ! 10^9

contains

   ! ---------------------------------------------------------------------------
   ! bg_set(a, na, val)
   !
   ! Set the big-int  a  to the small non-negative value  val  (val < BB).
   ! All limbs are zeroed first to guarantee a clean state; then limb 1 is set
   ! to val and na is set to 1 (a single-limb number).
   ! ---------------------------------------------------------------------------
   subroutine bg_set(a, na, val)
      integer(kind=int64), intent(out) :: a(ml)   ! limb array to initialise
      integer,             intent(out) :: na       ! number of active limbs
      integer(kind=int64), intent(in)  :: val      ! value to store (< BB)

      a = 0_int64   ! zero every limb so no stale data remains
      a(1) = val    ! least-significant limb carries the whole value
      na = 1        ! exactly one active limb
   end subroutine bg_set

   ! ---------------------------------------------------------------------------
   ! bg_copy(src, ns, dst, nd)
   !
   ! Copy big-int src (with ns active limbs) into dst, setting nd = ns.
   ! A whole-array copy is used because ML is fixed at compile time; it is
   ! fast and ensures no stale limbs survive in dst beyond position ns.
   ! ---------------------------------------------------------------------------
   subroutine bg_copy(src, ns, dst, nd)
      integer(kind=int64), intent(in)  :: src(ml)  ! source big-int
      integer,             intent(in)  :: ns        ! active limbs in src
      integer(kind=int64), intent(out) :: dst(ml)  ! destination big-int
      integer,             intent(out) :: nd        ! active limbs in dst

      dst = src   ! copies all ML limbs (unused ones are already zero)
      nd  = ns
   end subroutine bg_copy

   ! ---------------------------------------------------------------------------
   ! bg_add(a, na, b, nb, c, nc)
   !
   ! c = a + b
   !
   ! Standard schoolbook addition, limb by limb from least significant to most.
   ! At each position we accumulate  carry + a(i) + b(i),  write the result mod
   ! BB into c(i), and pass floor(result/BB) as carry to the next position.
   ! We iterate up to max(na,nb)+1 limbs to catch any final carry overflow.
   ! After the loop the leading-zero trimmer keeps nc as small as possible.
   !
   ! NOTE: c must be a distinct array from both a and b (no aliasing).
   ! ---------------------------------------------------------------------------
   subroutine bg_add(a, na, b, nb, c, nc)
      integer(kind=int64), intent(in)  :: a(ml), b(ml)  ! operands
      integer,             intent(in)  :: na, nb          ! their active limbs
      integer(kind=int64), intent(out) :: c(ml)          ! result  c = a + b
      integer,             intent(out) :: nc              ! active limbs in c

      integer(kind=int64) :: carry   ! carry propagating upward through limbs
      integer :: i, nm

      nm = max(na, nb)   ! highest limb index that might be non-zero in a or b
      c = 0_int64
      carry = 0_int64

      do i = 1, nm + 1   ! +1 to absorb a possible carry out of the top limb
         c(i) = carry
         if (i <= na) c(i) = c(i) + a(i)   ! add a's contribution (if present)
         if (i <= nb) c(i) = c(i) + b(i)   ! add b's contribution (if present)
         carry   = c(i) / bb                ! extract carry for next limb
         c(i)    = mod(c(i), bb)            ! keep only the low BB part
      end do

      ! Set nc to the index of the highest non-zero limb (at least 1).
      nc = nm + 1
      do while (nc > 1 .and. c(nc) == 0_int64)
         nc = nc - 1
      end do
   end subroutine bg_add

   ! ---------------------------------------------------------------------------
   ! bg_mul_s(a, na, k, c, nc)
   !
   ! c = a * k,  where k is a small non-negative integer (0 <= k <= ~100).
   !
   ! Multiplies each limb by k, accumulating carries upward.  The product
   !   a(i)*k + carry_in
   ! fits comfortably in int64 because:
   !   a(i) < 10^9,  k <= 100,  carry_in < 100
   !   => max product ~= 10^9 * 100 + 100 = 10^11,  well within 2^63 - 1.
   !
   ! k=0 is handled specially (result is zero stored as a single zero limb).
   !
   ! NOTE: c must be a distinct array from a (no aliasing).
   ! ---------------------------------------------------------------------------
   subroutine bg_mul_s(a, na, k, c, nc)
      integer(kind=int64), intent(in)  :: a(ml)   ! multiplicand
      integer,             intent(in)  :: na, k   ! active limbs; small factor
      integer(kind=int64), intent(out) :: c(ml)   ! result  c = a * k
      integer,             intent(out) :: nc       ! active limbs in c

      integer(kind=int64) :: carry, kk   ! carry between limbs; k as int64
      integer :: i

      c = 0_int64

      ! Special case: multiplying by zero yields zero (single zero limb).
      if (k == 0) then
         nc = 1
         return
      end if

      kk    = int(k, int64)   ! widen k to int64 to avoid 32-bit overflow
      carry = 0_int64

      do i = 1, na
         c(i)  = a(i) * kk + carry   ! partial product plus incoming carry
         carry = c(i) / bb            ! carry for the next limb
         c(i)  = mod(c(i), bb)        ! store only the low BB part
      end do

      ! If there is a leftover carry it forms a new most-significant limb.
      if (carry > 0_int64) then
         nc      = na + 1
         c(nc)   = carry
      else
         nc = na
      end if
   end subroutine bg_mul_s

   ! ---------------------------------------------------------------------------
   ! bg_sub(a, na, b, nb, c, nc)
   !
   ! c = a - b   (caller must ensure a >= b; the result is always >= 0)
   !
   ! Standard schoolbook subtraction with borrow.  We work on a local copy
   ! `tmp` so that c may safely alias a (c and a may point to the same array).
   ! After subtraction the leading-zero trimmer reduces nc.
   ! ---------------------------------------------------------------------------
   subroutine bg_sub(a, na, b, nb, c, nc)
      integer(kind=int64), intent(in)  :: a(ml), b(ml)  ! a >= b required
      integer,             intent(in)  :: na, nb
      integer(kind=int64), intent(out) :: c(ml)         ! c = a - b
      integer,             intent(out) :: nc

      integer(kind=int64) :: borrow    ! borrow propagating upward (0 or 1)
      integer(kind=int64) :: tmp(ml)   ! local working copy so c may alias a
      integer :: i

      tmp    = a            ! copy a into tmp so we can safely write into c
      borrow = 0_int64

      do i = 1, na
         tmp(i) = tmp(i) - borrow            ! subtract any borrow from above
         if (i <= nb) tmp(i) = tmp(i) - b(i) ! subtract b's limb (if present)
         ! If the result went negative, borrow from the next higher limb.
         if (tmp(i) < 0_int64) then
            tmp(i) = tmp(i) + bb   ! add one full base to restore positivity
            borrow = 1_int64       ! and carry that borrow upward
         else
            borrow = 0_int64       ! result non-negative: no borrow needed
         end if
      end do

      ! Trim leading zero limbs (but always keep at least one limb).
      nc = na
      do while (nc > 1 .and. tmp(nc) == 0_int64)
         nc = nc - 1
      end do
      c = tmp   ! write result out (safe even if c aliases a, since tmp is local)
   end subroutine bg_sub

   ! ---------------------------------------------------------------------------
   ! bg_cmp(a, na, b, nb)  ->  integer result
   !
   ! Compare two big-ints.  Returns:
   !   +1  if  a > b
   !    0  if  a = b
   !   -1  if  a < b
   !
   ! Strategy:
   !   1. A number with more active limbs is larger (they share the same base).
   !   2. If both have the same limb count, scan from most significant to least;
   !      the first differing limb determines the outcome.
   !   3. If every limb matches the numbers are equal.
   ! ---------------------------------------------------------------------------
   function bg_cmp(a, na, b, nb) result(res)
      integer(kind=int64), intent(in) :: a(ml), b(ml)
      integer,             intent(in) :: na, nb
      integer :: res, i

      ! More limbs => larger magnitude (no leading-zero limbs are ever stored).
      if (na > nb) then;  res =  1;  return;  end if
      if (na < nb) then;  res = -1;  return;  end if

      ! Same limb count: scan from the most significant limb downward.
      do i = na, 1, -1
         if (a(i) > b(i)) then;  res =  1;  return;  end if
         if (a(i) < b(i)) then;  res = -1;  return;  end if
      end do

      res = 0   ! every limb matched: the numbers are equal
   end function bg_cmp

   ! ---------------------------------------------------------------------------
   ! bg_div_dec(a, na, b, nb, ndec, rstr)
   !
   ! Compute  a / b  and write the result into the character variable  rstr
   ! in decimal notation with exactly  ndec  digits after the decimal point.
   ! Example:  a/b = golden ratio, ndec=32  =>  rstr = "1.61803398874989484820..."
   !
   ! The caller must ensure  a >= b  (i.e. ratio >= 1), and that  rstr  is
   ! at least  ndec + len(int_part) + 1  characters wide.
   !
   ! Algorithm
   ! ---------
   ! Step 1 -- Integer part
   !   Find the largest integer  q  such that  q*b <= a,  by incrementing q
   !   from 0 while  (q+1)*b <= a.  For our sequences q is at most b+1 (the
   !   metallic parameter), so at most ~11 multiplications are needed.
   !   Compute the remainder  r = a - q*b.
   !
   ! Step 2 -- Decimal digits (long division)
   !   The classical long-division digit loop:
   !     r  <-  r * 10                    (shift remainder left one decimal place)
   !     d  <-  floor(r / b),  0 <= d <= 9  (next decimal digit)
   !     r  <-  r - d*b                   (new remainder for the next digit)
   !   Finding d uses the same trial-multiplication trick: start d=0, increment
   !   while (d+1)*b <= r.  Since 0 <= d <= 9 this costs at most 10 multiplies.
   !
   ! Step 3 -- Assemble result string
   !   Write the integer part using Fortran's I0 format, append '.', then
   !   append the decimal digits one character at a time.
   ! ---------------------------------------------------------------------------
   subroutine bg_div_dec(a, na, b, nb, ndec, rstr)
      integer(kind=int64), intent(in)  :: a(ml), b(ml)  ! numerator, denominator
      integer,             intent(in)  :: na, nb         ! their active limbs
      integer,             intent(in)  :: ndec           ! decimal places wanted
      character(len=*),    intent(out) :: rstr           ! result string

      integer(kind=int64) :: r(ml), t1(ml), t2(ml)  ! remainder; temporaries
      integer :: nr, nt1, nt2                        ! their active limb counts
      integer :: int_part   ! integer part of a/b (small, ~0-10 for our use)
      integer :: d          ! one decimal digit (0-9)
      integer :: j          ! loop counter over decimal positions
      integer :: pos        ! current write position in rstr
      character(len=20) :: ips   ! scratch to format the integer part

      ! ------------------------------------------------------------------
      ! Step 1: find the integer part q = floor(a/b).
      ! We compute (q+1)*b and check whether it still <= a.  If so, q is
      ! not yet large enough and we increment.  When (q+1)*b > a we stop
      ! and q is correct.
      ! ------------------------------------------------------------------
      int_part = 0
      call bg_mul_s(b, nb, int_part + 1, t1, nt1)   ! t1 = 1*b
      do while (bg_cmp(t1, nt1, a, na) <= 0)         ! while (int_part+1)*b <= a
         int_part = int_part + 1
         call bg_mul_s(b, nb, int_part + 1, t1, nt1) ! next candidate
      end do

      ! Compute the remainder  r = a - int_part * b
      call bg_mul_s(b, nb, int_part, t1, nt1)   ! t1 = int_part * b
      call bg_sub(a, na, t1, nt1, r, nr)         ! r  = a - t1

      ! ------------------------------------------------------------------
      ! Assemble the start of the result string: "<integer_part>."
      ! I0 format suppresses leading zeros, e.g. 1 -> "1", 12 -> "12".
      ! After assignment to rstr (a fixed-length char variable) Fortran
      ! pads with spaces to the right; the subsequent character-by-character
      ! writes will overwrite those spaces with decimal digits.
      ! pos points to the first character position after the decimal point.
      ! ------------------------------------------------------------------
      write(ips, '(I0)') int_part
      rstr = trim(ips) // '.'
      pos  = len_trim(ips) + 2   ! 1-based: skip int_part digits and the '.'

      ! ------------------------------------------------------------------
      ! Step 2: long-division decimal digit loop.
      ! Each iteration extracts exactly one decimal digit and updates r.
      ! ------------------------------------------------------------------
      do j = 1, ndec

         ! Shift the remainder left by one decimal place (multiply by 10).
         ! We must use a temporary because bg_mul_s cannot alias input/output.
         call bg_mul_s(r, nr, 10, t1, nt1)
         r  = t1    ! update r in-place via copy (t1 is a local array)
         nr = nt1

         ! Find the digit d = floor(r / b),  0 <= d <= 9.
         ! We try d=0,1,2,... until (d+1)*b > r.
         d = 0
         call bg_mul_s(b, nb, d + 1, t1, nt1)   ! t1 = (d+1)*b = 1*b
         do while (bg_cmp(t1, nt1, r, nr) <= 0)  ! while (d+1)*b <= r
            d = d + 1
            call bg_mul_s(b, nb, d + 1, t1, nt1) ! next candidate digit
         end do
         ! Now d is the correct decimal digit.

         ! Store the digit as a character ('0'+d).
         rstr(pos:pos) = char(ichar('0') + d)
         pos = pos + 1

         ! Update the remainder:  r <- r - d*b
         call bg_mul_s(b, nb, d, t1, nt1)        ! t1 = d*b
         call bg_sub(r, nr, t1, nt1, t2, nt2)    ! t2 = r - t1
         r  = t2    ! copy back into r
         nr = nt2

      end do   ! next decimal digit

   end subroutine bg_div_dec

end module bigint_mod


! =============================================================================
! PROGRAM metallic_ratios
!
! Main driver.  For each of the first ten Metallic Ratios (b = 0 .. 9):
!   1. Generate and print the first 15 terms of the Lucas-like sequence
!        x(n) = b * x(n-1) + x(n-2),   x(1) = x(2) = 1.
!   2. Iterate the recurrence using arbitrary-precision big-ints, computing
!      the ratio x(n+1)/x(n) to 32 decimal places after each step.
!   3. Stop when two consecutive ratio strings are identical (convergence).
!   4. Report the converged ratio and the iteration index n.
!
! Stretch goal: repeat step 2-4 for b=1 (Golden Ratio) to 256 decimal places.
!
! Variable layout
! ---------------
!   xa, xb   -- the two most-recent sequence terms (big-ints)
!   xc       -- temporary holding  b * x(n),  the first part of the recurrence
!   xt       -- the new term  x(n+1) = b*x(n) + x(n-1)
!   na..nt_v -- active limb counts for xa..xt respectively
!   seq(15)  -- first 15 terms as ordinary int64 (they fit; largest is ~3e12)
!   prev_str -- ratio string from the previous iteration
!   curr_str -- ratio string from the current iteration
! =============================================================================
program metallic_ratios
   use bigint_mod
   use iso_fortran_env, only: int64
   implicit none

   ! Big-int working arrays for the two live sequence terms and two temporaries.
   integer(kind=int64) :: xa(ml), xb(ml), xc(ml), xt(ml)
   integer :: na    ! active limbs in xa  (the x(n-1) term)
   integer :: nb_v  ! active limbs in xb  (the x(n)   term)
   integer :: nc    ! active limbs in xc  (intermediate: b * x(n))
   integer :: nt_v  ! active limbs in xt  (the new x(n+1) term)

   integer :: bv        ! metallic parameter (the "b" in the task description)
   integer :: i         ! loop index for building seq()
   integer :: iter      ! iteration count (= n in the task description)
   integer :: ndec      ! how many decimal places we want

   ! First 15 terms stored as 64-bit integers (the largest, for b=9, is ~3e12,
   ! well within int64's range of ~9.2e18).
   integer(kind=int64) :: seq(15)

   ! Ratio strings.  Length 310 comfortably holds:
   !   2 chars  (integer part + '.') + 256 decimal digits + padding spaces
   character(len=310) :: prev_str   ! ratio string from the previous iteration
   character(len=310) :: curr_str   ! ratio string from the current iteration

   logical :: converged   ! .TRUE. once curr_str matches prev_str

   ! Human-readable names indexed by b.
   character(len=10) :: names(0:9)
   names(0) = 'Platinum'
   names(1) = 'Golden'
   names(2) = 'Silver'
   names(3) = 'Bronze'
   names(4) = 'Copper'
   names(5) = 'Nickel'
   names(6) = 'Aluminum'
   names(7) = 'Iron'
   names(8) = 'Tin'
   names(9) = 'Lead'

   ! -------------------------------------------------------------------------
   ! Header
   ! -------------------------------------------------------------------------
   write(*, '(A)') "Metallic Ratios via Lucas-like Sequences"
   write(*, '(A)') repeat('=', 62)
   write(*, *)

   ndec = 32   ! target precision for the main loop

   ! =========================================================================
   ! Main loop: b = 0, 1, ..., 9
   ! =========================================================================
   do bv = 0, 9

      ! -----------------------------------------------------------------------
      ! Part A: Generate the first 15 sequence terms using ordinary int64.
      !
      ! The sequence is  x(n) = bv * x(n-1) + x(n-2),  x(1) = x(2) = 1.
      ! For bv=0 every term is 1.
      ! For bv=1 we get the Fibonacci sequence: 1,1,2,3,5,8,...
      ! For bv=9 the 15th term is ~3.26e12, comfortably within int64.
      ! -----------------------------------------------------------------------
      seq(1) = 1_int64
      seq(2) = 1_int64
      do i = 3, 15
         seq(i) = int(bv, int64) * seq(i - 1) + seq(i - 2)
      end do

      write(*, '(A,I1,A,A)') "b=", bv, "  ", trim(names(bv))
      write(*, '(A)') "First 15 Lucas sequence terms:"
      write(*, '(5(1X,I18))') seq(1:5)    ! rows of 5 for compact display
      write(*, '(5(1X,I18))') seq(6:10)
      write(*, '(5(1X,I18))') seq(11:15)

      ! -----------------------------------------------------------------------
      ! Part B: Converge the ratio to ndec=32 decimal places.
      !
      ! Special case b=0: the sequence is identically 1, the ratio is
      ! exactly 1, and the formula (0 + sqrt(4)) / 2 = 1 confirms this.
      ! No iteration is needed.
      ! -----------------------------------------------------------------------
      if (bv == 0) then
         write(*, '(A)') "Ratio = 1  (exact, degenerate Platinum case)"
         write(*, *)
         cycle   ! skip to next b
      end if

      ! Seed the big-int sequence with x(1) = 1, x(2) = 1.
      call bg_set(xa, na,   1_int64)   ! xa = x(n-1) = 1  (previous-previous term)
      call bg_set(xb, nb_v, 1_int64)  ! xb = x(n)   = 1  (previous term)

      prev_str  = ''       ! empty string: guaranteed not to match first ratio
      converged = .false.
      iter      = 2        ! we already have terms 1 and 2 in xa and xb

      ! -----------------------------------------------------------------------
      ! Iteration loop.
      !
      ! Each pass:
      !   1. Compute  x(n+1) = bv * x(n) + x(n-1)   using big-int arithmetic.
      !   2. Compute the ratio  x(n+1) / x(n)  to ndec decimal places.
      !   3. Compare with the ratio from the previous iteration.
      !      If identical, the approximation has stabilised and we stop.
      !      Otherwise record the new ratio and advance the sequence.
      !
      ! Note on aliasing: bg_mul_s and bg_add require distinct arrays for
      ! input and output, so we use xc and xt as intermediates rather than
      ! overwriting xb in place.
      ! -----------------------------------------------------------------------
      do while (.not. converged)

         ! Step 1: compute the new sequence term.
         !   xc = bv * xb             (b * x(n))
         !   xt = xc + xa             (b * x(n) + x(n-1) = x(n+1))
         call bg_mul_s(xb, nb_v, bv, xc, nc)     ! xc = bv * x(n)
         call bg_add(xc, nc, xa, na, xt, nt_v)   ! xt = xc + x(n-1)

         ! Step 2: form the ratio x(n+1) / x(n) as a decimal string.
         call bg_div_dec(xt, nt_v, xb, nb_v, ndec, curr_str)

         ! Step 3: test for convergence.
         ! Fortran character comparison checks all len(curr_str) characters,
         ! including the trailing spaces that Fortran appends to fixed-length
         ! strings.  Both strings are built identically by bg_div_dec, so
         ! trailing-space counts match and the comparison is correct.
         if (curr_str == prev_str) then
            converged = .true.   ! ratio unchanged to ndec places -- done
         else
            ! Not yet converged: record this ratio and advance the sequence.
            !   xa <-- old xb  (old x(n) becomes the new x(n-1))
            !   xb <-- xt      (new x(n+1) becomes the new x(n))
            prev_str = curr_str
            call bg_copy(xb, nb_v, xa, na)    ! xa = old xb
            call bg_copy(xt, nt_v, xb, nb_v) ! xb = xt
            iter = iter + 1
         end if

      end do   ! while not converged

      write(*, '(A,I6)') "Converged at n =", iter
      write(*, '(2A)')   "Ratio = ", trim(curr_str)
      write(*, *)

   end do   ! bv = 0 .. 9

   ! ==========================================================================
   ! Stretch goal: Golden Ratio (b=1) to 256 decimal places.
   !
   ! The Golden Ratio has the SLOWEST convergence of all Metallic Ratios
   ! because its continued fraction [1;1,1,1,...] converges most slowly
   ! among all irrationals (a classical result in Diophantine approximation).
   ! Roughly 256 / log10(phi) ≈ 1225 iterations are needed.
   ! At that point each Fibonacci number has about 256 decimal digits, so
   ! the big-int arithmetic is the only practical approach.
   ! ==========================================================================
   write(*, '(A)') repeat('=', 62)
   write(*, '(A)') "STRETCH GOAL: Golden Ratio to 256 decimal places"
   write(*, '(A)') repeat('=', 62)
   write(*, *)

   ndec = 256
   bv   = 1      ! Golden Ratio: b=1

   call bg_set(xa, na,   1_int64)   ! seed x(1) = 1
   call bg_set(xb, nb_v, 1_int64)  ! seed x(2) = 1
   prev_str  = ''
   converged = .false.
   iter      = 2

   ! Same iteration loop as above; only ndec differs.
   do while (.not. converged)
      call bg_mul_s(xb, nb_v, bv, xc, nc)
      call bg_add(xc, nc, xa, na, xt, nt_v)
      call bg_div_dec(xt, nt_v, xb, nb_v, ndec, curr_str)
      if (curr_str == prev_str) then
         converged = .true.
      else
         prev_str = curr_str
         call bg_copy(xb, nb_v, xa, na)
         call bg_copy(xt, nt_v, xb, nb_v)
         iter = iter + 1
      end if
   end do

   write(*, '(A,I6)') "Converged at n =", iter
   write(*, '(A)')    "Golden Ratio ="

   ! Print the 258-character result ("1." + 256 digits) as 40-character lines
   ! so it is readable in a standard 80-column terminal.
   ! The slice boundaries step by 40 through the string, with the last slice
   ! picking up the remaining 18 characters (positions 241-258).
   write(*, '(A)') curr_str(  1: 42)
   write(*, '(A)') curr_str( 43: 82)
   write(*, '(A)') curr_str( 83:122)
   write(*, '(A)') curr_str(123:162)
   write(*, '(A)') curr_str(163:202)
   write(*, '(A)') curr_str(203:242)
   write(*, '(A)') curr_str(243:258)
   write(*, *)

end program metallic_ratios