Forbidden numbers

Forbidden numbers 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.

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

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 other than summing 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 less 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

• Wikipedia: Lagrange's theorem
• OEIS A004215 - Numbers that are the sum of 4 but no fewer nonzero squares

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

Python

""" 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[$p1]×4) ?? @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

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