Forbidden numbers
You are encouraged to solve this task according to the task description, using any language you may know.
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 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
ABC
HOW TO REPORT is.forbidden n:
WHILE n>1 AND n mod 4 = 0:
PUT floor (n/4) IN n
REPORT n mod 8 = 7
HOW TO RETURN next.forbidden n:
PUT n+1 IN n
WHILE NOT is.forbidden n: PUT n+1 IN n
RETURN n
PUT 0 IN n
PUT 0 IN count
PUT 500 IN limit
WHILE limit < 50000000:
PUT next.forbidden n IN n
IF count<50:
WRITE n>>5
IF count mod 10=9: WRITE/
IF n >= limit:
WRITE "Forbidden numbers <= ", limit>>10, ": ", count>>10/
PUT limit*10 IN limit
PUT count+1 IN count- 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 Forbidden numbers <= 500: 82 Forbidden numbers <= 5000: 831 Forbidden numbers <= 50000: 8330 Forbidden numbers <= 500000: 83331 Forbidden numbers <= 5000000: 833329
ALGOL 68
Based on the Wren version of the example at oeis.org/A004215
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
ALGOL W
Based on the Wren version of the example at oeis.org/A004215
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 %
logical procedure isForbidden ( integer value n ) ;
begin
integer m, p4;
m := n;
p4 := 1;
while m > 1 and m rem 4 = 0 do begin
m := m div 4;
p4 := p4 * 4
end while_m_gt_1_and_m_rem_4_eq_0 ;
( n div p4 ) rem 8 = 7
end isForbidden ;
begin % show the first 50 forbidden numbers and forbidden number counts %
integer fCount, nextToShow;
fCount := 0;
nextToShow := 500;
for i := 1 until 500000000 do begin
if isForbidden( i ) then begin
fCount := fCount + 1;
if fCount <= 50 then begin
writeon( i_w := 4, s_w := 0, " ", i );
if fCount rem 10 = 0 then write()
end if_fCount_le_50
end if_isForbidden ;
if i = nextToShow then begin
write( i_w := 8, s_w := 0, "There are ", fCount, " Forbidden numbers up to "
, i_w := 1, i
);
nextToShow := nextToShow * 10
end if_i_eq_nextToShow
end for_i
end
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 There are 83333328 Forbidden numbers up to 500000000
AppleScript
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
AWK
Based on the Wren version of the example at oeis.org/A004215
BEGIN \
{
fCount = 0;
nextToShow = 500;
for( i = 1; i <= 5000000; i ++ )
{
if( isForbidden( i ) )
{
fCount += 1;
if( fCount <= 50 )
{
printf( " %3d%s", i, ( ( fCount % 10 == 0 ) ? "\n" : "" ) );
}
}
if( i == nextToShow )
{
printf( "There are %8d Forbidden numbers up to %d\n", fCount, i );
nextToShow *= 10;
}
}
}
function isForbidden( n, m, p4 )
{
m = n;
p4 = 1;
while( m > 1 && m % 4 == 0 )
{
m = int( m / 4 );
p4 *= 4;
}
return int( n / p4 ) % 8 == 7;
}
- 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
BASIC
10 DEFINT A-Z
20 N=0: C=0: L=500
30 N=N+1
40 M=N
50 IF M>1 AND (M AND 3)=0 THEN M=M\4: GOTO 50
60 IF (M AND 7)<>7 THEN 30
70 C=C+1
80 IF C<=50 THEN PRINT USING "#####";N;: IF C MOD 10=0 THEN PRINT
90 IF N<L THEN 30
100 PRINT USING "Forbidden numbers up to #####: #####";L;C
110 IF L<5000 THEN L=L*10: GOTO 30
120 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 Forbidden numbers up to 500: 83 Forbidden numbers up to 5000: 832
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
CLU
forbidden_numbers = iter () yields (int)
n: int := 1
while true do
m: int := n
while m>1 cand m//4 = 0 do
m := m/4
end
if m//8 = 7 then yield(n) end
n := n+1
end
end forbidden_numbers
start_up = proc ()
po: stream := stream$primary_output()
count: int := 0
lim: int := 500
for n: int in forbidden_numbers() do
if count < 50 then
stream$putright(po, int$unparse(n), 5)
if count//10 = 9 then stream$putl(po, "") end
end
if n >= lim then
stream$puts(po, "Forbidden numbers <= ")
stream$putright(po, int$unparse(lim), 10)
stream$puts(po, ": ")
stream$putright(po, int$unparse(count), 10)
stream$putl(po, "")
if lim > 50000000 then break end
lim := lim * 10
end
count := count + 1
end
end start_up- 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 Forbidden numbers <= 500: 82 Forbidden numbers <= 5000: 831 Forbidden numbers <= 50000: 8330 Forbidden numbers <= 500000: 83331 Forbidden numbers <= 5000000: 833329 Forbidden numbers <= 50000000: 8333330 Forbidden numbers <= 500000000: 83333328
Cowgol
include "cowgol.coh";
sub is_forbidden(n: uint32): (r: uint8) is
r := 0;
while n>1 and n&3 == 0 loop
n := n>>2;
end loop;
if n&7 == 7 then
r := 1;
end if;
end sub;
sub next_forbidden(n: uint32): (r: uint32) is
loop
n := n + 1;
if is_forbidden(n) != 0 then
break;
end if;
end loop;
r := n;
end sub;
var n: uint32 := 0;
var count: uint32 := 0;
var lim: uint32 := 500;
loop
n := next_forbidden(n);
if count < 50 then
print_i32(n);
if count % 10 != 9 then
print_char('\t');
else
print_nl();
end if;
end if;
if n >= lim then
print("Forbidden numbers <= ");
print_i32(lim);
print(": ");
print_i32(count);
print_nl();
if lim == 50000000 then
break;
end if;
lim := lim * 10;
end if;
count := count + 1;
end loop;- 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 Forbidden numbers <= 500: 82 Forbidden numbers <= 5000: 831 Forbidden numbers <= 50000: 8330 Forbidden numbers <= 500000: 83331 Forbidden numbers <= 5000000: 833329 Forbidden numbers <= 50000000: 8333330
Draco
proc forbidden(word n) bool:
while n>1 and n&3 = 0 do
n := n >> 2
od;
n & 7 = 7
corp
proc main() void:
word n, count, lim;
n := 0;
count := 0;
lim := 500;
while
while n := n+1; not forbidden(n) do od;
count := count + 1;
if count <= 50 then
write(n:5);
if count % 10 = 0 then writeln() fi
fi;
if n >= lim then
writeln("Forbidden numbers <= ", lim:5, ": ", count-1:5);
if lim < 50000 then lim := lim * 10 fi
fi;
n < lim
do od
corp- 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 Forbidden numbers <= 500: 82 Forbidden numbers <= 5000: 831 Forbidden numbers <= 50000: 8330
EasyLang
Based on the Wren version of the example at oeis.org/A004215
func isForbidden n .
m = n
p4 = 1
while m > 1 and m mod 4 = 0
m /= 4
p4 *= 4
.
if (n / p4) mod 8 = 7
return 1
.
return 0
.
fCount = 0
nextToShow = 500
for i = 1 to 500000000
if isForbidden i = 1
fCount += 1
if fCount <= 50
write " " & i
if fCount mod 10 = 0 : print ""
.
.
if i = nextToShow
print "There are " & fCount & " Forbidden numbers up to " & i
nextToShow *= 10
.
.- 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
Forth
: forbidden? ( n -- f )
begin
dup 1 > if dup 3 and 0= else false then
while
2 rshift
repeat
7 and 7 =
;
: main
." The first 50 forbidden numbers are:" cr
0 0
begin
over 50 <
while
dup forbidden? if
dup 3 .r
swap 1+ swap
over 10 mod 0= if cr else space then
then
1+
repeat
cr
500 >r
begin
r@ 5000000 <=
while
dup forbidden? if
dup r@ > if
." There are " over . ." forbidden numbers <= " r@ . cr
r> 10 * >r
then
swap 1+ swap
then
1+
repeat
rdrop 2drop
;
main
bye
- 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 numbers <= 500 There are 831 forbidden numbers <= 5000 There are 8330 forbidden numbers <= 50000 There are 83331 forbidden numbers <= 500000 There are 833329 forbidden numbers <= 5000000
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.
FutureBasic
local fn isForbidden (num As uint64) As bool
uint64 fours = num, pow4 = 0
bool result = _false
While (fours > 1) And (fours Mod 4 = 0)
fours = fours \ 4
pow4 += 1
Wend
if (num \ 4 ^ pow4) Mod 8 = 7 then result = _true
End fn = result
uint64 cnt = 0
uint64 i = 0
Print "The first 50 forbidden numbers are:"
Do
i += 1
If fn isForbidden(i)
cnt += 1
If cnt <= 50 Then Print Using "####"; i; : If cnt Mod 10 = 0 Then Print
End If
If i = 500 Then Print "Forbidden number count " + str$(i) + str$(cnt)
If i = 5e3 Or i = 5e4 Or i = 5e5 Or i = 5e6 Then Print "Forbidden number count " + str$(i) + str$(cnt)
Until i => 5e6
handleevents- 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
Go
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
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
""" 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.
Miranda
main :: [sys_message]
main = [Stdout (table 10 5 (take 50 forbidden)),
Stdout "\n",
Stdout (lay (map showcount [500, 5000, 50000, 500000]))
]
where showcount n = show (forbidden_count n) ++
" forbidden numbers <= " ++
show n
table :: num->num->[num]->[char]
table w cw ls = lay (map (concat . map (rjustify cw . show)) (split w ls))
split :: num->[*]->[[*]]
split n [] = []
split n ls = take n ls:split n (drop n ls)
forbidden_count :: num->num
forbidden_count n = #takewhile (<=n) forbidden
forbidden :: [num]
forbidden = filter is_forbidden [1..]
is_forbidden :: num->bool
is_forbidden n = False, if d=0
= is_forbidden d, if m=0
= n mod 8 = 7, otherwise
where d = n div 4
m = n mod 4- 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 82 forbidden numbers <= 500 831 forbidden numbers <= 5000 8330 forbidden numbers <= 50000 83331 forbidden numbers <= 500000
Lua
Based on the Wren version of the example at oeis.org/A004215
Uses the // (floor division) operator of Lua versions 5.3.* and 5.4.* - replace a // b with math.floor( a / b ) for earlier versions.
do
local isForbidden = function( n )
local m, p4 = n, 1
while m > 1 and m % 4 == 0 do
m, p4 = ( m // 4 ), p4 * 4
end
return ( n // p4 ) % 8 == 7
end
local fCount, nextToShow = 0, 500
for i = 1, 50000000 do
if isForbidden( i ) then
fCount = fCount + 1
if fCount <= 50 then
io.write( string.format( " %3d", i ) )
if fCount % 10 == 0 then io.write( "\n" ) end
end
end
if i == nextToShow then
io.write( "There are "..string.format( "%8d", fCount ).." Forbidden numbers up to "..i.."\n" )
nextToShow = nextToShow * 10
end
end
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
Modula-2
MODULE ForbiddenNumbers;
FROM InOut IMPORT WriteCard, WriteString, WriteLn;
VAR n, count, limit: CARDINAL;
PROCEDURE forbidden(n: CARDINAL): BOOLEAN;
BEGIN
WHILE (n>0) AND (n MOD 4 = 0) DO
n := n DIV 4
END;
RETURN n MOD 8 = 7
END forbidden;
BEGIN
n := 0;
count := 0;
limit := 500;
LOOP
REPEAT INC(n) UNTIL forbidden(n);
IF count < 50 THEN
WriteCard(n, 5);
IF count MOD 10=9 THEN WriteLn END
END;
IF n >= limit THEN
WriteString("Forbidden numbers <= ");
WriteCard(limit, 5);
WriteString(": ");
WriteCard(count, 5);
WriteLn;
IF limit = 50000 THEN EXIT END;
limit := limit * 10
END;
INC(count)
END
END ForbiddenNumbers.
- 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 Forbidden numbers <= 500: 82 Forbidden numbers <= 5000: 831 Forbidden numbers <= 50000: 8330
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
##
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
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
PL/M
100H:
BDOS: PROCEDURE (FN,ARG); DECLARE FN BYTE, ARG ADDRESS; GO TO 5; END BDOS;
EXIT: PROCEDURE; GO TO 0; END EXIT;
PRINT: PROCEDURE (STR); DECLARE STR ADDRESS; CALL BDOS(9,STR); END PRINT;
PRINT$NUM: PROCEDURE (N, WIDTH);
DECLARE N ADDRESS, (WIDTH, I) BYTE;
DECLARE S (6) BYTE INITIAL ('.....$');
I = 5;
DIGIT:
S(I := I-1) = N MOD 10 + '0';
IF (N := N/10) > 0 THEN GO TO DIGIT;
DO WHILE I>0;
S(I := I-1) = ' ';
END;
CALL PRINT(.S(5-WIDTH));
END PRINT$NUM;
FORBIDDEN: PROCEDURE (N) BYTE;
DECLARE N ADDRESS;
DO WHILE N>1 AND (N AND 3)=0;
N = SHR(N, 2);
END;
RETURN (N AND 7) = 7;
END FORBIDDEN;
DECLARE (N, COUNT, LIM) ADDRESS INITIAL (0, 0, 500);
DO WHILE 1;
DO WHILE NOT FORBIDDEN(N := N+1); END;
COUNT = COUNT + 1;
IF COUNT<=50 THEN DO;
CALL PRINT$NUM(N, 4);
IF COUNT MOD 10=0 THEN CALL PRINT(.(13,10,'$'));
END;
IF N >= LIM THEN DO;
CALL PRINT(.'FORBIDDEN NUMBER COUNT <= $');
CALL PRINT$NUM(LIM, 5);
CALL PRINT(.': $');
CALL PRINT$NUM(COUNT-1, 5);
CALL PRINT(.(13,10,'$'));
IF LIM = 50000 THEN CALL EXIT;
LIM = LIM * 10;
END;
END;
EOF- 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 FORBIDDEN NUMBER COUNT <= 500: 82 FORBIDDEN NUMBER COUNT <= 5000: 831 FORBIDDEN NUMBER COUNT <= 50000: 8330
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.
Quackery
[ 0 over
[ dup 1 > while
dup 4 mod 0 = while
4 / dip 1+ again ]
drop
4 swap ** / 8 mod 7 = ] is forbidden ( n --> b )
say "First 50 forbidden numbers:"
0 []
[ over forbidden if
[ over number$ nested join ]
dip 1+
dup size 50 = until ]
nip 50 wrap$
cr cr
' [ 500 5000 50000 500000 5000000 ]
witheach
[ dup say "Forbidden numbers <= "
echo say ": "
0 swap times
[ i^ 1+ forbidden + ]
echo cr ]- 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 Forbidden numbers <= 500: 82 Forbidden numbers <= 5000: 831 Forbidden numbers <= 50000: 8330 Forbidden numbers <= 500000: 83331 Forbidden numbers <= 5000000: 833329
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
Refal
$ENTRY Go {
= <First50>
<Count 500 5000 50000 500000>;
};
First50 {
= <Each (Prout) <Grp 10 <Gen 50 Forbidden>>>;
};
Count {
(s.C s.N) = ;
(s.C s.N) s.L e.Ls,
<Add 1 s.C> <Next Forbidden s.N>: e.Next,
<Compare s.N s.L>: {
'-' = <Count (e.Next) s.L e.Ls>;
s.1 = <Prout "Forbidden numbers up to "
<Fmt 6 s.L> ": " <Fmt 6 s.C>>
<Count (e.Next) e.Ls>;
};
e.X = <Count (0 <Next Forbidden 1>) e.X>;
};
Forbidden {
s.N, <Divmod s.N 4>: {
(0) s.D = False;
(s.M) 0 = <Forbidden s.M>;
(s.M) s.D, <Mod s.N 8>: 7 = True;
(s.M) s.D = False;
};
};
Gen {
s.I s.F = <Gen s.I s.F 1>;
0 s.F s.N = ;
s.I s.F s.N, <Mu s.F s.N>: True =
s.N <Gen <Sub s.I 1> s.F <Add s.N 1>>;
s.I s.F s.N = <Gen s.I s.F <Add s.N 1>>;
};
Next {
s.F s.N, <Add 1 s.N>: s.G, <Mu s.F s.G>: {
True = s.G;
False = <Next s.F s.G>;
};
};
Rpt { 0 s.C = ; s.N s.C = s.C <Rpt <Sub s.N 1> s.C>; };
Fmt { s.W s.N, <Last s.W <Rpt s.W ' '> <Symb s.N>>: (e.1) e.2 = e.2; };
Grp { s.S = ; s.S e.X, <First s.S e.X>: (e.1) e.2 = (e.1) <Grp s.S e.2>; };
Each { (e.F) = ; (e.F) t.X e.Xs = <Mu e.F t.X> <Each (e.F) e.Xs>; };- 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 ) Forbidden numbers up to 500: 82 Forbidden numbers up to 5000: 831 Forbidden numbers up to 50000: 8330 Forbidden numbers up to 500000: 83331
RPL
≪ 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
SETL
program forbidden_numbers;
n := 0;
count := 0;
lim := 500;
loop until lim > 50000000 do
loop until forbidden(n) do n +:= 1; end loop;
count +:= 1;
if count <= 50 then
nprint(lpad(str n, 5));
if count mod 10=0 then print; end if;
end if;
if n >= lim then
print("Forbidden numbers up to " + lpad(str lim, 10) +
": " + lpad(str (count-1), 10));
lim *:= 10;
end if;
end loop;
proc forbidden(n);
loop while n>1 and n mod 4 = 0 do
n div:= 4;
end loop;
return n mod 8 = 7;
end proc;
end program;- 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 Forbidden numbers up to 500: 82 Forbidden numbers up to 5000: 831 Forbidden numbers up to 50000: 8330 Forbidden numbers up to 500000: 83331 Forbidden numbers up to 5000000: 833329 Forbidden numbers up to 50000000: 8333330
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
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