Lagrange's theorem tells us that every positive integer can be written as a sum of at most four squares.

Task
Forbidden numbers
You are encouraged to solve this task according to the task description, using any language you may know.

Many, (square numbers) can be written as a "sum" of a single square. E.G.: 4 == 22.

Some numbers require at least two squares to be summed. 5 == 22 + 12

Others require at least three squares. 6 == 22 + 12 + 12

Finally, some, (the focus of this task) require the sum at least four squares. 7 == 22 + 12 + 12 + 12. There is no way to reach 7 summing fewer than four squares.

These numbers show up in crystallography; x-ray diffraction patterns of cubic crystals depend on a length squared plus a width squared plus a height squared. Some configurations that require at least four squares are impossible to index and are colloquially known as forbidden numbers.


Note that some numbers can be made from the sum of four squares: 16 == 22 + 22 + 22 + 22, but since it can also be formed from fewer than four, 16 == 42, it does not count as a forbidden number.


Task
  • Find and show the first fifty forbidden numbers.
  • Find and show the count of forbidden numbers up to 500, 5,000.


Stretch
  • Find and show the count of forbidden numbers up to 50,000, 500,000.


See also


ALGOL 68

Uses the algorithm from the OEIS page, as used by the Wren, Python, etc. samples.

BEGIN # find some forbidden numbers: numbers that cannot be formed by        #
      # summing fewer than four squares                                      #
    # returns TRUE if n is a Forbidden numbr, FALSE otherwise                #
    #         based on the Wren version of the example at oeis.org/A004215   #
    PROC is forbidden = ( INT n )BOOL:
         BEGIN
             INT m  := n;
             INT p4 := 1;
             WHILE m > 1 AND m MOD 4 = 0 DO
                 m OVERAB 4;
                 p4   *:= 4
             OD;
             ( n OVER p4 ) MOD 8 = 7
         END # is forbidden # ;
    # show the first 50 forbidden numbers and counts of Forbidden numbers    #
    INT f count      :=   0;
    INT next to show := 500;
    FOR i TO 50 000 000 DO
        IF is forbidden( i ) THEN
            f count +:= 1;
            IF f count <= 50 THEN
                print( ( " ", whole( i, -4 ) ) );
                IF f count MOD 10 = 0 THEN print( ( newline ) ) FI
            FI
        FI;
        IF i = next to show THEN
            print( ( "There are ", whole( f count, -8 )
                   , " Forbidden numbers up to ", whole( i, 0 )
                   , newline
                   )
                 );
            next to show *:= 10
        FI
    OD
END
Output:
    7   15   23   28   31   39   47   55   60   63
   71   79   87   92   95  103  111  112  119  124
  127  135  143  151  156  159  167  175  183  188
  191  199  207  215  220  223  231  239  240  247
  252  255  263  271  279  284  287  295  303  311
There are       82 Forbidden numbers up to 500
There are      831 Forbidden numbers up to 5000
There are     8330 Forbidden numbers up to 50000
There are    83331 Forbidden numbers up to 500000
There are   833329 Forbidden numbers up to 5000000
There are  8333330 Forbidden numbers up to 50000000

AppleScript

Translation of: Phix
on isForbidden(n)
    repeat while ((n mod 4 = 0) and (n > 7))
        set n to n div 4
    end repeat
    
    return (n mod 8 = 7)
end isForbidden

on forbiddenCount(limit)
    set {counter, p4, n} to {0, 1, 7}
    repeat while (n  limit)
        set p4 to p4 * 4
        set counter to counter + (limit - n) div (p4 + p4) + 1
        set n to p4 * 7
    end repeat
    
    return counter
end forbiddenCount

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 intToText(int, separator)
    set groups to {}
    repeat while (int > 999)
        set groups's beginning to ((1000 + (int mod 1000 as integer)) as text)'s text 2 thru 4
        set int to int div 1000
    end repeat
    set groups's beginning to int as integer
    return join(groups, separator)
end intToText

on task()
    set output to {"First fifty forbidden numbers:"}
    set {counter, n, row} to {0, 0, {}}
    repeat until (counter = 50)
        set n to n + 1
        if (isForbidden(n)) then
            set counter to counter + 1
            set row's end to ("   " & n)'s text -4 thru -1
            if (counter mod 10 = 0) then
                set output's end to join(row, "")
                set row to {}
            end if
        end if
    end repeat
    set output's end to row
    repeat with target in {500, 5000, 50000, 500000, 5000000, 50000000, 500000000}
        set output's end to intToText(forbiddenCount(target), ",") & ¬
            " forbidden numbers ≤ " & intToText(target, ",")
    end repeat
    return join(output, linefeed)
end task

task()
Output:
"First fifty forbidden numbers:
   7  15  23  28  31  39  47  55  60  63
  71  79  87  92  95 103 111 112 119 124
 127 135 143 151 156 159 167 175 183 188
 191 199 207 215 220 223 231 239 240 247
 252 255 263 271 279 284 287 295 303 311

82 forbidden numbers ≤ 500
831 forbidden numbers ≤ 5,000
8,330 forbidden numbers ≤ 50,000
83,331 forbidden numbers ≤ 500,000
833,329 forbidden numbers ≤ 5,000,000
8,333,330 forbidden numbers ≤ 50,000,000
83,333,328 forbidden numbers ≤ 500,000,000"

Arturo

forbidden?: function [n][
    m: new n
    v: 0
    while -> and? m > 1 0 = m % 4 [
        'm / 4
        inc 'v
    ]
    7 = mod n / 4 ^ v 8
]

print "First 50 forbidden numbers:"
forbidden: split.every:10 select.first:50 0..∞ => forbidden?
loop forbidden 'row [
    loop row 'n -> prints pad ~"|n|" 4
    print ""
]

print ""
[target n count]: [500 0 0]
while -> target =< 5e6 [
    if forbidden? n -> inc 'count
    if n = target [
        print [count "forbidden numbers up to" target]
        'target * 10
    ]
    inc 'n
]
Output:
First 50 forbidden numbers:
   7  15  23  28  31  39  47  55  60  63
  71  79  87  92  95 103 111 112 119 124
 127 135 143 151 156 159 167 175 183 188
 191 199 207 215 220 223 231 239 240 247
 252 255 263 271 279 284 287 295 303 311

82 forbidden numbers up to 500 
831 forbidden numbers up to 5000 
8330 forbidden numbers up to 50000 
83331 forbidden numbers up to 500000 
833329 forbidden numbers up to 5000000

C

A translation of the Python code in the OEIS link. Runtime around 5.6 seconds.

#include <stdio.h>
#include <stdbool.h>
#include <math.h>
#include <locale.h>

bool isForbidden(int n) {
    int m = n, v = 0, p;
    while (m > 1 && !(m % 4)) {
        m /= 4;
        ++v;
    }
    p = (int)pow(4.0, (double)v);
    return n / p % 8 == 7;
}

int main() {
    int i = 0, count = 0, limit = 500;
    printf("The first 50 forbidden numbers are:\n");
    for ( ; count < 50; ++i) {
        if (isForbidden(i)) {
            printf("%3d ", i);
            ++count;
            if (!(count+1)%10) printf("\n");
        }
    }
    printf("\n\n");
    setlocale(LC_NUMERIC, "");
    for (i = 1, count = 0; ; ++i) {
        if (isForbidden(i)) ++count;
        if (i == limit) {
            printf("Forbidden number count <= %'11d: %'10d\n", limit, count);
            if (limit == 500000000) break;
            limit *= 10;
        }
    }
    return 0;
}
Output:
The first 50 forbidden numbers are:
  7  15  23  28  31  39  47  55  60  63  71  79  87  92  95 103 111 112 119 124 127 135 143 151 156 159 167 175 183 188 191 199 207 215 220 223 231 239 240 247 252 255 263 271 279 284 287 295 303 311 

Forbidden number count <=         500:         82
Forbidden number count <=       5,000:        831
Forbidden number count <=      50,000:      8,330
Forbidden number count <=     500,000:     83,331
Forbidden number count <=   5,000,000:    833,329
Forbidden number count <=  50,000,000:  8,333,330
Forbidden number count <= 500,000,000: 83,333,328

C++

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

bool is_forbidden(const uint32_t& number) {
	uint32_t copy_number = number;
	uint32_t power_of_4 = 0;
	while ( copy_number > 1 && copy_number % 4 == 0 ) {
		copy_number /= 4;
		power_of_4++;
	}
	return ( number / static_cast<uint32_t>(std::pow(4, power_of_4)) ) % 8 == 7;
}

int main() {
	std::vector<uint32_t> forbiddens = { };
	for ( uint32_t i = 1; i <= 500'000; ++i ) {
		if ( is_forbidden(i) ) {
			forbiddens.emplace_back(i);
		}
	}

	std::cout << "The first 50 forbidden numbers are:" << "\n";
	for ( uint32_t n = 0; n < 50; ++n ) {
		std::cout << std::setw(3) << forbiddens[n] << ( n % 10 == 9 ? "\n" : " " );
	}
	std::cout << "\n";

	for ( uint32_t limit : { 500, 5'000, 50'000, 500'000 } ) {
		const auto iter = std::lower_bound(forbiddens.begin(), forbiddens.end(), limit);
		const uint32_t count = std::distance(forbiddens.begin(), iter);
	    std::cout << "There are " << count << " forbidden number count <= " << limit << "\n";
	}
}
Output:
The first 50 forbidden numbers are:
  7  15  23  28  31  39  47  55  60  63
 71  79  87  92  95 103 111 112 119 124
127 135 143 151 156 159 167 175 183 188
191 199 207 215 220 223 231 239 240 247
252 255 263 271 279 284 287 295 303 311

There are 82 forbidden number count <= 500
There are 831 forbidden number count <= 5000
There are 8330 forbidden number count <= 50000
There are 83331 forbidden number count <= 500000

FreeBASIC

Function isForbidden (num As Uinteger) As Uinteger
    Dim As Uinteger fours = num, pow4 = 0
    While (fours > 1) And (fours Mod 4 = 0)
        fours \= 4
        pow4 += 1
    Wend
    Return (num \ 4 ^ pow4) Mod 8 = 7
End Function

Dim As Integer i = 0, cnt = 0
Print "The first 50 forbidden numbers are:"
Do
    i += 1
    If isForbidden(i) Then 
        cnt += 1
        If cnt <= 50 Then Print Using "####"; i; : If cnt Mod 10 = 0 Then Print
    End If
    If i = 500 Then Print Using !"\nForbidden number count <= #,###,###: ###,###"; i; cnt
    If i = 5e3 Or i = 5e4 Or i = 5e5 Or i = 5e6 Then Print Using "Forbidden number count <= #,###,###: ###,###"; i ; cnt
Loop Until i = 5e6

Sleep
Output:
Same as Wren entry.

Go

Library: Go-rcu

A translation of the Python code in the OEIS link. Runtime around 7 seconds.

package main

import (
    "fmt"
    "math"
    "rcu"
)

func isForbidden(n int) bool {
    m := n
    v := 0
    for m > 1 && m%4 == 0 {
        m /= 4
        v++
    }
    pow := int(math.Pow(4, float64(v)))
    return n/pow%8 == 7
}

func main() {
    forbidden := make([]int, 50)
    for i, count := 0, 0; count < 50; i++ {
        if isForbidden(i) {
            forbidden[count] = i
            count++
        }
    }
    fmt.Println("The first 50 forbidden numbers are:")
    rcu.PrintTable(forbidden, 10, 3, false)
    fmt.Println()
    limit := 500
    count := 0
    for i := 1; ; i++ {
        if isForbidden(i) {
            count++
        }
        if i == limit {
            slimit := rcu.Commatize(limit)
            scount := rcu.Commatize(count)
            fmt.Printf("Forbidden number count <= %11s: %10s\n", slimit, scount)
            if limit == 500_000_000 {
                return
            }
            limit *= 10
        }
    }
}
Output:
The first 50 forbidden numbers are:
  7  15  23  28  31  39  47  55  60  63 
 71  79  87  92  95 103 111 112 119 124 
127 135 143 151 156 159 167 175 183 188 
191 199 207 215 220 223 231 239 240 247 
252 255 263 271 279 284 287 295 303 311 

Forbidden number count <=         500:         82
Forbidden number count <=       5,000:        831
Forbidden number count <=      50,000:      8,330
Forbidden number count <=     500,000:     83,331
Forbidden number count <=   5,000,000:    833,329
Forbidden number count <=  50,000,000:  8,333,330
Forbidden number count <= 500,000,000: 83,333,328

J

For this task we can find forbidden numbers up to some limiting value y:

forbid=: {{
   s1=. *:i.1+<.%:y
   s2=. ~.(#~ y>:]),s1+/s1
   s3=. ~.(#~ y>:]),s2+/s1
   (1+i.y)-.s1,s2,s3
}}

In other words: s1 is square numbers up through y, s2 is unique sums of those squares up through y, s3 is unique sums of members of those two sequences, up through y, and our result is numbers up through y which do not appear in s1, s2 or s3.

The task then becomes:

   5 10$forbid 500
  7  15  23  28  31  39  47  55  60  63
 71  79  87  92  95 103 111 112 119 124
127 135 143 151 156 159 167 175 183 188
191 199 207 215 220 223 231 239 240 247
252 255 263 271 279 284 287 295 303 311
   #forbid 500
82
   #forbid 5000
831
   #forbid 50000
8330
   #forbid 500000
83331

Java

import java.util.List;
import java.util.stream.IntStream;

public final class ForbiddenNumbers {

	public static void main(String[] args) {
		List<Integer> forbiddens = 
            IntStream.rangeClosed(1, 500_000).filter( i -> isForbidden(i) ).boxed().toList();
		
		System.out.println("The first 50 forbidden numbers are:");
		for ( int n = 0; n < 50; n++ ) {
			System.out.print(String.format("%3d%s", forbiddens.get(n), ( n % 10 == 9 ? "\n" : " " )));
		}
		System.out.println();

		for ( int limit : List.of( 500, 5_000, 50_000, 500_000 ) ) {
		     final long count = forbiddens.stream().filter( i -> i <= limit ).count();
		     System.out.println("There are " + count + " forbidden number count <= " + limit);
		}
	}
	
	private static boolean isForbidden(int number) {
		int copyNumber = number;
	    int powerOf4 = 0;
	    while ( copyNumber > 1 && copyNumber % 4 == 0 ) {
	        copyNumber /= 4;
	        powerOf4 += 1;
	    }
	    return ( number / Math.pow(4, powerOf4) ) % 8 == 7;
	} 

}
Output:
The first 50 forbidden numbers are:
  7  15  23  28  31  39  47  55  60  63
 71  79  87  92  95 103 111 112 119 124
127 135 143 151 156 159 167 175 183 188
191 199 207 215 220 223 231 239 240 247
252 255 263 271 279 284 287 295 303 311

There are 82 forbidden number count <= 500
There are 831 forbidden number count <= 5000
There are 8330 forbidden number count <= 50000
There are 83331 forbidden number count <= 500000

jq

Works with: jq

The following also works with gojq, the Go implementation of jq, except that beyond forbidden(5000), gojq's speed and memory requirements might become a problem.

def count(s): reduce s as $x (0; .+1);

def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;

# The def of _nwise can 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;

def forbidden($max):
  def ub($a;$b):
    if $b < 0 then 0 else [$a, ($b|sqrt)] | min end;

  [false, range(1; 1 + $max)]
  | reduce range(1; 1 + ($max|sqrt)) as $i (.;
      ($i*$i) as $s1
      | .[$s1] = false
      | reduce range(1; 1 + ub($i; ($max - $s1))) as $j (.;
          ($s1 + $j*$j) as $s2
	  | .[$s2] = false
          | reduce range(1; 1 + ub($j; ($max - $s2))) as $k (.;
              .[$s2 + $k*$k] = false ) ) )
  | map( select(.) ) ;

forbidden(500) as $f
| "First fifty forbidden numbers:",
  ( $f[:50] | _nwise(10) | map(lpad(3)) | join(" ") ),
  "\nForbidden number count up to 500: \(count($f[]))",
  ((5000, 50000, 500000) | "\nForbidden number count up to \(.): \(count(forbidden(.)[])) ")
Output:
First fifty forbidden numbers:
  7  15  23  28  31  39  47  55  60  63
 71  79  87  92  95 103 111 112 119 124
127 135 143 151 156 159 167 175 183 188
191 199 207 215 220 223 231 239 240 247
252 255 263 271 279 284 287 295 303 311

Forbidden number count up to 500: 82

Forbidden number count up to 5000: 831 

Forbidden number count up to 50000: 8330

Forbidden number count up to 500000: 83331

Julia

Translation of: Python
""" true if num is a forbidden number """
function isforbidden(num)
	fours, pow4 = num, 0
	while fours > 1 && fours % 4 == 0
		fours ÷= 4
		pow4 += 1
	end
	return (num ÷ 4^pow4) % 8 == 7
end

const f500M = filter(isforbidden, 1:500_000_000)

for (idx, fbd) in enumerate(f500M[begin:begin+49])
	print(lpad(fbd, 4), idx % 10 == 0 ? '\n' : "")
end

for fbmax in [500, 5000, 50_000, 500_000, 500_000_000]
	println("\nThere are $(sum(x <= fbmax for x in f500M)) forbidden numbers <= $fbmax.")
end
Output:
   7  15  23  28  31  39  47  55  60  63
  71  79  87  92  95 103 111 112 119 124
 127 135 143 151 156 159 167 175 183 188
 191 199 207 215 220 223 231 239 240 247
 252 255 263 271 279 284 287 295 303 311

There are 82 forbidden numbers <= 500.

There are 831 forbidden numbers <= 5000.

There are 8330 forbidden numbers <= 50000.

There are 83331 forbidden numbers <= 500000.

There are 83333328 forbidden numbers <= 500000000.

Nim

Uses the algorithm from the OEIS page.

import std/[math, strformat, strutils]

const Max = 500_000

func isForbidden(num: Positive): bool =
  ## Return "true" is "n" is a forbidden number.
  var fours = num
  var pow4 = 0
  while fours > 1 and (fours and 3) == 0:
    fours = fours shr 2
    inc pow4
  result = (num div 4^pow4 and 7) == 7

iterator forbiddenNumbers(): int =
  var n = 1
  while true:
    if n.isForbidden:
      yield n
    inc n

var count = 0
var lim = 500
for n in forbiddenNumbers():
  inc count
  if count <= 50:
    stdout.write &"{n:>3}"
    stdout.write if count mod 10 == 0: '\n' else: ' '
    if count == 50: echo()
  elif n > lim:
    echo &"Numbers of forbidden numbers up to {insertSep($lim)}: {insertSep($(count - 1))}"
    lim *= 10
    if lim > Max:
      break
Output:
  7  15  23  28  31  39  47  55  60  63
 71  79  87  92  95 103 111 112 119 124
127 135 143 151 156 159 167 175 183 188
191 199 207 215 220 223 231 239 240 247
252 255 263 271 279 284 287 295 303 311

Numbers of forbidden numbers up to 500: 82
Numbers of forbidden numbers up to 5_000: 831
Numbers of forbidden numbers up to 50_000: 8_330
Numbers of forbidden numbers up to 500_000: 83_331

Pascal

Free Pascal

modified formula to calc count. Using count as gag to get next forbidden number.
No runtime.

program ForbiddenNumbers;
{$IFDEF FPC}{$MODE DELPHI}{$OPTIMIZATION ON,ALL}{$ENDIF}
{$IFDEF WINDOWS}{$APPTYPE CONSOLE}{$ENDIF}
uses
  sysutils,strutils;

function isForbidden(n:NativeUint):boolean;inline;
// no need for power or div Only shr & AND when using Uint
// n > 7 => if n <= 7 -> only 4/0 would div 4 -> no forbidden number
Begin
  while (n > 7) AND (n MOD 4 = 0) do
    n := n DIV 4;
  result := n MOD 8 = 7;
end;

function CntForbiddenTilLimit(lmt:NativeUint):NativeUint;
//forNmb = 4^i * (8*j + 7) | i,j >= 0
//forNmb = Power4 *  8*j + Power4 * 7
//forNmb =  delta* j     + n
var
  Power4,delta,n : NativeUint;
begin
 result := 0;
 power4 := 1;
 repeat
   delta := Power4*8;// j = 1
   n := Power4*7;
   if n > lmt then
     Break;
   //max j to reach limit
   inc(result,(lmt-n) DIV delta+1);
   Power4 *= 4;
 until false;
end;

var
  lmt,n,cnt: NativeUint;
BEGIN
  writeln('First fifty forbidden numbers:');
  n := 1;
  lmt := 0;
  repeat;
    if isForbidden(n) then
    Begin
      write(n:4);
      inc(lmt);
      if LMT MOD 20 = 0 THEN
        writeln;
    end;
    n +=1;
  until lmt >= 50;
  writeln;
  writeln;

  writeln('count of forbidden numbers below iterative');
  n := 1;
  cnt := 0;
  lmt := 5;
  repeat
    repeat;
      //if isForbidden(n) then cnt+=1 takes to long  100% -> 65% of time
      inc(cnt,ORD(isForbidden(n)));
      n += 1;
    until n >= lmt;
    writeln(Numb2USA(IntToStr(lmt)):30,Numb2USA(IntToStr(Cnt)):25);
    lmt *= 10;
  until lmt > 500*1000*1000;
  writeln;

  writeln('count of forbidden numbers below ');
  lmt := 5;
  repeat
    writeln(Numb2USA(IntToStr(lmt)):30,Numb2USA(IntToStr(CntForbiddenTilLimit(lmt))):25);
    lmt *= 10;
  until lmt > High(lmt) DIV 4;
END.
@TIO.RUN:

First fifty forbidden numbers:
   7  15  23  28  31  39  47  55  60  63  71  79  87  92  95 103 111 112 119 124
 127 135 143 151 156 159 167 175 183 188 191 199 207 215 220 223 231 239 240 247
 252 255 263 271 279 284 287 295 303 311

count of forbidden numbers below iterative
                             5                        0
                            50                        7
                           500                       82
                         5,000                      831
                        50,000                    8,330
                       500,000                   83,331
                     5,000,000                  833,329
                    50,000,000                8,333,330
                   500,000,000               83,333,328

count of forbidden numbers below
                             5                        0
                            50                        7
                           500                       82
                         5,000                      831
                        50,000                    8,330
                       500,000                   83,331
                     5,000,000                  833,329
                    50,000,000                8,333,330
                   500,000,000               83,333,328
                 5,000,000,000              833,333,330
                50,000,000,000            8,333,333,327
               500,000,000,000           83,333,333,328
             5,000,000,000,000          833,333,333,327
            50,000,000,000,000        8,333,333,333,327
           500,000,000,000,000       83,333,333,333,326
         5,000,000,000,000,000      833,333,333,333,327
        50,000,000,000,000,000    8,333,333,333,333,325
       500,000,000,000,000,000   83,333,333,333,333,323
Real time: 1.392 s User time: 1.343 s Sys. time: 0.042 s CPU share: 99.52 %

PascalABC.NET

Translation of: Python
##
function isforbidden(num: integer): boolean;
begin
  //true if num is a forbidden number 
  var fours := num;
  var pow4 := 0;
  while (fours > 1) and (fours mod 4 = 0) do
  begin
    fours := fours div 4;
    pow4 += 1;
  end;
  result := (num div int64(4 ** pow4)) mod 8 = 7;
end;

var f500k := (0..500_001).Where(n -> isforbidden(n));

f500k.Take(50).println; 
println;
foreach var fbmax in |500, 5000, 50_000, 500_000| do
  println('There are', f500k.Where(x -> x <= fbmax).Count, 'forbidden numbers <=', fbmax)
Output:
7 15 23 28 31 39 47 55 60 63 71 79 87 92 95 103 111 112 119 124 127 135 143 151 156 159 167 175 183 188 191 199 207 215 220 223 231 239 240 247 252 255 263 271 279 284 287 295 303 311

There are 82 forbidden numbers <= 500
There are 831 forbidden numbers <= 5000
There are 8330 forbidden numbers <= 50000
There are 83331 forbidden numbers <= 500000

Perl

use strict;
use warnings;
use List::AllUtils 'firstidx';
use Lingua::EN::Numbers qw(num2en);

my $limit = 1 + int 5e6 / 8;
my @f = map { $_*8 - 1 } 1..$limit;

my($p0,$p1, @F) = (1,0, $f[0]);
do {
    push @F, ($f[$p0] < $F[$p1]*4) ? $f[$p0++] : $F[$p1++]*4;
} until $p0 == $limit or $p1 == $limit;

printf "First %s forbidden numbers:\n", num2en 50;
print sprintf(('%4d')x50, @F[0..49]) =~ s/.{40}\K(?=.)/\n/gr;
print "\n\n";

for my $x (5e2, 5e3, 5e4, 5e5, 5e6) {
    printf "%6d = forbidden number count up to %s\n", (firstidx { $_ > $x } @F), num2en($x);
}
Output:
First fifty forbidden numbers:
   7  15  23  28  31  39  47  55  60  63
  71  79  87  92  95 103 111 112 119 124
 127 135 143 151 156 159 167 175 183 188
 191 199 207 215 220 223 231 239 240 247
 252 255 263 271 279 284 287 295 303 311

    82 = forbidden number count up to five hundred
   831 = forbidden number count up to five thousand
  8330 = forbidden number count up to fifty thousand
 83331 = forbidden number count up to five hundred thousand
833329 = forbidden number count up to five million

Phix

Translation of: Pascal
function forbidden(integer n)
    while n>7 and remainder(n,4)=0 do
        n /= 4
    end while
    return remainder(n,8) = 7
end function

function forbidden_le_count(atom lmt)
    atom result = 0, Power4 = 1, n = 7
    while n<=lmt do
        atom delta := Power4*8
        result += floor((lmt-n)/delta)+1;
        Power4 *= 4;
        n = Power4*7
    end while
    return result
end function

sequence f50 = {}
integer i = 0
while length(f50)<50 do
    if forbidden(i) then f50 &= i end if
    i += 1
end while
printf(1,"The first 50 forbidden numbers are:\n%s\n",
         {join_by(f50,1,10," ",fmt:="%3d")})
for t=2 to iff(machine_bits()=32?16:17) do
    atom lmt = 5*power(10,t), count = forbidden_le_count(lmt)
    printf(1,"Forbidden numbers up to %,d: %,d\n",{lmt,count})
end for
Output:
The first 50 forbidden numbers are:
  7  15  23  28  31  39  47  55  60  63
 71  79  87  92  95 103 111 112 119 124
127 135 143 151 156 159 167 175 183 188
191 199 207 215 220 223 231 239 240 247
252 255 263 271 279 284 287 295 303 311

Forbidden numbers up to 500: 82
Forbidden numbers up to 5,000: 831
Forbidden numbers up to 50,000: 8,330
Forbidden numbers up to 500,000: 83,331
Forbidden numbers up to 5,000,000: 833,329
Forbidden numbers up to 50,000,000: 8,333,330
Forbidden numbers up to 500,000,000: 83,333,328
Forbidden numbers up to 5,000,000,000: 833,333,330
Forbidden numbers up to 50,000,000,000: 8,333,333,327
Forbidden numbers up to 500,000,000,000: 83,333,333,328
Forbidden numbers up to 5,000,000,000,000: 833,333,333,327
Forbidden numbers up to 50,000,000,000,000: 8,333,333,333,327
Forbidden numbers up to 500,000,000,000,000: 83,333,333,333,326
Forbidden numbers up to 5,000,000,000,000,000: 833,333,333,333,327
Forbidden numbers up to 50,000,000,000,000,000: 8,333,333,333,333,325
Forbidden numbers up to 500,000,000,000,000,000: 83,333,333,333,333,323

Python

From Michael S. Branicky's example at oeis.org/A004215

""" rosettacode.org/wiki/Forbidden_numbers """

def isforbidden(num):
    """ true if num is a forbidden number """
    fours, pow4 = num, 0
    while fours > 1 and fours % 4 == 0:
        fours //= 4
        pow4 += 1
    return (num // 4**pow4) % 8 == 7


f500k = list(filter(isforbidden, range(500_001)))

for idx, fbd in enumerate(f500k[:50]):
    print(f'{fbd: 4}', end='\n' if (idx + 1) % 10 == 0 else '')

for fbmax in [500, 5000, 50_000, 500_000]:
    print(
        f'\nThere are {sum(x <= fbmax for x in f500k):,} forbidden numbers <= {fbmax:,}.')
Output:
   7  15  23  28  31  39  47  55  60  63
  71  79  87  92  95 103 111 112 119 124
 127 135 143 151 156 159 167 175 183 188
 191 199 207 215 220 223 231 239 240 247
 252 255 263 271 279 284 287 295 303 311

There are 82 forbidden numbers <= 500.

There are 831 forbidden numbers <= 5,000.

There are 8,330 forbidden numbers <= 50,000.

There are 83,331 forbidden numbers <= 500,000.

Raku

use Lingua::EN::Numbers;
use List::Divvy;

my @f = (1..*).map:8-1;

my @forbidden = lazy flat @f[0], gather for ^∞ {
    state ($p0, $p1) = 1, 0;
    take (@f[$p0] < @forbidden[$p14) ?? @f[$p0++] !! @forbidden[$p1++]×4;
}

put "First fifty forbidden numbers: \n" ~
  @forbidden[^50].batch(10)».fmt("%3d").join: "\n";

put "\nForbidden number count up to {.Int.&cardinal}: " ~
  comma +@forbidden.&upto: $_ for 5e2, 5e3, 5e4, 5e5, 5e6;
Output:
First fifty forbidden numbers: 
  7  15  23  28  31  39  47  55  60  63
 71  79  87  92  95 103 111 112 119 124
127 135 143 151 156 159 167 175 183 188
191 199 207 215 220 223 231 239 240 247
252 255 263 271 279 284 287 295 303 311

Forbidden number count up to five hundred: 82

Forbidden number count up to five thousand: 831

Forbidden number count up to fifty thousand: 8,330

Forbidden number count up to five hundred thousand: 83,331

Forbidden number count up to five million: 833,329

RPL

Translation of: XPLo
≪ R→B
   WHILE #0d OVER ≠ LAST #3d AND == AND REPEAT SR SR END
   #7d AND B→R 7 ==
≫ 'FORB?' STO

≪ 0 1 3 PICK FOR j IF j FORB? THEN 1 + END NEXT 
≫ 'NFORB' STO
≪ { } 0 WHILE OVER SIZE 50 < REPEAT 1 + IF DUP FORB? THEN SWAP OVER + SWAP END END DROP ≫
500 NFORB
5000 NFORB
50000 NFORB
Output:
4: { 7 15 23 28 31 39 47 55 60 63 71 79 87 92 95 103 111 112 119 124 127 135 143 151 156 159 167 175 183 188 191 199 207 215 220 223 231 239 240 247 252 255 263 271 279 284 287 295 303 311 }
3: 82
2: 831
1: 8330

Sidef

say ("First 50 terms: ", 50.by { .squares_r(3) == 0 })

for n in (500, 5_000, 50_000, 500_000) {
    var v = (1..n -> count {|n|
        idiv(n, ipow(4, n.valuation(4))).is_congruent(7, 8)
    })
    say "There are #{v} forbidden numbers up to #{n.commify}."
}
Output:
First 50 terms: [7, 15, 23, 28, 31, 39, 47, 55, 60, 63, 71, 79, 87, 92, 95, 103, 111, 112, 119, 124, 127, 135, 143, 151, 156, 159, 167, 175, 183, 188, 191, 199, 207, 215, 220, 223, 231, 239, 240, 247, 252, 255, 263, 271, 279, 284, 287, 295, 303, 311]
There are 82 forbidden numbers up to 500.
There are 831 forbidden numbers up to 5,000.
There are 8330 forbidden numbers up to 50,000.
There are 83331 forbidden numbers up to 500,000.

Wren

Version 1

Library: Wren-fmt

This uses a sieve to filter out those numbers which are the sums of one, two or three squares. Works but very slow (c. 52 seconds).

import "./fmt" for Fmt

var forbidden = Fn.new { |limit, countOnly|
    var sieve = List.filled(limit+1, false)
    var ones
    var twos
    var threes
    var i = 0
    while((ones = i*i) <= limit) {
        sieve[ones] = true
        var j = i
        while ((twos = ones + j*j) <= limit) {
            sieve[twos] = true
            var k = j
            while ((threes = twos + k*k) <= limit) {
                sieve[threes] = true
                k = k + 1
            }
            j = j + 1
        }
        i = i + 1
    }
    if (countOnly) return sieve.count { |b| !b }
    var forbidden = []
    for (i in 1..limit) {
        if (!sieve[i]) forbidden.add(i)
    }
    return forbidden
}

System.print("The first 50 forbidden numbers are:")
Fmt.tprint("$3d", forbidden.call(400, false).take(50), 10)
System.print()
for (limit in [500, 5000, 50000, 500000, 5000000]) {
     var count = forbidden.call(limit, true)
     Fmt.print("Forbidden number count <= $,9d: $,7d", limit, count)
}
Output:
The first 50 forbidden numbers are:
  7  15  23  28  31  39  47  55  60  63
 71  79  87  92  95 103 111 112 119 124
127 135 143 151 156 159 167 175 183 188
191 199 207 215 220 223 231 239 240 247
252 255 263 271 279 284 287 295 303 311

Forbidden number count <=       500:      82
Forbidden number count <=     5,000:     831
Forbidden number count <=    50,000:   8,330
Forbidden number count <=   500,000:  83,331
Forbidden number count <= 5,000,000: 833,329

Version 2

This is a translation of the formula-based Python code in the OEIS link which at around 1.1 seconds is almost 50 times faster than Version 1 and is also about 3 times faster than the PARI code in that link.

import "./fmt" for Fmt

var isForbidden = Fn.new { |n|
    var m = n
    var v = 0
    while (m > 1 && m % 4 == 0) {
        m = (m/4).floor
        v = v + 1
    }
    return (n/4.pow(v)).floor % 8 == 7
} 

var f400 = (1..400).where { |i| isForbidden.call(i) }
System.print("The first 50 forbidden numbers are:")
Fmt.tprint("$3d", f400.take(50), 10)
System.print()
for (limit in [500, 5000, 50000, 500000, 5000000]) {
     var count = (1..limit).count { |i| isForbidden.call(i) }
     Fmt.print("Forbidden number count <= $,9d: $,7d", limit, count)
}
Output:
Same as Version 1

XPL0

Runtime is around 7.5 seconds on a Raspberry Pi4.

include xpllib; \for RlOutC

func Forbidden(N);      \Return 'true' if N is a forbidden number
int  N;
[while (N&3) = 0 and N do N:= N>>2;
return (N&7) = 7;
];

int N, Count, Limit;
[Text(0, "The first 50 forbidden numbers are:^m^j");
Format(4, 0);
N:= 0;  Count:= 0;
while Count < 50 do
    [if Forbidden(N) then
        [RlOut(0, (float(N)));
        Count:= Count+1;
        if rem(Count/10) = 0 then CrLf(0);
        ];
    N:= N+1;
    ];
CrLf(0);
Format(9, 0);
N:= 1;  Count:= 0;  Limit:= 500;
loop    [if Forbidden(N) then Count:= Count+1;
        if N = Limit then
            [Text(0, "Forbidden number count <= ");
            RlOutC(0, float(Limit));
            RlOutC(0, float(Count));
            CrLf(0);
            if Limit = 500_000_000 then quit;
            Limit:= Limit * 10;
            ];
        N:= N+1;
        ];
]
Output:
The first 50 forbidden numbers are:
   7  15  23  28  31  39  47  55  60  63
  71  79  87  92  95 103 111 112 119 124
 127 135 143 151 156 159 167 175 183 188
 191 199 207 215 220 223 231 239 240 247
 252 255 263 271 279 284 287 295 303 311

Forbidden number count <=         500         82
Forbidden number count <=       5,000        831
Forbidden number count <=      50,000      8,330
Forbidden number count <=     500,000     83,331
Forbidden number count <=   5,000,000    833,329
Forbidden number count <=  50,000,000  8,333,330
Forbidden number count <= 500,000,000 83,333,328