Jump to content

Greatest prime dividing the n-th cubefree number

From Rosetta Code
Task
Greatest prime dividing the n-th cubefree number
You are encouraged to solve this task according to the task description, using any language you may know.
Definitions

A cubefree number is a positive integer whose prime factorization does not contain any third (or higher) power factors. If follows that all primes are trivially cubefree and the first cubefree number is 1 because it has no prime factors.

Let a[n] be the greatest prime dividing the n-th cubefree number for n >= 2. By convention, let a[1] = 1 even though the first cubefree number, 1, has no prime factors.


Examples

a[2] is clearly 2 because it is the second cubefree number and also prime. The fourth cubefree number is 4 and it's highest prime factor is 2, hence a[4] = 2.


Task

Compute and show on this page the first 100 terms of a[n]. Also compute and show the 1,000th, 10,000th and 100,000th members of the sequence.


Stretch

Compute and show the 1 millionth and 10 millionth terms of the sequence.

This may take a while for interpreted languages.


Reference

ALGOL 68

Translation of: EasyLang – this is a line-by-line translation with a minor optimisation in the prime factorisation routine

Relative to the EasyLang original, the Algol 68 version has declarations, more punctuation and keywords instead of ".". Also, a routine cam't modify its parameters, so a temporary is needed for the "num" parameter of maxprimecubefree (the contents of reference parameters can be changed, but that would modify the the thing that the routine was called with).
The outer BEGIN and END can be omitted if your Algol 68 compiler/interpreter allows a particular program to be a serial clause (e.g., ALGOL 68 Genie)

BEGIN # find the largest prime factor of some cubefree numbers - translation of EasyLang #
   PROC maxprimcubefree = ( INT numin )INT: BEGIN
      INT num := numin, t := 2, count := 0;
      WHILE t * t <= num DO
         IF num MOD t = 0 THEN
            count +:= 1;
            IF count = 3 THEN
               num := 0
            FI;
            num OVERAB t
         ELSE
            count := 0;
            t +:= IF ODD t THEN 2 ELSE 1 FI
         FI
      OD;
      IF count = 2 AND t = num THEN
         num := 0
      FI;
      num
   END;
   INT i := 1, a count := 0;
   WHILE a count < 100 DO
      INT h = maxprimcubefree( i );
      IF h /= 0 THEN
         a count +:= 1;
         print( ( whole( h, 0 ), " " ) )
      FI;
      i +:= 1
   OD;
   print( ( newline ) );
   WHILE a count < 100 000 DO
      INT h = maxprimcubefree( i );
      IF h /= 0 THEN
         a count +:= 1;
         IF a count = 1 000 OR a count = 10 000 OR a count = 100 000 THEN
            print( ( whole( h, 0 ), newline ) )
         FI
      FI;
      i +:= 1
   OD
END
Output:
1 2 3 2 5 3 7 3 5 11 3 13 7 5 17 3 19 5 7 11 23 5 13 7 29 5 31 11 17 7 3 37 19 13 41 7 43 11 5 23 47 7 5 17 13 53 11 19 29 59 5 61 31 7 13 11 67 17 23 7 71 73 37 5 19 11 13 79 41 83 7 17 43 29 89 5 13 23 31 47 19 97 7 11 5 101 17 103 7 53 107 109 11 37 113 19 23 29 13 59
109
101
1693

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>

std::vector<uint32_t> prime_factors(uint32_t n) {
    std::vector<uint32_t> factors = { };
    while ( n % 2 == 0 ) {
    	factors.emplace_back(2);
    	n /= 2;
    }
    for ( uint32_t i = 3; i <= std::sqrt(n); i += 2 ) {
    	while ( n % i == 0 ) {
    		factors.emplace_back(i);
    		n /= i;
    	}
    }
    if ( n > 2 ) {
    	factors.emplace_back(n);
    }
    return factors;
}

int main() {
	const uint32_t maximum = 10000000;
	uint32_t count = 1;
	uint32_t i = 2;
	const uint32_t lower_limit = 100;
	uint32_t upper_limit = 1000;
	std::vector<uint32_t> first_hundred = { 1 };

	while ( count < maximum ) {
	    bool cube_free = false;
	    std::vector<uint32_t> factors = prime_factors(i);

	    if ( factors.size() < 3 ) {
	        cube_free = true;
	    } else {
	    	cube_free = true;
	    	for ( uint32_t i = 2; i < factors.size(); ++i ) {
	    		if ( factors[i - 2] == factors[i - 1] && factors[i - 1] == factors[i] ) {
	    			cube_free = false;
	    			break;
	    		}
	    	}
	    }

	    if ( cube_free ) {
	    	if ( count < lower_limit ) {
	    		first_hundred.emplace_back(factors.back());
	    	}
	    	count += 1;
	    	if ( count == lower_limit ) {
	    		std::cout << "The first " << lower_limit << " terms of a370833 are:" << "\n";
	    		for ( uint32_t i = 0; i < 100; ++i ) {
	    			std::cout << std::setw(3) << first_hundred[i] << ( i % 10 == 9 ? "\n" : " " );
	    		}
	    		std::cout << "\n";
	    	} else if ( count == upper_limit ) {
	    		std::cout << "The " << count << "th term of a370833 is " << factors.back() << "\n";
	    		upper_limit *= 10;
	    	}
	    }

	    i += 1;
	 }
}
Output:
The first 100 terms of a370833 are:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The 1000th term of a370833 is 109
The 10000th term of a370833 is 101
The 100000th term of a370833 is 1693
The 1000000th term of a370833 is 1202057
The 10000000th term of a370833 is 1202057

Dart

Translation of: Wren
import 'dart:math';

void main() {
  List<int> res = [1];
  int count = 1;
  int i = 2;
  int lim1 = 100;
  int lim2 = 1000;
  double max = 1e7;
  var t0 = DateTime.now();

  while (count < max) {
    bool cubeFree = false;
    List<int> factors = primeFactors(i);
    if (factors.length < 3) {
      cubeFree = true;
    } else {
      cubeFree = true;
      for (int i = 2; i < factors.length; i++) {
        if (factors[i - 2] == factors[i - 1] && factors[i - 1] == factors[i]) {
          cubeFree = false;
          break;
        }
      }
    }
    if (cubeFree) {
      if (count < lim1) res.add(factors.last);
      count += 1;
      if (count == lim1) {
        print("First $lim1 terms of a[n]:");
        print(res.take(lim1).join(', '));
        print("");
      } else if (count == lim2) {
        print("The $count term of a[n] is ${factors.last}");
        lim2 *= 10;
      }
    }
    i += 1;
  }
  print("${DateTime.now().difference(t0).inSeconds} sec.");
}

List<int> primeFactors(int n) {
  List<int> factors = [];
  while (n % 2 == 0) {
    factors.add(2);
    n ~/= 2;
  }
  for (int i = 3; i <= sqrt(n); i += 2) {
    while (n % i == 0) {
      factors.add(i);
      n ~/= i;
    }
  }
  if (n > 2) {
    factors.add(n);
  }
  return factors;
}
Output:
First 100 terms of a[n]:
1, 2, 3, 2, 5, 3, 7, 3, 5, 11, 3, 13, 7, 5, 17, 3, 19, 5, 7, 11, 23, 5, 13, 7, 29, 5, 31, 11, 17, 7, 3, 37, 19, 13, 41, 7, 43, 11, 5, 23, 47, 7, 5, 17, 13, 53, 11, 19, 29, 59, 5, 61, 31, 7, 13, 11, 67, 17, 23, 7, 71, 73, 37, 5, 19, 11, 13, 79, 41, 83, 7, 17, 43, 29, 89, 5, 13, 23, 31, 47, 19, 97, 7, 11, 5, 101, 17, 103, 7, 53, 107, 109, 11, 37, 113, 19, 23, 29, 13, 59

The 1000 term of a[n] is 109
The 10000 term of a[n] is 101
The 100000 term of a[n] is 1693
The 1000000 term of a[n] is 1202057
The 10000000 term of a[n] is 1202057

EasyLang

func maxprimcubefree num .
   t = 2
   while t * t <= num
      if num mod t = 0
         cnt += 1
         if cnt = 3
            num = 0
         .
         num = num / t
      else
         cnt = 0
         t += 1
      .
   .
   if cnt = 2 and t = num
      num = 0
   .
   return num
.
i = 1
while cnt < 100
   h = maxprimcubefree i
   if h <> 0
      cnt += 1
      write h & " "
   .
   i += 1
.
print ""
while cnt < 100000
   h = maxprimcubefree i
   if h <> 0
      cnt += 1
      if cnt = 1000 or cnt = 10000 or cnt = 100000
         print h & " "
      .
   .
   i += 1
.
Output:
1 2 3 2 5 3 7 3 5 11 3 13 7 5 17 3 19 5 7 11 23 5 13 7 29 5 31 11 17 7 3 37 19 13 41 7 43 11 5 23 47 7 5 17 13 53 11 19 29 59 5 61 31 7 13 11 67 17 23 7 71 73 37 5 19 11 13 79 41 83 7 17 43 29 89 5 13 23 31 47 19 97 7 11 5 101 17 103 7 53 107 109 11 37 113 19 23 29 13 59 
109 
101 
1693 

FreeBASIC

Translation of: Wren

Without using external libraries, it takes about 68 seconds to run on my system (Core i5).

Dim Shared As Integer factors()
Dim As Integer res(101)
res(0) = 1
Dim As Integer count = 1
Dim As Integer j, i = 2
Dim As Integer lim1 = 100
Dim As Integer lim2 = 1000
Dim As Integer max = 1e7
Dim As Integer cubeFree = 0

Sub primeFactors(n As Integer, factors() As Integer)
    Dim As Integer i = 2, cont = 2
    While (i * i <= n)
        While (n Mod i = 0)
            n /= i
            cont += 1
            Redim Preserve factors(1 To cont)
            factors(cont) = i
        Wend
        i += 1
    Wend
    If (n > 1) Then
        cont += 1
        Redim Preserve factors(1 To cont)
        factors(cont) = n
    End If
End Sub

While (count < max)
    primeFactors(i, factors())
    If (Ubound(factors) < 3) Then
        cubeFree = 1
    Else
        cubeFree = 1
        For j = 2 To Ubound(factors)
            If (factors(j-2) = factors(j-1) And factors(j-1) = factors(j)) Then
                cubeFree = 0
                Exit For
            End If
        Next j
    End If
    
    If (cubeFree = 1) Then
        If (count < lim1) Then
            res(count) = factors(Ubound(factors))
        End If
        count += 1
        If (count = lim1) Then
            Print "First "; lim1; " terms of a[n]:"
            For j = 1 To lim1
                Print Using "####"; res(j-1);
                If (j Mod 10 = 0) Then Print
            Next j
            Print
        Elseif (count = lim2) Then
            Print "The"; count; " term of a[n] is"; factors(Ubound(factors))
            lim2 *= 10
        End If
    End If
    i += 1
Wend

Sleep
Output:
First 100 terms of a[n]:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The 1000th term of a[n] is 109
The 10000th term of a[n] is 101
The 100000th term of a[n] is 1693
The 1000000th term of a[n] is 1202057
The 10000000th term of a[n] is 1202057

Java

import java.util.Arrays;
import java.util.List;
import java.util.stream.Collectors;
import java.util.stream.IntStream;
import java.util.stream.Stream;

public final class GreatestPrimeDividingTheNthCubeFreeNumber {

	public static void main(String[] args) {
		final int limit = 10_000_000;
		createPrimes(limit);
		createCubeFree((int) (limit * 1.25));
		
		System.out.println("The first 100 terms of a370833 are:");
		for ( int i = 1; i <= 100; i++ ) {
			System.out.print(String.format("%3d%s", a370833(i).factor, ( i % 10 == 0 ? "\n" : " " )));			
		}
		System.out.println();
			
		int n = 1_000;
		while ( n <= limit ) {
		    NumberAndFactor result = a370833(n);
		    System.out.println("The " + n + "th term of a370833 is " + result.factor 
		    	               + " for the cube free number " + result.number);
		    n *= 10;
		} 
	}
	
	private static NumberAndFactor a370833(int n) {
	    if ( n == 1 ) { 
	    	return new NumberAndFactor(1, 1);
	    }
	    
	    final int nth = cubeFree.get(n - 1);
	    return new NumberAndFactor(nth, greatestPrimeFactor(nth));
	}
	
	private static void createCubeFree(int n) {
	    List<Boolean> indicators = Stream.generate( () -> true ).limit(n + 1).collect(Collectors.toList());
	    
	    for ( int i = 0, prime = primes[i]; prime < Math.cbrt(n) + 1; prime = primes[++i] ) {
	        int p3 = prime * prime * prime;
	        int k = 1;
	        while ( p3 * k <= n ) {
	            indicators.set(p3 * k, false);
	            k += 1;
	        }
	    }
	     
	    cubeFree = IntStream.range(1, indicators.size()).filter( i -> indicators.get(i) ).boxed().toList();
	}
	
	private static int greatestPrimeFactor(int number) {
		int largest = 0;
		for ( int i = 0, prime = primes[i]; prime <= number; prime = primes[++i] ) {
			while ( number % prime == 0 ) {
				number /= prime;
				largest = prime;
			}
		}
		return largest;
	}
	
	private static void createPrimes(int limit) {
		final int halfLimit = ( limit + 1 ) / 2;
		boolean[] composite = new boolean[halfLimit];
		for ( int i = 1, p = 3; i < halfLimit; p += 2, i++ ) {
			if ( ! composite[i] ) {
				for ( int a = i + p; a < halfLimit; a += p ) {
					composite[a] = true;
				}
			}
		}
		int[] tempPrimes = new int[composite.length];
		tempPrimes[0] = 2;
		int primesIndex = 1;
		for ( int i = 1, p = 3; i < halfLimit; p += 2, i++ ) {
			if ( ! composite[i] ) {
				tempPrimes[primesIndex++] = p;
			}
		}
		primes = Arrays.copyOfRange(tempPrimes, 0, primesIndex);
	}
	
	private static record NumberAndFactor(int number, int factor) { }
	
	private static int[] primes;
	private static List<Integer> cubeFree;

}
Output:
The first 100 terms of a370833 are:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The 1000th term of a370833 is 109 for the cube free number 1199
The 10000th term of a370833 is 101 for the cube free number 12019
The 100000th term of a370833 is 1693 for the cube free number 120203
The 1000000th term of a370833 is 1202057 for the cube free number 1202057
The 10000000th term of a370833 is 1202057 for the cube free number 12020570

jq

Works with jq, the C implementation of jq

Works with gojq, the Go implementation of jq

# The following may be omitted if using the C implementation of jq
def _nwise($n):
  def n: if length <= $n then . else .[0:$n] , (.[$n:] | n) end;
  n;

### Generic functions
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;

# tabular print
def tprint($columns; $width):
  reduce _nwise($columns) as $row ("";
     . + ($row|map(lpad($width)) | join(" ")) + "\n" );

# like while/2 but emit the final term rather than the first one
def whilst(cond; update):
     def _whilst:
         if cond then update | (., _whilst) else empty end;
     _whilst;

## Prime factors

# Emit an array of the prime factors of 'n' in order using a wheel with basis [2, 3, 5]
# e.g. 44 | primeFactors => [2,2,11]
def primeFactors:
  def out($i): until (.n % $i != 0; .factors += [$i] | .n = ((.n/$i)|floor) );
  if . < 2 then []
  else [4, 2, 4, 2, 4, 6, 2, 6] as $inc
    | { n: .,
        factors: [] }
    | out(2)
    | out(3)
    | out(5)
    | .k = 7
    | .i = 0
    | until(.k * .k > .n;
        if .n % .k == 0
        then .factors += [.k]
        | .n = ((.n/.k)|floor)
        else .k += $inc[.i]
        | .i = ((.i + 1) % 8)
        end)
    | if .n > 1 then .factors += [ .n ] else . end
  | .factors
  end;

### Cube-free numbers
# If cubefree then emit the largest prime factor, else emit null
def cubefree:
  if . % 8 == 0 or . % 27 == 0 then false
  else  primeFactors as $factors
  | ($factors|length) as $n
  | {i: 2, cubeFree: true}
  | until (.cubeFree == false or .i >= $n;
      $factors[.i-2] as $f
      | if $f == $factors[.i-1] and $f == $factors[.i]
        then .cubeFree = false
        else .i += 1
        end)
  | if .cubeFree then $factors[-1] else null end
  end;

## The tasks
  { res:    [1],  # by convention
    count:    1,  # see the previous line
    i:        2,
    lim1:   100,
    lim2:  1000,
     max: 10000 }
  | whilst (.count <= .max;
      .emit = null
      | (.i|cubefree) as $result
      | if $result
        then .count += 1
        | if .count <= .lim1 then .res += [$result] end
        | if .count == .lim1
          then .emit = ["First \(.lim1) terms of a[n]:"]
          | .emit += [.res | tprint(10; 3)]
          elif .count == .lim2
          then .lim2 *= 10
          | .emit = ["The \(.count) term of a[n] is \($result)"]
          end
        end
      | .i += 1
      | if .i % 8 == 0 or .i % 27 == 0
        then .i += 1
        end
    )
  | select(.emit) | .emit[]
Output:
First 100 terms of a[n]:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The 1000 term of a[n] is 109
The 10000 term of a[n] is 101
The 100000 term of a[n] is 1693
The 1000000 term of a[n] is 1202057

Julia

using Formatting
using Primes
using ResumableFunctions

const MAXINMASK = 10_000_000_000 # memory on test machine could not spare a bitmask much larger than this

""" return a bitmask containing at least max_wanted cubefreenumbers """
function cubefreemask(max_wanted)
    size_wanted = Int(round(max_wanted * 1.21))
    mask = trues(size_wanted)
    p = primes(Int(floor(size_wanted^(1/3))))
    for i in p
        interval = i^3
        for j in interval:interval:size_wanted
            mask[j] = false
        end
    end
    return mask
end

""" generator for cubefree numbers up to max_wanted in number """
@resumable function nextcubefree(max_wanted = MAXINMASK)
    cfmask = cubefreemask(max_wanted)
    @yield 1
    for i in firstindex(cfmask)+1:lastindex(cfmask)
        if cfmask[i]
            @yield i
        end
    end
    @warn "past end of allowable size of sequence A370833"
end

""" various task output with OEIS sequence A370833 """
function testA370833(toprint)
    println("First 100 terms of a[n]:")
    for (i, a) in enumerate(nextcubefree())
        if i < 101
            f = factor(a).pe # only factor the ones we want to print 
            highestprimefactor = isempty(f) ? 1 : f[end][begin]
            print(rpad(highestprimefactor, 4), i % 10 == 0 ? "\n" : "")
        elseif i  toprint
            highestprimefactor = (factor(a).pe)[end][begin]
            println("\n The ", format(i, commas = true), "th term of a[n] is ",
               format(highestprimefactor, commas = true))
        end
        i >= toprint[end] && break
    end
end

testA370833(map(j -> 10^j, 3:Int(round(log10(MAXINMASK)))))
Output:
First 100 terms of a[n]:
1   2   3   2   5   3   7   3   5   11  
3   13  7   5   17  3   19  5   7   11  
23  5   13  7   29  5   31  11  17  7   
3   37  19  13  41  7   43  11  5   23  
47  7   5   17  13  53  11  19  29  59  
5   61  31  7   13  11  67  17  23  7   
71  73  37  5   19  11  13  79  41  83  
7   17  43  29  89  5   13  23  31  47  
19  97  7   11  5   101 17  103 7   53
107 109 11  37  113 19  23  29  13  59

 The 1,000th term of a[n] is 109

 The 10,000th term of a[n] is 101

 The 100,000th term of a[n] is 1,693

 The 1,000,000th term of a[n] is 1,202,057

 The 10,000,000th term of a[n] is 1,202,057

 The 100,000,000th term of a[n] is 20,743

 The 1,000,000,000th term of a[n] is 215,461

 The 10,000,000,000th term of a[n] is 1,322,977

Mathematica / Wolfram Language

ClearAll[cubeFreeQ, largestPrimeFactor, cubeFreeInts];
cubeFreeQ[n_Integer] := Max@FactorInteger[n][[All, 2]] < 3;
largestPrimeFactor[n_Integer] := Last@FactorInteger[n][[All, 1]];

cubeFreeInts = Select[Range[Ceiling[10^7 Zeta[3]]], cubeFreeQ]; 

Print["The first 100 terms are: "];
Print /@ StringJoin /@ Partition[StringPadLeft[ToString /@ largestPrimeFactor /@ cubeFreeInts[[;; 100]], 4], UpTo[10]];

Print[];

Print["The ", #, "th term is: ", largestPrimeFactor[cubeFreeInts[[#]]]] & /@ (10^Range[3, 7]);
Output:
The first 100 terms are: 
   1   2   3   2   5   3   7   3   5  11
   3  13   7   5  17   3  19   5   7  11
  23   5  13   7  29   5  31  11  17   7
   3  37  19  13  41   7  43  11   5  23
  47   7   5  17  13  53  11  19  29  59
   5  61  31   7  13  11  67  17  23   7
  71  73  37   5  19  11  13  79  41  83
   7  17  43  29  89   5  13  23  31  47
  19  97   7  11   5 101  17 103   7  53
 107 109  11  37 113  19  23  29  13  59

The 1000th term is: 109
The 10000th term is: 101
The 100000th term is: 1693
The 1000000th term is: 1202057
The 10000000th term is: 1202057

Pascal

Free Pascal

Uses sieving with cube powers of primes.
Nearly linear runtime. 0.8s per Billion. Highest Prime for 1E9 is 997

program CubeFree2;
{$IFDEF FPC}
  {$MODE DELPHI}
  {$OPTIMIZATION ON,ALL}
//{$CODEALIGN proc=16,loop=8} //TIO.RUN (Intel Xeon 2.3 Ghz) loves it 
  {$COPERATORS ON}
{$ELSE}
  {$APPTYPE CONSOLE}
{$ENDIF}
uses
  sysutils
{$IFDEF WINDOWS},Windows{$ENDIF}
  ;
const
  SizeCube235 =4* (2*2*2* 3*3*3 *5*5*5);//2*27000 <= 64kb level I
type
  tPrimeIdx = 0..65535;
  tPrimes = array[tPrimeIdx] of Uint32;
  tSv235IDx = 0..SizeCube235-1;
  tSieve235 = array[tSv235IDx] of boolean;
  tDl3 = record
               dlPr3 : UInt64;
               dlSivMod,
               dlSivNum : Uint32;
            end;
  tDelCube = array[tPrimeIdx] of tDl3;

var
  {$ALIGN 8}
  SmallPrimes: tPrimes;
  Sieve235,
  Sieve : tSieve235;
  DelCube : tDelCube;

procedure InitSmallPrimes;
//get primes. #0..65535.Sieving only odd numbers
const
  MAXLIMIT = (821641-1) shr 1;
var
  pr : array[0..MAXLIMIT] of byte;
  p,j,d,flipflop :NativeUInt;
Begin
  SmallPrimes[0] := 2;
  fillchar(pr[0],SizeOf(pr),#0);
  p := 0;
  repeat
    repeat
      p +=1
    until pr[p]= 0;
    j := (p+1)*p*2;
    if j>MAXLIMIT then
      BREAK;
    d := 2*p+1;
    repeat
      pr[j] := 1;
      j += d;
    until j>MAXLIMIT;
  until false;

  SmallPrimes[1] := 3;
  SmallPrimes[2] := 5;
  j := 3;
  d := 7;
  flipflop := (2+1)-1;//7+2*2,11+2*1,13,17,19,23
  p := 3;
  repeat
    if pr[p] = 0 then
    begin
      SmallPrimes[j] := d;
      inc(j);
    end;
    d += 2*flipflop;
    p+=flipflop;
    flipflop := 3-flipflop;
  until (p > MAXLIMIT) OR (j>High(SmallPrimes));
end;

procedure Init235(var Sieve235:tSieve235);
var
  i,j,k : NativeInt;
begin
  fillchar(Sieve235,SizeOf(Sieve235),Ord(true));
  Sieve235[0] := false;
  for k in [2,3,5] do
  Begin
    j := k*k*k;
    i := j;
    while i < SizeCube235 do
    begin
      Sieve235[i] := false;
      inc(i,j);
    end;
  end;
end;

procedure InitDelCube(var DC:tDelCube);
var
  i,q,r : Uint64;
begin
  for i in tPrimeIdx do
  begin
    q := SmallPrimes[i];
    q *= sqr(q);
    with DC[i] do
    begin
      dlPr3 := q;
      r := q div SizeCube235;
      dlSivNum := r;
      dlSivMod := q-r*SizeCube235;
    end;
  end;
end;
function Numb2USA(n:Uint64):Ansistring;
var
  pI :pChar;
  i,j : NativeInt;
Begin
  str(n,result);
  i := length(result);
 //extend s by the count of comma to be inserted
  j := i+ (i-1) div 3;
  if i<> j then
  Begin
    setlength(result,j);
    pI := @result[1];
    dec(pI);
    while i > 3 do
    Begin
       //copy 3 digits
       pI[j] := pI[i];pI[j-1] := pI[i-1];pI[j-2] := pI[i-2];
       // insert comma
       pI[j-3] := ',';
       dec(i,3);
       dec(j,4);
    end;
  end;
end;

function highestDiv(n: uint64):Uint64;
//can test til 821641^2 ~ 6,75E11
var
  pr : Uint64;
  i : integer;
begin
  result := n;
  i := 0;
  repeat
    pr := Smallprimes[i];
    if pr*pr>result then
      BREAK;
    while (result > pr) AND (result MOD pr = 0) do
      result := result DIV pr;
    inc(i);
  until i > High(SmallPrimes);
end;

procedure OutNum(lmt,n:Uint64);
begin
  writeln(Numb2Usa(lmt):18,Numb2Usa(n):19,Numb2Usa(highestDiv(n)):18);
end;


var
  sieveNr,minIdx,maxIdx : Uint32;
procedure SieveOneSieve;
var
  j : Uint64;
  i : Uint32;
begin
 // sieve with previous primes
  Sieve := Sieve235;
  For i := minIdx to MaxIdx do
    with DelCube[i] do
      if dlSivNum = sievenr then
      begin
        j :=  dlSivMod;
        repeat
          sieve[j] := false;
          inc(j,dlPr3);
        until j >= SizeCube235;
        dlSivMod := j Mod SizeCube235;
        dlSivNum += j div SizeCube235;
      end;
  //sieve with new primes
  while DelCube[maxIdx+1].dlSivNum = sieveNr do
  begin
    inc(maxIdx);
    with DelCube[maxIdx] do
    begin
      j :=  dlSivMod;
      repeat
        sieve[j] := false;
        inc(j,dlPr3);
      until j >= SizeCube235;
      dlSivMod := j Mod SizeCube235;
      dlSivNum := sieveNr + j div SizeCube235;
    end;
  end;
end;

var
  T0:Int64;
  cnt,lmt : Uint64;
  i : integer;

Begin
  T0 := GetTickCount64;
  InitSmallPrimes;
  Init235(Sieve235);
  InitDelCube(DelCube);

  sieveNr := 0;
  minIdx := low(tPrimeIdx);
  while Smallprimes[minIdx]<=5 do
    inc(minIdx);
  MaxIdx := minIdx;
  while DelCube[maxIdx].dlSivNum <= sieveNr do
    inc(maxIdx);

  SieveOneSieve;
  i := 1;
  cnt := 0;
  lmt := 100;
  repeat
    if sieve[i] then
    begin
      inc(cnt);
      write(highestDiv(i):4);
      if cnt mod 10 = 0 then
        Writeln;
    end;
    inc(i);
  until cnt = lmt;
  Writeln;

  cnt := 0;
  lmt *=10;
  repeat
    For i in tSv235IDx do
    begin
      if sieve[i] then
      begin
        inc(cnt);
        if cnt = lmt then
        Begin
          OutNum(lmt,i+sieveNr*SizeCube235);
          lmt*= 10;
        end;
      end;
    end;
    inc(sieveNr);
    SieveOneSieve;
  until lmt > 1*1000*1000*1000;

  T0 := GetTickCount64-T0;
  writeln('runtime ',T0/1000:0:3,' s');
end.
@home:
   1   2   3   2   5   3   7   3   5  11
   3  13   7   5  17   3  19   5   7  11
  23   5  13   7  29   5  31  11  17   7
   3  37  19  13  41   7  43  11   5  23
  47   7   5  17  13  53  11  19  29  59
   5  61  31   7  13  11  67  17  23   7
  71  73  37   5  19  11  13  79  41  83
   7  17  43  29  89   5  13  23  31  47
  19  97   7  11   5 101  17 103   7  53
 107 109  11  37 113  19  23  29  13  59

             1,000              1,199               109
            10,000             12,019               101
           100,000            120,203             1,693
         1,000,000          1,202,057         1,202,057
        10,000,000         12,020,570         1,202,057
       100,000,000        120,205,685            20,743
     1,000,000,000      1,202,056,919           215,461 // TIO.RUN 2.4 s 
    10,000,000,000     12,020,569,022         1,322,977
   100,000,000,000    120,205,690,298           145,823
 1,000,000,000,000  1,202,056,903,137   400,685,634,379
runtime 805.139 s
real    13m25,140s

recursive alternative

Using Apéry's Constant, which is a quite good estimate.
Only checking powers of 10.Not willing to test prime factors > 2,642,246-> 0

program CubeFree3;
{$IFDEF FPC}
  {$MODE DELPHI}{$OPTIMIZATION ON,ALL}  {$COPERATORS ON}
  {$CODEALIGN proc=16,loop=8} //TIO.RUN loves loop=8
{$ELSE}
  {$APPTYPE CONSOLE}
{$ENDIF}
uses
  sysutils
{$IFDEF WINDOWS},Windows{$ENDIF}
  ;
const
   //Apéry's Constant
  Z3 : extended  = 1.20205690315959428539973816151144999;
  RezZ3 = 0.831907372580707468683126278821530734417;

type
  tPrimeIdx = 0..192724;// primes til 2,642,246 ^3 ~ 2^64-1= High(Uint64)
  tPrimes = array[tPrimeIdx] of Uint32;
  tDl3 = UInt64;
  tPrmCubed = array[tPrimeIdx] of tDl3;

var
  SmallPrimes: tPrimes;
  {$ALIGN 32}
  PrmCubed : tPrmCubed;

procedure InitSmallPrimes;
//get primes. #0..192724.Sieving only odd numbers
const
  MAXLIMIT = (2642246-1) shr 1;
var
  pr : array[0..MAXLIMIT] of byte;
  p,j,d,flipflop :NativeUInt;
Begin
  SmallPrimes[0] := 2;
  fillchar(pr[0],SizeOf(pr),#0);
  p := 0;
  repeat
    repeat
      p +=1
    until pr[p]= 0;
    j := (p+1)*p*2;
    if j>MAXLIMIT then
      BREAK;
    d := 2*p+1;
    repeat
      pr[j] := 1;
      j += d;
    until j>MAXLIMIT;
  until false;

  SmallPrimes[1] := 3;
  SmallPrimes[2] := 5;
  j := 3;
  d := 7;
  flipflop := (2+1)-1;//7+2*2,11+2*1,13,17,19,23
  p := 3;
  repeat
    if pr[p] = 0 then
    begin
      SmallPrimes[j] := d;
      inc(j);
    end;
    d += 2*flipflop;
    p+=flipflop;
    flipflop := 3-flipflop;
  until (p > MAXLIMIT) OR (j>High(SmallPrimes));
end;

procedure InitPrmCubed(var DC:tPrmCubed);
var
  i,q : Uint64;
begin
  for i in tPrimeIdx do
  begin
    q := SmallPrimes[i];
    q *= sqr(q);
    DC[i] := q;
  end;
end;

function Numb2USA(n:Uint64):Ansistring;
var
  pI :pChar;
  i,j : NativeInt;
Begin
  str(n,result);
  i := length(result);
 //extend s by the count of comma to be inserted
  j := i+ (i-1) div 3;
  if i<> j then
  Begin
    setlength(result,j);
    pI := @result[1];
    dec(pI);
    while i > 3 do
    Begin
       //copy 3 digits
       pI[j] := pI[i];pI[j-1] := pI[i-1];pI[j-2] := pI[i-2];
       // insert comma
       pI[j-3] := ',';
       dec(i,3);
       dec(j,4);
    end;
  end;
end;

function highestDiv(n: uint64):Uint64;
//can test til 2642246^2 ~ 6,98E12
var
  pr : Uint64;
  i : integer;
begin
  result := n;
  for i in tPrimeIdx do
  begin
    pr := Smallprimes[i];
    if pr*pr>result then
      BREAK;
    while (result > pr) AND (result MOD pr = 0) do
      result := result DIV pr;
  end;
end;

procedure OutNum(lmt,n:Uint64);
var
  MaxPrimeFac : Uint64;
begin
  MaxPrimeFac := highestDiv(lmt);
  if  MaxPrimeFac > sqr(SmallPrimes[high(tPrimeIdx)]) then
    MaxPrimeFac := 0;
  writeln(Numb2Usa(lmt):26,'|',Numb2Usa(n):26,'|',Numb2Usa(MaxPrimeFac):15);
end;
//##########################################################################
var
  cnt : Int64;

procedure check(lmt:Uint64;i:integer;flip :Boolean);
var
  p : Uint64;
begin
  For i := i to high(tPrimeIdx) do
  begin
    p := PrmCubed[i];
    if lmt < p then
      BREAK;
    p := lmt DIV p;
    if flip then
      cnt -= p
    else
      cnt += p;
    if p >= PrmCubed[i+1] then
      check(p,i+1,not(flip));
  end;
end;

function GetLmtfromCnt(inCnt:Uint64):Uint64;
begin
  result := trunc(Z3*inCnt);
  repeat
    cnt := result;
    check(result,0,true);
    //new approximation
    inc(result,trunc(Z3*(inCnt-Cnt)));
  until cnt = inCnt;
  //maybe lmt is not cubefree, like 1200 for cnt 1000
  //faster than checking for cubefree of lmt for big lmt
  repeat
    dec(result);
    cnt := result;
    check(result,0,true);
  until cnt < inCnt;
  inc(result);
end;
//##########################################################################

var
  T0,lmt:Int64;
  i : integer;
Begin
  InitSmallPrimes;
  InitPrmCubed(PrmCubed);

  For i := 1 to 100 do
  Begin
    lmt := GetLmtfromCnt(i);
    write(highestDiv(lmt):4);
    if i mod 10 = 0 then
      Writeln;
  end;
  Writeln;

  Writeln('Tested with Apéry´s Constant approximation of ',Z3:17:15);
  write('                   ');
  writeln('Limit  |       cube free numbers  |max prim factor');
  T0 := GetTickCount64;
  lmt := 1;
  For i := 0 to 18 do
  Begin
    OutNum(GetLmtfromCnt(lmt),lmt);
    lmt *= 10;
  end;
  T0 := GetTickCount64-T0;
  writeln(' runtime ',T0/1000:0:3,' s');

end.
@home:
   1   2   3   2   5   3   7   3   5  11
   3  13   7   5  17   3  19   5   7  11
  23   5  13   7  29   5  31  11  17   7
   3  37  19  13  41   7  43  11   5  23
  47   7   5  17  13  53  11  19  29  59
   5  61  31   7  13  11  67  17  23   7
  71  73  37   5  19  11  13  79  41  83
   7  17  43  29  89   5  13  23  31  47
  19  97   7  11   5 101  17 103   7  53
 107 109  11  37 113  19  23  29  13  59

Tested with Apéry´s Constant approximation of 1.202056903159594
                   Limit  |       cube free numbers  |max prim factor|   divs
                         1|                         1|              1|      0
                        11|                        10|             11|      1
                       118|                       100|             59|      2
                     1,199|                     1,000|            109|      6<
                    12,019|                    10,000|            101|     14
                   120,203|                   100,000|          1,693|     30
                 1,202,057|                 1,000,000|      1,202,057|     65<
                12,020,570|                10,000,000|      1,202,057|    141
               120,205,685|               100,000,000|         20,743|    301
             1,202,056,919|             1,000,000,000|        215,461|    645<
            12,020,569,022|            10,000,000,000|      1,322,977|  1,392
           120,205,690,298|           100,000,000,000|        145,823|  3,003
         1,202,056,903,137|         1,000,000,000,000|400,685,634,379|  6,465<
        12,020,569,031,641|        10,000,000,000,000|      1,498,751| 13,924
       120,205,690,315,927|       100,000,000,000,000|         57,349| 30,006
     1,202,056,903,159,489|     1,000,000,000,000,000| 74,509,198,733| 64,643<
    12,020,569,031,596,003|    10,000,000,000,000,000|              0|139,261
   120,205,690,315,959,316|   100,000,000,000,000,000|              0|300,023
 1,202,056,903,159,593,905| 1,000,000,000,000,000,000|         89,387|646,394<
 runtime 0.008 s //best value.
real    0m0,013s
Tested with Apéry´s Constant approximation of   1.000000000000000
 runtime 0.065 s
real    0m0,071s

Perl

Translation of: Python
# 20240918 Perl programming solution

use strict;
use warnings;
use Math::Prime::Util qw(factor);

my ($count, $i, $lim1, $lim2, $max, @result) = (1, 2, 100, 1000, 1e7, (1));

while ($count < $max) {
   my ($cubeFree, @factors) = (1, factor($i));
   if (scalar @factors >= 3) {
      for my $j (2..$#factors) {
         if ($factors[$j-2]==$factors[$j-1] && $factors[$j-1]==$factors[$j]) {
            $cubeFree = 0;
            last;
         }
      }
   }
   if ($cubeFree) {
      if ($count < $lim1) { push @result, $factors[-1] }
      $count++;
      if ($count == $lim1) {
         print "The first terms of A370833 are:\n";
         for (my $k = 0; $k < @result; $k += 10) {
            printf "%3d " x 10 . "\n", @result[$k..$k+9];
         }
         print "\n";
      }
      if ($count == $lim2) {
         printf "The %8dth term of A370833 is %7d\n", $count, $factors[-1];
         $lim2 *= 10;
      }
   }
   $i++;
}
Output:
The first terms of A370833 are:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The     1000th term of A370833 is     109
The    10000th term of A370833 is     101
The   100000th term of A370833 is    1693
The  1000000th term of A370833 is 1202057
The 10000000th term of A370833 is 1202057

Phix

Completes to 1e9 in the blink of an eye, and 93s for 1e11. Admittedly 1e12 is well out of reach, and this is far from the best way to get the full series.

with javascript_semantics
function cubes_before(atom n)
    -- nb: if n is /not/ cube-free it /is/ included in the result.
    -- nth := n-cubes_before(n) means n is the nth cube-free integer,
    --                but only if cubes_before(n-1)==cubes_before(n),
    -- otherwise cubes_before(cubicate) isn't really very meaningful.
    atom r = 0
    bool xpm = true -- extend prior multiples
    sequence pm = {}
    for i=1 to floor(power(n,1/3)) do
        atom p3 = power(get_prime(i),3)
        if p3>n then exit end if
        integer k = floor(n/p3)
        for mask=1 to power(2,length(pm))-1 do
            integer m = mask, mi = 0, bc = count_bits(mask)
            atom kpm = p3  
            while m do  
                mi += 1
                if odd(m) then
                    kpm *= pm[mi]
                end if
                m = floor(m/2)
            end while
            if kpm>n then
                if bc=1 then
                    xpm = false
                    pm = pm[1..mi-1]
                    exit
                end if
            else
                -- account for already counted multiples.
                integer l = floor(n/kpm)
                -- -pairs +triples -quads ... as per link above
                if odd(bc) then
                    k -= l
                else
                    k += l
                end if
            end if
        end for
        r += k
        if xpm then
            pm &= p3
        end if
    end for
    return r
end function

function cube_free(atom nth)
    -- get the nth cube-free integer
    atom lo = nth, hi = lo*2, mid, cb, k
    while hi-cubes_before(hi)<nth do
        lo = hi
        hi = lo*2
    end while
    -- bisect until we have a possible...
    atom t1 = time()+1
    while true do
        mid = floor((lo+hi)/2)
        cb = cubes_before(mid)
        k = mid-cb
        if k=nth then
            -- skip omissions
            while cubes_before(mid-1)!=cb do
                mid -= 1
                cb -= 1
            end while
            exit
        elsif k<nth then
            lo = mid
        else
            hi = mid
        end if
        if time()>t1 then
            progress("bisecting %,d..%,d (diff %,d)...\r",{lo,hi,hi-lo})
            t1 = time()+1
        end if
    end while 
    return mid
end function

function A370833(atom nth)
    if nth=1 then return {1,1,1} end if
    atom n = cube_free(nth)
    sequence f = prime_powers(n)
    return {nth,n,f[$][1]}
end function

atom t0 = time()
sequence f100 = vslice(apply(tagset(100),A370833),3)
printf(1,"First 100 terms of a[n]:\n%s\n",join_by(f100,1,10," ",fmt:="%3d"))
for n in sq_power(10,tagset(11,3)) do
    printf(1,"The %,dth term of a[n] is %,d with highest divisor %,d.\n",A370833(n))
end for
?elapsed(time()-t0)
Output:
First 100 terms of a[n]:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The 1,000th term of a[n] is 1,199 with highest divisor 109.
The 10,000th term of a[n] is 12,019 with highest divisor 101.
The 100,000th term of a[n] is 120,203 with highest divisor 1,693.
The 1,000,000th term of a[n] is 1,202,057 with highest divisor 1,202,057.
The 10,000,000th term of a[n] is 12,020,570 with highest divisor 1,202,057.
The 100,000,000th term of a[n] is 120,205,685 with highest divisor 20,743.
The 1,000,000,000th term of a[n] is 1,202,056,919 with highest divisor 215,461.
The 10,000,000,000th term of a[n] is 12,020,569,022 with highest divisor 1,322,977.
The 100,000,000,000th term of a[n] is 120,205,690,298 with highest divisor 145,823.
"1 minute and 33s"

Python

Library: SymPy
Works with: Python version 3.x
Translation of: FreeBASIC

This takes under 12 minutes to run on my system (Core i5).

#!/usr/bin/python

from sympy import primefactors

res = [1]
count = 1
i = 2
lim1 = 100
lim2 = 1000
max = 1e7

while count < max:
    cubeFree = False
    factors = primefactors(i)
    if len(factors) < 3:
        cubeFree = True
    else:
        cubeFree = True
        for j in range(2, len(factors)):
            if factors[j-2] == factors[j-1] and factors[j-1] == factors[j]:
                cubeFree = False
                break
    if cubeFree:
        if count < lim1:
            res.append(factors[-1])
        count += 1
        if count == lim1:
            print("First {} terms of a[n]:".format(lim1))
            for k in range(0, len(res), 10):
                print(" ".join(map(str, res[k:k+10])))
            print("")
        elif count == lim2:
            print("The {} term of a[n] is {}".format(count, factors[-1]))
            lim2 *= 10
    i += 1
Output:
Similar to FreeBASIC entry.

Raku

use Prime::Factor;

sub max_factor_if_cubefree ($i) {
    my @f = prime-factors($i);
    return @f.tail if @f.elems          < 3
                   or @f.Bag.values.all < 3;
}

constant @Aₙ = lazy flat 1, map &max_factor_if_cubefree, 2..*;

say 'The first terms of A370833 are:';
say .fmt('%3d') for @Aₙ.head(100).batch(10);

say '';

for 10 X** (3..6) -> $k {
    printf "The %8dth term of A370833 is %7d\n", $k, @Aₙ[$k-1];
}
Output:
The first terms of A370833 are:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The     1000th term of A370833 is     109
The    10000th term of A370833 is     101
The   100000th term of A370833 is    1693
The  1000000th term of A370833 is 1202057

RPL

I confirm it does take a while for an interpreted language like RPL. Getting the 100,000th term is likely to be a question of hours, even on an emulator.

Works with: HP version 49
« 1 { } DUP2 + → n res2 res1
  « 2
    WHILE 'n' INCR 10000 ≤ REPEAT
      WHILE DUP FACTORS DUP 1 « 3 < NSUB 2 MOD NOT OR » DOSUBS ΠLIST NOT 
      REPEAT DROP 1 + END
      HEAD
      CASE 
         n 100 ≤      THEN 'res1' OVER STO+ END
         n LOG FP NOT THEN 'res2' OVER STO+ END
      END
      DROP 1 +
    END
    DROP res1 res2
» » 'TASK' STO 
Output:
2: { 1 2 3 2 5 3 7 3 5 11 3 13 7 5 17 3 19 5 7 11 23 5 13 7 29 5 31 11 17 7 3 37 19 13 41 7 43 11 5 23 47 7 5 17 13 53 11 19 29 59 5 61 31 7 13 11 67 17 23 7 71 73 37 5 19 11 13 79 41 83 7 17 43 29 89 5 13 23 31 47 19 97 7 11 5 101 17 103 7 53 107 109 11 37 113 19 23 29 13 59 }
1: { 109 101 }

Sidef

Version 1

Built-in as Number#nth_cubefree():

say "The first 100 terms of A370833 are:"
cubefree(100.nth_cubefree).map{.gpf}.each_slice(10, {|*a|
    a.map { '%3s' % _ }.join(' ').say
}); say ''

for n in (1..20) {
    say ("a(10^#{n}) = ", nth_cubefree(10**n).gpf)
}
Output:
The first 100 terms of A370833 are:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

a(10^1) = 11
a(10^2) = 59
a(10^3) = 109
a(10^4) = 101
a(10^5) = 1693
a(10^6) = 1202057
a(10^7) = 1202057
a(10^8) = 20743
a(10^9) = 215461
a(10^10) = 1322977
a(10^11) = 145823
a(10^12) = 400685634379
a(10^13) = 1498751
a(10^14) = 57349
a(10^15) = 74509198733
a(10^16) = 632661527978737
a(10^17) = 3018045307
a(10^18) = 89387
a(10^19) = 14563696231181
a(10^20) = 13237054323968663

Version 2

Explicit implementation:

func powerfree_count(n, k=2) {
    var sum = 0
    n.iroot(k).each_squarefree {|v|
        sum += (moebius(v) * idiv(n, v**k))
    }
    return sum
}

func nth_powerfree(n,k=2) {

    n <= 0 && return NaN
    n == 1 && return 1

    var v = int(zeta(k)*n)
    var c = powerfree_count(v, k)

    while (!v.is_powerfree(k)) {
        --v
    }

    while (c != n) {
        var j = (n <=> c)
        v += j
        c += j
        v += j while !v.is_powerfree(k)
    }

    return v
}

say {|n| nth_powerfree(n, 3).gpf }.map(1..100)

for n in (1..20) {
    say ("a(10^#{n}) = ", nth_powerfree(10**n, 3).gpf)
}

Try it online!

Wren

Library: Wren-math
Library: Wren-fmt

Version 1

Simple brute force approach though skipping numbers which are multiples of 8 or 27 as these can't be cubefree.

Runtime about 3 minutes 36 seconds on my system (Core i7).

import "./math" for Int
import "./fmt" for Fmt

var res = [1]
var count = 1
var i = 2
var lim1 = 100
var lim2 = 1000
var max = 1e7
while (count < max) {
    var cubeFree = false
    var factors = Int.primeFactors(i)
    if (factors.count < 3) {
        cubeFree = true
    } else {
        cubeFree = true
        for (i in 2...factors.count) {
            if (factors[i-2] == factors[i-1] && factors[i-1] == factors[i]) {
                cubeFree = false
                break
            }
        }
    }
    if (cubeFree) {
        if (count < lim1) res.add(factors[-1])
        count = count + 1
        if (count == lim1) {
            System.print("First %(lim1) terms of a[n]:")
            Fmt.tprint("$3d", res, 10)
        } else if (count == lim2) {
            Fmt.print("\nThe $,r term of a[n] is $,d.", count, factors[-1])
            lim2 = lim2 * 10
        }
    }
    i = i + 1
    if (i % 8 == 0 || i % 27 == 0) i = i + 1
}
Output:
First 100 terms of a[n]:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The 1,000th term of a[n] is 109.

The 10,000th term of a[n] is 101.

The 100,000th term of a[n] is 1,693.

The 1,000,000th term of a[n] is 1,202,057.

The 10,000,000th term of a[n] is 1,202,057.

Version 2

This uses a simple sieve to find cubefree numbers up to a given limit which means we only now need to factorize the numbers of interest to find the greatest prime factor.

The 10 millionth term is now found in 0.85 seconds and the 1 billionth in about 94 seconds.

However, a lot of memory is needed for the sieve since all values in Wren (including bools) need 8 bytes of storage each.

We could use only 1/32nd as much memory by importing the BitArray class from Wren-array (see program comments for changes needed). However, unfortunately this is much slower to index than a normal List of booleans and the 10 millionth term would now take just over 2 seconds to find and the 1 billionth just under 4 minutes.

import "./math" for Int
//import "./array" for BitArray
import "./fmt" for Fmt

var cubeFreeSieve = Fn.new { |n|
    var cubeFree = List.filled(n+1, true) // or BitArray.new(n+1, true)
    var primes = Int.primeSieve(n.cbrt.ceil)
    for (p in primes) {
        var p3 = p * p * p
        var k = 1
        while (p3 * k <= n) {
            cubeFree[p3 * k] = false
            k = k + 1
        }
    }
    return cubeFree
}

var al = [1]
var count = 1
var i = 2
var lim1 = 100
var lim2 = 1000
var max = 1e9
var cubeFree = cubeFreeSieve.call(max * 1.25)
while (count < max) {
    if (cubeFree[i]) {
        count = count + 1
        if (count <= lim1) {
            var factors = Int.primeFactors(i)
            al.add(factors[-1])
            if (count == lim1) {
                System.print("First %(lim1) terms of a[n]:")
                Fmt.tprint("$3d", al, 10)
            }
        } else if (count == lim2) {
            var factors = Int.primeFactors(i)
            Fmt.print("\nThe $,r term of a[n] is $,d.", count, factors[-1])
            lim2 = lim2 * 10
        }
    }
    i = i + 1
}
Output:
First 100 terms of a[n]:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The 1,000th term of a[n] is 109.

The 10,000th term of a[n] is 101.

The 100,000th term of a[n] is 1,693.

The 1,000,000th term of a[n] is 1,202,057.

The 10,000,000th term of a[n] is 1,202,057.

The 100,000,000th term of a[n] is 20,743.

The 1,000,000,000th term of a[n] is 215,461.

Version 3

Translation of: Phix
Library: Wren-iterate

Even greater speed-up, now taking only 0.04 seconds to reach 10 million and 36.5 seconds to reach 10 billion.

import "./math" for Int
import "./fmt" for Conv, Fmt
import "./iterate" for Stepped

var start = System.clock

var primes = Int.primeSieve(3000)

var countBits = Fn.new { |n| Conv.bin(n).count { |c| c == "1" } }

var cubesBefore = Fn.new { |n|
    var r = 0
    var xpm = true
    var pm = []
    for (i in 1..n.cbrt.floor) {
        var p3 = primes[i-1].pow(3)
        if (p3 > n) break
        var k = (n/p3).floor
        for (mask in Stepped.ascend(1..2.pow(pm.count)-1)) {
            var m = mask
            var mi = 0
            var bc = countBits.call(mask)
            var kpm = p3
            while (m != 0) {
                mi = mi + 1
                if (m % 2 == 1) kpm = kpm * pm[mi-1]
                m = (m/2).floor
            }
            if (kpm > n) {
                if (bc == 1) {
                    xpm = false
                    pm = pm[0...mi-1]
                    break
                }
            } else {
                var l = (n/kpm).floor
                if (bc % 2 == 1) {
                    k = k - l
                } else {
                    k = k + l
                }
            }
        }
        r = r + k
        if (xpm) pm.add(p3)
    }
    return r
}

var cubeFree = Fn.new { |nth|
    var lo = nth
    var hi = lo * 2
    var mid
    var cb
    var k
    while (hi - cubesBefore.call(hi) < nth) {
        lo = hi
        hi = lo * 2
    }
    while (true) {
        mid = ((lo + hi)/2).floor
        cb = cubesBefore.call(mid)
        k = mid - cb
        if (k == nth) {
            while (cubesBefore.call(mid-1) != cb) {
                mid = mid - 1
                cb = cb - 1
            }
            break
        } else if (k < nth) {
            lo = mid
        } else {
            hi = mid
        }
    }
    return mid
}

var a370833 = Fn.new { |n|
    if (n == 1) return [1, 1]
    var nth = cubeFree.call(n)
    return [Int.primeFactors(nth)[-1], nth]
}

System.print("First 100 terms of a[n]:")
Fmt.tprint("$3d", (1..100).map { |i| a370833.call(i)[0] }, 10)
System.print()
var n = 1000
while (n <= 1e10) {
    var res = a370833.call(n)
    Fmt.print("The $,r term of a[n] is $,d for cubefree number $,d.", n, res[0], res[1])
    n = n * 10
}
System.print("\nTook %(System.clock - start) seconds.")
Output:
First 100 terms of a[n]:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The 1,000th term of a[n] is 109 for cubefree number 1,199.
The 10,000th term of a[n] is 101 for cubefree number 12,019.
The 100,000th term of a[n] is 1,693 for cubefree number 120,203.
The 1,000,000th term of a[n] is 1,202,057 for cubefree number 1,202,057.
The 10,000,000th term of a[n] is 1,202,057 for cubefree number 12,020,570.
The 100,000,000th term of a[n] is 20,743 for cubefree number 120,205,685.
The 1,000,000,000th term of a[n] is 215,461 for cubefree number 1,202,056,919.
The 10,000,000,000th term of a[n] is 1,322,977 for cubefree number 12,020,569,022.

Took 36.458647 seconds.

XPL0

Simple brute force takes 43 seconds on Raspberry Pi 4.

func    MaxFactor(N);           \Return maximum factor of N that's cube-free
int     N, Div, Count, Max;
[Max:= N;  Div:= 2;
while N >= Div*Div do
        [Count:= 0;
        while rem(N/Div) = 0 do
                [Count:= Count+1;
                if Count = 3 then return 0;
                Max:= Div;
                N:= N/Div;
                ];
        Div:= Div+1;
        ];
if N > 1 then Max:= N;          \prime
return Max;
];

int     I, N, Pow10, Fac;
[Format(4, 0);
I:= 1;  N:= 0;  Pow10:= 1000;
loop    [Fac:= MaxFactor(I);
        if Fac # 0 then
                [N:= N+1;
                if N <= 100 then
                        [RlOut(0, float(Fac));
                        if rem(N/10) = 0 then CrLf(0);
                        ]
                else if N = Pow10 then
                        [IntOut(0, Pow10);  Text(0, "th term of a[n] is ");
                        IntOut(0, Fac);  CrLf(0);
                        if Pow10 = 10_000_000 then quit;
                        Pow10:= Pow10*10;
                        ];
                ];
        I:= I+1;
        ];
]
Output:
   1   2   3   2   5   3   7   3   5  11
   3  13   7   5  17   3  19   5   7  11
  23   5  13   7  29   5  31  11  17   7
   3  37  19  13  41   7  43  11   5  23
  47   7   5  17  13  53  11  19  29  59
   5  61  31   7  13  11  67  17  23   7
  71  73  37   5  19  11  13  79  41  83
   7  17  43  29  89   5  13  23  31  47
  19  97   7  11   5 101  17 103   7  53
 107 109  11  37 113  19  23  29  13  59
1000th term of a[n] is 109
10000th term of a[n] is 101
100000th term of a[n] is 1693
1000000th term of a[n] is 1202057
10000000th term of a[n] is 1202057
Cookies help us deliver our services. By using our services, you agree to our use of cookies.