C Determine coefficients k of linear combinations
C of the square roots of 3 integers i, j, k > 1,
C such that the sum ki*sqrt(i) + kj*sqrt(j) + kk*sqrt(k)
C is minimized for a growing maximum m of the absolute
C value of the coefficients.
C The 3 integers are read from the command line.
C An optional 4th parameter is a search limit mmax.
C Without this parameter the default (and maximum) 10^6 is used.
C
C To produce A379918: Call with parameters 3 5 7
C A379919: Call with 2 3 5
C
C Tested with gfortran, other Fortan compilers not available.
C
C Author: Hugo Pfoertner http://www.pfoertner.org
C
C Version history
C 2025-01-10 Checks for the command parameters
C 2025-01-09 Initial version
C
C Implicit type declarations are assumed for integers.
C (i-n).
C
      parameter (mmax=1000000)
      implicit real*16 (q)
      doubleprecision s, s1(mmax), s2(mmax), s3(mmax),
     &   sr1, sr2, sr3,
     &   zero, one, vlarge,
     &   s1i, s2i, s3i, d, dd, ddd, eps
      character cline*10
C
      parameter (zero=0.0D0, one=1.0D0, vlarge=1.0D100)
C limit for switching to quadruple precision recomputation
C of sum.
      parameter (eps=1.0D-8)
C
C Parse command parameters
      nparm = command_argument_count ()
      if ( nparm .lt. 3 )
     & stop 'Please provide 3 distinct positive integers'
      call get_command_argument (1, cline)
      read (cline,*) ii1
      call get_command_argument (2, cline)
      read (cline,*) ii2
      call get_command_argument (3, cline)
      read (cline,*) ii3
      if (nparm .eq. 4) then
        call get_command_argument (4, cline)
        read (cline,*) mlimit
        mlimit = max(10000,min(mmax,mlimit))
      else
        mlimit = mmax
      endif
      ii1 = abs(ii1)
      ii2 = abs(ii2)
      ii3 = abs(ii3)
      isum = ii1 + ii2 + ii3
      i1 = min (ii1, ii2, ii3)
      i3 = max (ii1, ii2, ii3)
      i2 = isum - i1 - i3
      if ( i1 .eq. 0 .or. i1 .eq. i2 .or. i2 .eq. i3)
     & stop 'Parameters must be distinct and != 0'
C
C Quadruple precision square roots for sums < eps
      q1 = real(i1,kind=16)
      q2 = real(i2,kind=16)
      q3 = real(i3,kind=16)
      qsr1 = sqrt(q1)
      qsr2 = sqrt(q2)
      qsr3 = sqrt(q3)
C default double precision square roots
      sr1 = sqrt(dble(i1))
      sr2 = sqrt(dble(i2))
      sr3 = sqrt(dble(i3))
C
C Auxiliary arrays of constant multiples
      do 10 k = 1, mmax
      s1(k) = k * sr1
      s2(k) = k * sr2
      s3(k) = k * sr3
10    continue
C
C Inverses to avoid divisions
      s1i = one / sr1
      s2i = one / sr2
      s3i = one / sr3
C
      d = vlarge
      locf = 0
C
C Main loop setting one of the weights to a given maximum m
C
      do 100 m = 1, mlimit
C Print to show progress
      if (mod(m,10000) .eq. 0) write(*,*) -m
      s = zero
      dd = vlarge
C
C separate loops keeping one of the coefficients at m;
C 3rd coefficient must not exceed m
      do 12 j = 1, m
      s = s1(m) + s2(j)
      k = nint(s*s3i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr3)
        if (ddd .lt. dd) then
          kk1 = m
          kk2 = j
          kk3 = -k
          dd = ddd
          locf = 121
        endif
      endif
      s = s1(m) - s2(j)
      k = nint(s*s3i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr3)
        if (ddd .lt. dd) then
          kk1 = m
          kk2 = -j
          kk3 = -k
          dd = ddd
          locf = 122
        endif
      endif
12    continue
C
      do 13 j = 1, m
      s = s1(m) + s3(j)
      k = nint(s*s2i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr2)
        if (ddd .lt. dd) then
          kk1 = m
          kk2 = -k
          kk3 = j
          dd = ddd
          locf = 131
        endif
      endif
      s = s1(m) - s3(j)
      k = nint(s*s2i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr2)
        if (ddd .lt. dd) then
          kk1 = m
          kk2 = -k
          kk3 = -j
          dd = ddd
          locf = 132
        endif
      endif
13    continue
C
      do 21 j = 1, m
      s = s1(j) + s2(m)
      k = nint(s*s3i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr3)
        if (ddd .lt. dd) then
          kk1 = j
          kk2 = m
          kk3 = -k
          dd = ddd
          locf = 211
        endif
      endif
      s = - s1(j) + s2(m)
      k = nint(s*s3i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr3)
        if (ddd .lt. dd) then
          kk1 = -j
          kk2 = m
          kk3 = -k
          dd = ddd
          locf = 212
        endif
      endif
21    continue
C
      do 23 j = 1, m
      s = s2(m) + s3(j)
      k = nint(s*s1i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr1)
        if (ddd .lt. dd) then
          kk1 = -k
          kk2 = m
          kk3 = j
          dd = ddd
          locf = 231
        endif
      endif
      s = s2(m) - s3(j)
      k = nint(s*s1i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr1)
        if (ddd .lt. dd) then
          kk1 = -k
          kk2 = m
          kk3 = -j
          dd = ddd
          locf = 232
        endif
      endif
23    continue
C
      do 31 j = 1, m
      s = s1(j) + s3(m)
      k = nint(s*s2i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr2)
        if (ddd .lt. dd) then
          kk1 = j
          kk2 = -k
          kk3 = m
          dd = ddd
          locf = 311
        endif
      endif
      s = -s1(j) + s3(m)
      k = nint(s*s2i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr2)
        if (ddd .lt. dd) then
          kk1 = -j
          kk2 = -k
          kk3 = m
          dd = ddd
          locf = 312
        endif
      endif
31    continue
C
      do 32 j = 1, m
      s = s2(j) + s3(m)
      k = nint(s*s1i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr1)
        if (ddd .lt. dd) then
          kk1 = -k
          kk2 = j
          kk3 = m
          dd = ddd
          locf = 321
        endif
      endif
      s = - s2(j) + s3(m)
      k = nint(s*s1i)
      if (abs(k) .le. m) then
        ddd = abs(s - k*sr1)
        if (ddd .lt. dd) then
          kk1 = -k
          kk2 = -j
          kk3 = m
          dd = ddd
          locf = 322
        endif
      endif
32    continue
C
      if (dd .lt. eps) then
C Recompte weighted sum with extended precision
        dd = abs(kk1*qsr1 + kk2*qsr2 + kk3*qsr3)
      endif
C Check for decrease of minimum
      if (dd .lt. d) then
        write (*,*) m, dd, kk1, kk2, kk3, locf
        d = dd
      endif
100   continue
      end