Klarner-Rado sequence

From Rosetta Code
Revision as of 10:36, 25 August 2022 by Hout (talk | contribs) (→‎{{header Haskell}}: Added a Haskell draft)
Klarner-Rado sequence is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Klarner-Rado sequences are a class of similar sequences that were studied by the mathematicians David Klarner and Richard Rado.

The most well known is defined as the thinnest strictly ascending sequence K which starts 1, then, for each element n, it will also contain somewhere in the sequence, 2 × n + 1 and 3 × n + 1.


So, the sequence K starts with 1.

Set n equal to the first element 1; the sequence will also contain 2 × n + 1 and 3 × n + 1, or 3 and 4.

Set n equal to the next element: 3, somewhere in the sequence it will contain 2 × n + 1 and 3 × n + 1, or 7 and 10.

Continue setting n equal to each element in turn to add to the sequence.


Task
  • Find and display the first one hundred elements of the sequence.
  • Find and display the one thousandth and ten thousandth elements of the sequence.

Preferably without needing to find an over abundance and sorting.


Stretch
  • Find and display the one hundred thousandth and one millionth elements of the sequence.


See also


ALGOL 68

Generates the sequence in order. Note that to run this with ALGOL 68G under Windows (and probably Linux), a large heap size must be specified on the command line, e.g.: -heap 1024M. <lang algol68>BEGIN # find elements of the Klarner-Rado sequence #

     #    - if n is an element, so are 2n + 1 and 3n + 1                    #
   INT max element = 100 000 000;
   [ 0 : max element ]BOOL kr; FOR i FROM LWB kr TO UPB kr DO kr[ i ] := FALSE OD;
   INT n21      := 3;                    # next 2n+1 value                  #
   INT n31      := 4;                    # next 3n+1 value                  #
   INT p2       := 1;                    # the n for the next 2n+1 value    #
   INT p3       := 1;                    # the n for the next 3n+1 value    #
   kr[ 1 ]      := TRUE;
   INT kr count := 0;
   INT max count = 1 000 000;
   FOR i WHILE kr count < max count DO
       IF i = n21 THEN                   # found the next 2n+1 value        #
           IF kr[ p2 ] THEN
               kr[ i ] := TRUE
           FI;
           p2  +:= 1;
           n21 +:= 2
       FI;
       IF i = n31 THEN                   # found the next 3n+1 value        #
           IF kr[ p3 ] THEN
               kr[ i ] := TRUE
           FI;
           p3  +:= 1;
           n31 +:= 3
       FI;
       IF kr[ i ] THEN
           kr count +:= 1;
           IF   kr count <= 100 THEN
               print( ( " ", whole( i, -3 ) ) );
               IF kr count MOD 20 = 0 THEN print( ( newline ) ) FI
           ELIF kr count =     1 000 THEN
               print( ( "The     thousandth element is ", whole( i, -8 ), newline ) )
           ELIF kr count =    10 000 THEN
               print( ( "The ten thousandth element is ", whole( i, -8 ), newline ) )
           ELIF kr count =   100 000 THEN
               print( ( "The 100-thousandth element is ", whole( i, -8 ), newline ) )
           ELIF kr count = 1 000 000 THEN
               print( ( "The      millionth element is ", whole( i, -8 ), newline ) )
           FI
       FI
   OD

END</lang>

Output:
   1   3   4   7   9  10  13  15  19  21  22  27  28  31  39  40  43  45  46  55
  57  58  63  64  67  79  81  82  85  87  91  93  94 111 115 117 118 121 127 129
 130 135 136 139 159 163 165 166 171 172 175 183 187 189 190 193 202 223 231 235
 237 238 243 244 247 255 256 259 261 262 271 273 274 279 280 283 319 327 331 333
 334 343 345 346 351 352 355 364 367 375 379 381 382 387 388 391 405 406 409 418
The     thousandth element is     8487
The ten thousandth element is   157653
The 100-thousandth element is  2911581
The      millionth element is 54381285

AppleScript

One way to test numbers for membership of the sequence is to feed them to a recursive handler which determines whether or not there's a Klarner-Rado route from them down to 0. It makes finding the elements in order simple, but takes about five and a half minutes to get to the millionth one.

<lang applescript>-- Is n in the Klarner-Rado sequence? -- Fully recursive: (* on isKlarnerRado(n)

   set n to n - 1
   return ((n = 0) or ((n mod 2 = 0) and (isKlarnerRado(n div 2))) or ((n mod 3 = 0) and (isKlarnerRado(n div 3))))

end isKlarnerRado *)

-- Optimised with tail call elimination. 90 seconds faster than the above in this script. on isKlarnerRado(n)

   set n to n - 1
   repeat
       if ((n = 0) or ((n mod 3 = 0) and (isKlarnerRado(n div 3)))) then
           return true
       else if (n mod 2 = 0) then
           set n to n div 2 - 1
       else
           return false
       end if
   end repeat

end isKlarnerRado

on join(lst, delim)

   set astid to AppleScript's text item delimiters
   set AppleScript's text item delimiters to delim
   set txt to lst as text
   set AppleScript's text item delimiters to astid
   return txt

end join

on task()

   set output to {"First 100 elements:"}
   set n to 0
   set |count| to 0
   set K to {}
   repeat until (|count| = 100)
       set n to n + 1
       if (isKlarnerRado(n)) then
           set K's end to n
           set |count| to |count| + 1
       end if
   end repeat
   repeat with i from 1 to 100 by 20
       set output's end to join(K's items i thru (i + 19), "  ")
   end repeat
   
   repeat with this in {{1000, "1,000th element: "}, {10000, "10,000th element: "}, {100000, "100,000th element: "}, ¬
       {1000000, "1,000,000th element: "}}
       set {target, spiel} to this
       repeat until (|count| = target)
           set n to n + 1
           if (isKlarnerRado(n)) then set |count| to |count| + 1
       end repeat
       set output's end to spiel & n
   end repeat
   
   return join(output, linefeed)

end task

task()</lang>

Output:

<lang applescript>"First 100 elements: 1 3 4 7 9 10 13 15 19 21 22 27 28 31 39 40 43 45 46 55 57 58 63 64 67 79 81 82 85 87 91 93 94 111 115 117 118 121 127 129 130 135 136 139 159 163 165 166 171 172 175 183 187 189 190 193 202 223 231 235 237 238 243 244 247 255 256 259 261 262 271 273 274 279 280 283 319 327 331 333 334 343 345 346 351 352 355 364 367 375 379 381 382 387 388 391 405 406 409 418 1,000th element: 8487 10,000th element: 157653 100,000th element: 2911581 1,000,000th element: 54381285"</lang>

Another way is to produce a list with Klarner-Rado elements in the slots indexed by those numbers and 'missing values' in the other slots. If a number being tested is exactly divisible by 2 or by 3, and the slot whose index is the result of the division contains a number instead of 'missing value', the tested number plus 1 is a Klarner-Rado element and should go into the slot it indexes. The list will contain vastly more 'missing values' than Klarner-Rado elements and it — or portions of it — ideally need to exist before the sequence elements are inserted. So while an overabundance and sorting of sequence elements isn't needed, an overabundance of 'missing values' is! The script below carries out the task in about 75 seconds on my current machine and produces the same output as above.

<lang applescript>on KlarnerRadoSequence(target)

   -- To find a million KR numbers with this method nominally needs a list with at least 54,381,285
   -- slots. But the number isn't known in advance, "growing" a list to the required length would
   -- take forever, and accessing its individual items could also take a while. Instead, K will
   -- here contain reasonably sized and quickly produced sublists which are added as needed.
   -- The 1-based referencing calculations will be ((index - 1) div sublistLen + 1) for the sublist
   -- and ((index - 1) mod sublistLen + 1) for the slot within it.
   set sublistLen to 10000
   script o
       property spare : makeList(sublistLen, missing value)
       property K : {spare's items}
   end script
   
   -- Insert the initial 1, start the KR counter at 1, start the number-to-test variable at 2.
   set {o's K's beginning's 1st item, |count|, n} to {1, 1, 2}
   -- Test the first, second, third, and fifth of every six consecutive numbers starting at 2.
   -- These are known to be divisible by 2 or by 3 and the fourth and sixth known not to be.
   -- If the item at index (n div 2) or index (n div 3) isn't 'missing value', it's a number,
   -- so insert (n + 1) at index (n + 1).
   if (|count| < target) then ¬
       repeat -- Per increase of n by 6.
           -- The first of the six numbers in this range is divisible by 2.
           -- Precalculate (index - 1) for index (n div 2) to reduce code clutter and halve calculation time.
           set pre to n div 2 - 1
           if (o's K's item (pre div sublistLen + 1)'s item (pre mod sublistLen + 1) is missing value) then
           else
               -- Insert (n + 1) at index (n + 1). The 'n's in the reference calculations are for ((n + 1) - 1)!
               set o's K's item (n div sublistLen + 1)'s item (n mod sublistLen + 1) to n + 1
               set |count| to |count| + 1
               if (|count| = target) then exit repeat
           end if
           -- The second number of the six is divisible by 3.
           set n to n + 1
           set pre to n div 3 - 1
           if (o's K's item (pre div sublistLen + 1)'s item (pre mod sublistLen + 1) is missing value) then
           else
               set o's K's item (n div sublistLen + 1)'s item (n mod sublistLen + 1) to n + 1
               set |count| to |count| + 1
               if (|count| = target) then exit repeat
           end if
           -- The third is divisible by 2.
           set n to n + 1
           set pre to n div 2 - 1
           if (o's K's item (pre div sublistLen + 1)'s item (pre mod sublistLen + 1) is missing value) then
           else
               set o's K's item (n div sublistLen + 1)'s item (n mod sublistLen + 1) to n + 1
               set |count| to |count| + 1
               if (|count| = target) then exit repeat
           end if
           -- The fifth is divisible by both 2 and 3.
           set n to n + 2
           set pre2 to n div 2 - 1
           set pre3 to n div 3 - 1
           if ((o's K's item (pre2 div sublistLen + 1)'s item (pre2 mod sublistLen + 1) is missing value) and ¬
               (o's K's item (pre3 div sublistLen + 1)'s item (pre3 mod sublistLen + 1) is missing value)) then
           else
               set o's K's item (n div sublistLen + 1)'s item (n mod sublistLen + 1) to n + 1
               set |count| to |count| + 1
               if (|count| = target) then exit repeat
           end if
           
           -- Advance to the first number of the next six.
           set n to n + 2
           -- If another sublist is about to be needed, append one to the end of K.
           if ((n + 6) mod sublistLen < 6) then copy o's spare to o's K's end
       end repeat
   
   -- Once the target's reached, replace the sublists with lists containing just the numbers.
   set sublistCount to (count o's K)
   repeat with i from 1 to sublistCount
       set o's K's item i to o's K's item i's numbers
   end repeat
   -- Concatenate the number lists to leave K as a single list of numbers.
   repeat while (sublistCount > 1)
       set o's spare to o's K
       set o's K to {}
       repeat with i from 2 to sublistCount by 2
           set end of o's K to o's spare's item (i - 1) & o's spare's item i
       end repeat
       if (i < sublistCount) then set o's K's last item to o's K's end & o's spare's end
       set sublistCount to sublistCount div 2
   end repeat
   set o's K to o's K's beginning
   
   return o's K

end KlarnerRadoSequence

on makeList(len, content)

   script o
       property lst : {}
   end script
   if (len > 0) then
       set o's lst's end to content
       set |count| to 1
       repeat until (|count| + |count| > len)
           set o's lst to o's lst & o's lst
           set |count| to |count| + |count|
       end repeat
       if (|count| < len) then set o's lst to o's lst & o's lst's items 1 thru (len - |count|)
   end if
   
   return o's lst

end makeList

on join(lst, delim)

   set astid to AppleScript's text item delimiters
   set AppleScript's text item delimiters to delim
   set txt to lst as text
   set AppleScript's text item delimiters to astid
   return txt

end join

on task()

   script o
       property K : KlarnerRadoSequence(1000000)
   end script
   set output to {"First 100 elements:"}
   repeat with i from 1 to 100 by 20
       set output's end to join(o's K's items i thru (i + 19), "  ")
   end repeat
   set output's end to "1,000th element: " & o's K's item 1000
   set output's end to "10,000th element: " & o's K's item 10000
   set output's end to "100,000th element: " & o's K's item 100000
   set output's end to "1,000,000th element: " & o's K's item 1000000
   
   return join(output, linefeed)

end task

task()</lang>

F#

<lang fsharp> // Klarner-Rado sequence. Nigel Galloway: August 19th., 20 let kr()=Seq.unfold(fun(n,g)->Some(g|>Seq.filter((>)n)|>Seq.sort,(n*2L+1L,seq[for n in g do yield n; yield n*2L+1L; yield n*3L+1L]|>Seq.filter((<=)n)|>Seq.distinct)))(3L,seq[1L])|>Seq.concat let n=kr()|>Seq.take 1000000|>Array.ofSeq in n|>Array.take 100|>Array.iter(printf "%d "); printfn "\nkr[999] is %d\nkr[9999] is %d\nkr[99999] is %d\nkr[999999] is %d" n.[999] n.[9999] n.[99999] n.[999999] </lang>

Output:
1 3 4 7 9 10 13 15 19 21 22 27 28 31 39 40 43 45 46 55 57 58 63 64 67 79 81 82 85 87 91 93 94 111 115 117 118 121 127 129 130 135 136 139 159 163 165 166 171 172 175 183 187 189 190 193 202 223 231 235 237 238 243 244 247 255 256 259 261 262 271 273 274 279 280 283 319 327 331 333 334 343 345 346 351 352 355 364 367 375 379 381 382 387 388 391 405 406 409 418
kr[999] is 8487
kr[9999] is 157653
kr[99999] is 2911581
kr[999999] is 54381285

FreeBASIC

Translation of: Phix

<lang freebasic>#include "string.bi"

Dim As Integer limit = 1e8 Dim As Boolean kr(limit) kr(1) = True Dim As Integer n21 = 3, n31 = 4, p2 = 1, p3 = 1, count = 0, np10 = 1000 Print "First 100 Klarner-Rado sequence numbers:" For i As Integer = 1 To limit

   If i = n21 Then
       If kr(p2) Then kr(i) = True
       p2  += 1
       n21 += 2
   End If
   If i = n31 Then
       If kr(p3) Then kr(i) = True
       p3  += 1
       n31 += 3
   End If
   If kr(i) Then
       count += 1
       If count <= 100 Then
           Print Using "####"; i;
           If(count Mod 20) = 0 Then Print
       Elseif count = np10 Then
           Print "The "; Format(count, "#,###,###"); "th Klarner-Rado number is "; Format(i, "#,###,###")
           np10 *= 10
       End If
   End If

Next i Sleep</lang>

Output:
Same as Phix entry.

Haskell

import Data.List (intercalate)
import Data.List.Ordered (union)
import Data.List.Split (chunksOf)
import Text.Printf

----------------------- KLARNER-RADO ---------------------

klarnerRado :: (Num a, Ord a, Enum a) => [a]
klarnerRado =
  1 :
  union
    (succ . (* 2) <$> klarnerRado)
    (succ . (* 3) <$> klarnerRado)


--------------------------- TEST -------------------------
main :: IO ()
main = do
  putStrLn "First one hundred elements of the sequence:\n"
  mapM_
    putStrLn
    ( intercalate "  "
        <$> chunksOf
          10
          (printf "%3d" <$> take 100 klarnerRado)
    )

  putStrLn "\nKth and 10Kth elements of the sequence:\n"
  mapM_
    ( putStrLn
        . (<*>)
          (flip (printf "%5dth %s %7d") " ->")
          (klarnerRado !!)
    )
    [1000, 10000]
Output:
First one hundred elements of the sequence:

  1    3    4    7    9   10   13   15   19   21
 22   27   28   31   39   40   43   45   46   55
 57   58   63   64   67   79   81   82   85   87
 91   93   94  111  115  117  118  121  127  129
130  135  136  139  159  163  165  166  171  172
175  183  187  189  190  193  202  223  231  235
237  238  243  244  247  255  256  259  261  262
271  273  274  279  280  283  319  327  331  333
334  343  345  346  351  352  355  364  367  375
379  381  382  387  388  391  405  406  409  418

Kth and 10Kth elements of the sequence:

 1000th  ->    8488
10000th  ->  157654

J

Implementation:<lang J>krsieve=: {{

 for_i. i.<.-:#b=. (1+y){.0 1 do.
   if. i{b do. b=. 1 ((#~ y&>)1+2 3*i)} b end.
 end.
 I.b

}}</lang>

Task examples (including stretch):<lang J> #kr7e7=: krsieve 7e7 1215307

  100{.kr7e7

1 3 4 7 9 10 13 15 19 21 22 27 28 31 39 40 43 45 46 55 57 58 63 64 67 79 81 82 85 87 91 93 94 111 115 117 118 121 127 129 130 135 136 139 159 163 165 166 171 172 175 183 187 189 190 193 202 223 231 235 237 238 243 244 247 255 256 259 261 262 271 273 274 279 280 283 319 327 331 333 334 343 345 346 351 352 355 364 367 375 379 381 382 387 388 391 405 406 409 418

  (1e3-1){kr7e7

8487

  (1e4-1){kr7e7

157653

  (1e5-1){kr7e7

2911581

  (1e6-1){kr7e7

54381285</lang>

Julia

<lang julia>using Formatting

function KlarnerRado(N)

   K = [1]
   for i in 1:N
       j = K[i]
       firstadd, secondadd = 2j + 1, 3j + 1
       if firstadd < K[end]
           pos = findlast(<(firstadd), K) + 1
           K[pos] != firstadd && insert!(K, pos, firstadd)
       elseif K[end] != firstadd
           push!(K, firstadd)
       end
       if secondadd < K[end]
           pos = findlast(<(secondadd), K) + 1
           K[pos] != secondadd && insert!(K, pos, secondadd)
       elseif K[end] != secondadd
           push!(K, secondadd)
       end
   end
   return K

end

kr1m = KlarnerRado(1_000_000)

println("First 100 Klarner-Rado numbers:") foreach(p -> print(rpad(p[2], 4), p[1] % 20 == 0 ? "\n" : ""), enumerate(kr1m[1:100])) foreach(n -> println("The ", format(n, commas=true), "th Klarner-Rado number is ",

  format(kr1m[n], commas=true)), [1000, 10000, 100000, 1000000])

</lang>

Output:
First 100 Klarner-Rado numbers:
1   3   4   7   9   10  13  15  19  21  22  27  28  31  39  40  43  45  46  55  
57  58  63  64  67  79  81  82  85  87  91  93  94  111 115 117 118 121 127 129
130 135 136 139 159 163 165 166 171 172 175 183 187 189 190 193 202 223 231 235
237 238 243 244 247 255 256 259 261 262 271 273 274 279 280 283 319 327 331 333 
334 343 345 346 351 352 355 364 367 375 379 381 382 387 388 391 405 406 409 418
The 1,000th Klarner-Rado number is 8,487
The 10,000th Klarner-Rado number is 157,653
The 100,000th Klarner-Rado number is 2,911,581
The 1,000,000th Klarner-Rado number is 54,381,285

Faster version

Probably does get an overabundance, but no sorting. `falses()` is implemented as a bit vector, so huge arrays are not needed. From ALGOL. <lang ruby>using Formatting

function KlamerRado(N)

   kr = falses(100 * N)
   kr[1] = true
   for i in 1:30N
       if kr[i]
           kr[2i + 1] = true
           kr[3i + 1] = true
       end
   end
   return [i for i in eachindex(kr) if kr[i]]

end

kr1m = KlamerRado(1000000)

println("First 100 Klarner-Rado numbers:") foreach(p -> print(rpad(p[2], 4), p[1] % 20 == 0 ? "\n" : ""), enumerate(kr1m[1:100])) foreach(n -> println("The ", format(n, commas=true), "th Klarner-Rado number is ",

  format(kr1m[n], commas=true)), [1000, 10000, 100000, 1000000])

</lang>

Output:

same as above version.

Phix

Translation of: ALGOL_68
with javascript_semantics
constant limit = iff(platform()=JS?10_000_000:100_000_000)
sequence kr = repeat(false,limit); kr[1] = true
integer n21 = 3, n31 = 4, p2 = 1, p3 = 1, count = 0, np10 = 1_000
for i=1 to limit do
    if i = n21 then
        if kr[p2] then
            kr[i] := true
        end if
        p2  += 1
        n21 += 2
    end if
    if i = n31 then
        if kr[p3] then
            kr[i] := true
        end if
        p3  += 1
        n31 += 3
    end if
    if kr[i] then
        count += 1
        if count <= 100 then
            printf(1,"%4d%s",{i,iff(mod(count,20)=0?"\n":"")})
        elsif count = np10 then
            printf(1,"The %,d%s Klarner-Rado number is %,d\n",{count,ord(count),i})
            np10 *= 10
        end if
    end if
end for
Output:
   1   3   4   7   9  10  13  15  19  21  22  27  28  31  39  40  43  45  46  55
  57  58  63  64  67  79  81  82  85  87  91  93  94 111 115 117 118 121 127 129
 130 135 136 139 159 163 165 166 171 172 175 183 187 189 190 193 202 223 231 235
 237 238 243 244 247 255 256 259 261 262 271 273 274 279 280 283 319 327 331 333
 334 343 345 346 351 352 355 364 367 375 379 381 382 387 388 391 405 406 409 418
The 1,000th Klarner-Rado number is 8,487
The 10,000th Klarner-Rado number is 157,653
The 100,000th Klarner-Rado number is 2,911,581
The 1,000,000th Klarner-Rado number is 54,381,285

Unfortunately JavaScript can't quite cope/runs out of memory with a 100 million element sieve, so it only goes to the 10^5th under p2js.

PL/M

Works with: 8080 PL/M Compiler

... under CP/M (or an emulator)


As PL/M only handles unsigned 8 and 16 bit integers, this only finds the first 1000 elements. This is based on the Algol 68 sample, but as with the VB.NET sample, uses a "bit vector" (here an array of bytes) - as suggested by the Julia sample. <lang pli>100H: /* EIND ELEMENTS OF THE KLARNER-RADO SEQUENCE */

     /*    - IF N IS AN ELEMENT, SO ARE 2N+1 AND 3N+!, 1 IS AN ELEMENT     */
  /* CP/M SYSTEM CALL AND I/O ROUTINES                                     */
  BDOS:      PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
  PR$CHAR:   PROCEDURE( C ); DECLARE C BYTE;    CALL BDOS( 2, C );  END;
  PR$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S );  END;
  PR$NL:     PROCEDURE;   CALL PR$CHAR( 0DH ); CALL PR$CHAR( 0AH ); END;
  PR$NUMBER: PROCEDURE( N ); /* PRINTS A NUMBER IN THE MINIMUN FIELD WIDTH */
     DECLARE N ADDRESS;
     DECLARE V ADDRESS, N$STR ( 6 )BYTE, W BYTE;
     V = N;
     W = LAST( N$STR );
     N$STR( W ) = '$';
     N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
     DO WHILE( ( V := V / 10 ) > 0 );
        N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
     END;
     CALL PR$STRING( .N$STR( W ) );
  END PR$NUMBER;
  /* TASK                                                                  */
  DECLARE MAX$COUNT LITERALLY '1$000';
  DECLARE BIT ( 8 )BYTE INITIAL( 128, 1, 2, 4, 8, 16, 32, 64 );
  DECLARE ( N21, N31, P2, P3, KR$COUNT, I, I$3 ) ADDRESS;
  DECLARE KR ( 1251 )BYTE; DO I = 0 TO LAST( KR ); KR( I ) = 0; END;
  KR( 0 ) = BIT( 1 );
  KR$COUNT = 0; N21 = 3; N31 = 4; P2, P3 = 1;
  I = 0;
  DO WHILE( KR$COUNT < MAX$COUNT );
     I = I + 1;
     I$3 = SHR( I, 3 );
     IF I = N21 THEN DO;                         /* I IS 2N+1 WHERE N = P2 */
        IF ( KR( SHR( P2, 3 ) ) AND BIT( P2 AND 7 ) ) <> 0 THEN DO;
           KR( I$3 ) = KR( I$3 ) OR BIT( I AND 7 );
        END;
        P2  = P2  + 1;
        N21 = N21 + 2;
     END;
     IF I = N31 THEN DO;                         /* I IS 3N+1 WHERE N = P3 */
        IF ( KR( SHR( P3, 3 ) ) AND BIT( P3 AND 7 ) ) <> 0 THEN DO;
           KR( I$3 ) = KR( I$3 ) OR BIT( I AND 7 );
        END;
        P3  = P3  + 1;
        N31 = N31 + 3;
     END;
     IF ( KR( I$3 ) AND BIT( I AND 7 ) ) <> 0 THEN DO;
        KR$COUNT = KR$COUNT + 1;
        IF KR$COUNT <= 100 THEN DO;
           CALL PR$CHAR( ' ' );
           IF I <  10 THEN CALL PR$CHAR( ' ' );
           IF I < 100 THEN CALL PR$CHAR( ' ' );
           CALL PR$NUMBER( I );
           IF KR$COUNT MOD 20 = 0 THEN CALL PR$NL;
           END;
        ELSE IF KR$COUNT = MAX$COUNT THEN DO;
           CALL PR$STRING( .'ELEMENT $' );
           CALL PR$NUMBER( MAX$COUNT );
           CALL PR$STRING( .' IS: $' );
           CALL PR$NUMBER( I );
        END;
     END;
  END;

EOF</lang>

Output:
   1   3   4   7   9  10  13  15  19  21  22  27  28  31  39  40  43  45  46  55
  57  58  63  64  67  79  81  82  85  87  91  93  94 111 115 117 118 121 127 129
 130 135 136 139 159 163 165 166 171 172 175 183 187 189 190 193 202 223 231 235
 237 238 243 244 247 255 256 259 261 262 271 273 274 279 280 283 319 327 331 333
 334 343 345 346 351 352 355 364 367 375 379 381 382 387 388 391 405 406 409 418
ELEMENT 1000 IS: 8487

Python

<lang python>def KlarnerRado(N):

   K = [1]
   for i in range(N):
       j = K[i]
       firstadd, secondadd = 2 * j + 1, 3 * j + 1
       if firstadd < K[-1]:
           for pos in range(len(K)-1, 1, -1):
               if K[pos] < firstadd < K[pos + 1]:
                   K.insert(pos + 1, firstadd)
                   break
       elif firstadd > K[-1]:
           K.append(firstadd)
       if secondadd < K[-1]:
           for pos in range(len(K)-1, 1, -1):
               if K[pos] < secondadd < K[pos + 1]:
                   K.insert(pos + 1, secondadd)
                   break
       elif secondadd > K[-1]:
           K.append(secondadd)
   return K

kr1m = KlarnerRado(100_000)

print('First 100 Klarner-Rado sequence numbers:') for idx, v in enumerate(kr1m[:100]):

   print(f'{v: 4}', end='\n' if (idx + 1) % 20 == 0 else )

for n in [1000, 10_000, 100_000]:

   print(f'The {n :,}th Klarner-Rado number is {kr1m[n-1] :,}')

</lang>

Output:
First 100 Klarner-Rado sequence numbers:
   1   3   4   7   9  10  13  15  19  21  22  27  28  31  39  40  43  45  46  55
  57  58  63  64  67  79  81  82  85  87  91  93  94 111 115 117 118 121 127 129
 130 135 136 139 159 163 165 166 171 172 175 183 187 189 190 193 202 223 231 235
 237 238 243 244 247 255 256 259 261 262 271 273 274 279 280 283 319 327 331 333
 334 343 345 346 351 352 355 364 367 375 379 381 382 387 388 391 405 406 409 418
The 1,000th Klarner-Rado number is 8,487
The 10,000th Klarner-Rado number is 157,653
The 100,000th Klarner-Rado number is 2,911,581

faster version

<lang python>from numpy import ndarray

def KlarnerRado(N):

   kr = ndarray(100 * N, dtype=bool)
   kr[1] = True
   for i in range(30 * N):
       if kr[i]:
           kr[2 * i + 1] = True
           kr[3 * i + 1] = True
   return [i for i in range(100 * N) if kr[i]]

kr1m = KlarnerRado(1_000_000)

print('First 100 Klarner-Rado sequence numbers:') for idx, v in enumerate(kr1m[:100]):

   print(f'{v: 4}', end='\n' if (idx + 1) % 20 == 0 else )

for n in [1000, 10_000, 100_000, 1_000_000]:

   print(f'The {n :,}th Klarner-Rado number is {kr1m[n-1] :,}')

</lang>

Output:

Same as previous version.

Raku

<lang perl6>sub Klarner-Rado ($n) {

   my @klarner-rado = 1;
   my @next = x2, x3;
   loop {
       @klarner-rado.push: my $min = @next.min;
       @next[0] = x2 if @next[0] == $min;
       @next[1] = x3 if @next[1] == $min;
       last if +@klarner-rado > $n.max;
   }
   sub x2 { state $i = 0; @klarner-rado[$i++] × 2 + 1 }
   sub x3 { state $i = 0; @klarner-rado[$i++] × 3 + 1 }
   @klarner-rado[|$n]

}

use Lingua::EN::Numbers;

put "First 100 elements of Klarner-Rado sequence:"; put Klarner-Rado(^100).batch(10)».fmt("%3g").join: "\n"; put ; put (($_».Int».&ordinal».tc »~» ' element: ').List Z~ |(List.new: (Klarner-Rado ($_ »-» 1))».&comma)

   given $(1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6)).join: "\n";</lang>
Output:
First 100 elements of Klarner-Rado sequence:
  1   3   4   7   9  10  13  15  19  21
 22  27  28  31  39  40  43  45  46  55
 57  58  63  64  67  79  81  82  85  87
 91  93  94 111 115 117 118 121 127 129
130 135 136 139 159 163 165 166 171 172
175 183 187 189 190 193 202 223 231 235
237 238 243 244 247 255 256 259 261 262
271 273 274 279 280 283 319 327 331 333
334 343 345 346 351 352 355 364 367 375
379 381 382 387 388 391 405 406 409 418

First element: 1
Tenth element: 21
One hundredth element: 418
One thousandth element: 8,487
Ten thousandth element: 157,653
One hundred thousandth element: 2,911,581
One millionth element: 54,381,285

Visual Basic .NET

Based on the ALGOL 68 sample, but using a "bit vector" (array of integers), as suggested by the Julia sample. Simplifies printing "the powers of ten" elements, as in the Phix sample. <lang vbnet>Option Strict On Option Explicit On

Imports System.IO

Module KlarnerRado

   Private Const bitsWidth As Integer = 31
   Private bitMask() As Integer = _
       New Integer(){ 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024 _
                    , 2048, 4096, 8192, 16384, 32768, 65536, 131072 _
                    , 262144, 524288, 1048576, 2097152, 4194304, 8388608 _
                    , 16777216, 33554432, 67108864, 134217728, 268435456 _
                    , 536870912, 1073741824 _
                    }
   Private Const maxElement As Integer = 1100000000


   Private Function BitSet(bit As Integer, v As Integer) As Boolean
       Return (v And bitMask(bit - 1)) <> 0
   End Function
   Private Function SetBit(ByVal bit As Integer, ByVal b As Integer) As Integer
       Return b Or bitMask(bit - 1)
   End Function
   Public Sub Main
       Dim  kr(maxElement \ bitsWidth) As Integer
       For i As Integer = 0 To kr.Count() - 1
           kr(i) = 0
       Next  i
       Dim krCount As Integer =  0    ' count of sequende elements
       Dim n21 As Integer = 3         ' next 2n+1 value
       Dim n31 As Integer = 4         ' next 3n+1 value
       Dim p10 As Integer = 1000      ' next power of ten to show the count/element at
       Dim iBit As Integer = 0        ' i Mod bitsWidth  - used to set the ith bit
       Dim iOverBw As Integer = 0     ' i \ bitsWidth    - used to set the ith bit
       '   p2                         ' the n for the 2n+1 value
       Dim p2Bit As Integer = 1       ' p2 Mod bitsWidth  - used to set the p2th bit
       Dim p2OverBw As Integer = 0    ' p2 \ bitsWidth    - used to set the p2th bit
       '   p3                         ' the n for the 3n+1 value
       Dim p3Bit As Integer = 1       ' p3 Mod bitsWidth  - used to set the p3th bit
       Dim p3OverBw As Integer = 0    ' p3 \ bitsWidth    - used to set the p3th bit
       kr(0) = SetBit(1, kr(0))
       Dim kri As Boolean = True
       Dim lastI As Integer = 0
       For i As Integer = 1 To  maxElement
           iBit += 1
           If iBit > bitsWidth Then
               iBit = 1
               iOverBw += 1
           End If
           If i = n21 Then            ' found the next 2n+1 value
               If BitSet(p2Bit, kr(p2OverBw)) Then
                   kri = True
               End If
               p2Bit += 1
               If p2Bit > bitsWidth Then
                   p2Bit = 1
                   p2OverBw += 1
               End If
               n21 += 2
           End If
           If i = n31 Then            ' found the next 3n+1 value
               If BitSet(p3Bit, kr(p3OverBw)) Then
                   kri = True
               End If
               p3Bit += 1
               If p3Bit > bitsWidth Then
                   p3Bit = 1
                   p3OverBw += 1
               End If
               n31 += 3
           End If
           If kri Then
               lastI = i
               kr(iOverBw) = SetBit(iBit, kr(iOverBw))
               krCount += 1
               If krCount <= 100 Then
                   Console.Out.Write(" " & i.ToString().PadLeft(3))
                   If krCount Mod 20 = 0 Then
                       Console.Out.WriteLine()
                   End If
               ElseIf krCount = p10 Then
                   Console.Out.WriteLine("Element " & p10.ToString().PadLeft(10) & " is " & i.ToString().PadLeft(10))
                   p10 *= 10
               End If
               kri = False
           End If
       Next  i
       Console.Out.WriteLine("Element " & krCount.ToString().PadLeft(10) & " is " & lastI.ToString().PadLeft(10))
   End Sub

End Module</lang>

Output:
   1   3   4   7   9  10  13  15  19  21  22  27  28  31  39  40  43  45  46  55
  57  58  63  64  67  79  81  82  85  87  91  93  94 111 115 117 118 121 127 129
 130 135 136 139 159 163 165 166 171 172 175 183 187 189 190 193 202 223 231 235
 237 238 243 244 247 255 256 259 261 262 271 273 274 279 280 283 319 327 331 333
 334 343 345 346 351 352 355 364 367 375 379 381 382 387 388 391 405 406 409 418
Element       1000 is       8487
Element      10000 is     157653
Element     100000 is    2911581
Element    1000000 is   54381285
Element   10000000 is 1031926801
Element   10543878 is 1099640002

Wren

Version 1

Library: Wren-sort
Library: Wren-fmt

There's no actual sorting here. The Find class (and its binary search methods) just happen to be in Wren-sort. <lang ecmascript>import "./sort" for Find import "./fmt" for Fmt

var klarnerRado = Fn.new { |limit|

   var kr = [1]
   for (i in 0...limit) {
       var n = kr[i]
       for (e in [2*n + 1, 3*n + 1]) {
           if (e > kr[-1]) {
               kr.add(e)
           } else {
               var ix = Find.nearest(kr, e)  // binary search
               if (kr[ix] != e) kr.insert(ix, e)
           }
       }
   }
   return kr[0...limit]

}

var kr = klarnerRado.call(1e6) System.print("First 100 elements of Klarner-Rado sequence:") Fmt.tprint("$3d ", kr[0..99], 10) System.print() var limits = [1, 10, 1e2, 1e3, 1e4, 1e5, 1e6] for (limit in limits) {

   Fmt.print("The $,r element: $,d", limit, kr[limit-1])

} </lang>

Output:
First 100 elements of Klarner-Rado sequence:
  1    3    4    7    9   10   13   15   19   21  
 22   27   28   31   39   40   43   45   46   55  
 57   58   63   64   67   79   81   82   85   87  
 91   93   94  111  115  117  118  121  127  129  
130  135  136  139  159  163  165  166  171  172  
175  183  187  189  190  193  202  223  231  235  
237  238  243  244  247  255  256  259  261  262  
271  273  274  279  280  283  319  327  331  333  
334  343  345  346  351  352  355  364  367  375  
379  381  382  387  388  391  405  406  409  418  

The 1st element: 1
The 10th element: 21
The 100th element: 418
The 1,000th element: 8,487
The 10,000th element: 157,653
The 100,000th element: 2,911,581
The 1,000,000th element: 54,381,285

Version 2 (faster)

Translation of: ALGOL 68
Library: Wren-array

Takes around 13.8 seconds to find the millionth element compared to 48 seconds for the first version.

As in the case of some other examples, uses a bit array to save memory. However, whilst the BitArray class crams 32 bools into a 64 bit float (Wren doesn't have integers as such), it is significantly slower to index than the native List class. If the latter is used instead, the time comes down to 5.3 seconds.

Although not shown here, if the size of the BitArray is increased to 1.1 billion and 'max' to 1e7, then the 10 millionth element (1,031,926,801) will eventually be found but takes 4 minutes 50 seconds to do so. <lang ecmascript>import "./array" for BitArray import "./fmt" for Fmt

var kr = BitArray.new(1e8) var first100 = List.filled(100, 0) var pow10 = {} kr[1] = true var n21 = 3 var n31 = 4 var p2 = 1 var p3 = 1 var count = 0 var max = 1e6 var i = 1 var limit = 1 while (count < max) {

   if (i == n21) {
       if (kr[p2]) kr[i] = true
       p2 = p2 + 1
       n21 = n21 + 2
   }
   if (i == n31) {
       if (kr[p3]) kr[i] = true
       p3 = p3 + 1
       n31 = n31 + 3
   }
   if (kr[i]) {
       count = count + 1
       if (count <= 100) first100[count-1] = i
       if (count == limit) {
           pow10[count] = i
           if (limit == max) break
           limit = 10 * limit
       }
   }
   i = i + 1

}

System.print("First 100 elements of Klarner-Rado sequence:") Fmt.tprint("$3d ", first100, 10) System.print() limit = 1 while (limit <= max) {

   Fmt.print("The $,r element: $,d", limit, pow10[limit])
   limit = 10 * limit

}</lang>

Output:
Same as Version 1.