// illustration of the prime number theorem

PrimeNumberTheorem:=procedure(limit)
   n := 10; 
   final_prime := 7;
   true_count := 4;
   // low-precision fields for convenience of printing
   R8 := RealField(8); R4 := RealField(4);

   // the heading for the output
   print "n\tpi(n)\tli(n)\tpi(n)/li(n)\tR(n)\tpi(n)/R(n)";
   print "-"^58;

   while n le limit do
      p := NextPrime(final_prime);
      while p le n do
         true_count +:= 1;
         final_prime := p;
         p := NextPrime(final_prime);
      end while;
      li_count   := LogIntegral(R8!n);
      true_on_li := true_count / li_count;
      rie_count  :=&+[LogIntegral(Root(R8!n, k))*MoebiusMu(k)/k : k in [1..15]];
      true_on_rie:= true_count / rie_count;

      printf "%o\t%o\t%o\t%o\t\t%o\t%o\n", 
        n, true_count,
        Round(li_count), (true_on_li lt 1 select R4 else R8)!true_on_li,
        Round(rie_count), (true_on_rie lt 1 select R4 else R8)!true_on_rie;

      n *:= 10;
   end while;
end procedure;

PrimeNumberTheorem(1000000);