Forbidden numbers: Difference between revisions
Drkameleon (talk | contribs) added Arturo |
→{{header|Free Pascal}}: inserted iterative counting with optimized isForbidden. 500,000,000 in 1,4 s @TIO.RUN |
||
Line 346: | Line 346: | ||
uses |
uses |
||
sysutils,strutils; |
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; |
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 |
var |
||
Power4, |
Power4,delta,n : NativeUint; |
||
begin |
begin |
||
result := 0; |
result := 0; |
||
power4 := 1; |
power4 := 1; |
||
repeat |
repeat |
||
delta := Power4*8; |
delta := Power4*8;// j = 1 |
||
n := Power4*7; |
n := Power4*7; |
||
if n > lmt then |
if n > lmt then |
||
Break; |
Break; |
||
//max j to reach limit |
|||
inc(result,(lmt-n) DIV delta+1); |
inc(result,(lmt-n) DIV delta+1); |
||
Power4 *= 4; |
Power4 *= 4; |
||
until false; |
until false; |
||
end; |
end; |
||
var |
var |
||
lmt,n: NativeUint; |
lmt,n,cnt: NativeUint; |
||
BEGIN |
BEGIN |
||
writeln('First fifty forbidden numbers:'); |
writeln('First fifty forbidden numbers:'); |
||
n := 1; |
n := 1; |
||
lmt := |
lmt := 0; |
||
repeat; |
|||
if isForbidden(n) then |
|||
Begin |
|||
⚫ | |||
inc(lmt); |
|||
if LMT MOD 20 = 0 THEN |
|||
⚫ | |||
end; |
|||
n +=1; |
|||
until lmt >= 50; |
|||
⚫ | |||
writeln; |
|||
writeln('count of forbidden numbers below iterative'); |
|||
n := 1; |
|||
cnt := 0; |
|||
lmt := 5; |
|||
repeat |
repeat |
||
repeat |
repeat; |
||
//if isForbidden(n) then cnt+=1 takes to long 100% -> 65% of time |
|||
⚫ | |||
inc(cnt,ORD(isForbidden(n))); |
|||
until CntForbiddenTilLimit(n) >= lmt; |
|||
n += 1; |
|||
until n >= lmt; |
|||
writeln(Numb2USA(IntToStr(lmt)):30,Numb2USA(IntToStr(Cnt)):25); |
|||
⚫ | |||
lmt |
lmt *= 10; |
||
until lmt > |
until lmt > 500*1000*1000; |
||
writeln; |
writeln; |
||
⚫ | |||
writeln('count of forbidden numbers below '); |
writeln('count of forbidden numbers below '); |
||
lmt := 5; |
lmt := 5; |
||
repeat |
repeat |
||
writeln(Numb2USA(IntToStr(lmt)):30,Numb2USA(IntToStr(CntForbiddenTilLimit(lmt))):25); |
writeln(Numb2USA(IntToStr(lmt)):30,Numb2USA(IntToStr(CntForbiddenTilLimit(lmt))):25); |
||
lmt *= 10; |
lmt *= 10; |
||
Line 390: | Line 419: | ||
{{out|@TIO.RUN}} |
{{out|@TIO.RUN}} |
||
<pre> |
<pre> |
||
First fifty forbidden numbers: |
First fifty forbidden numbers: |
||
7 15 23 28 31 39 47 55 60 63 71 79 87 92 95 103 111 112 119 124 |
7 15 23 28 31 39 47 55 60 63 71 79 87 92 95 103 111 112 119 124 |
||
Line 395: | Line 425: | ||
252 255 263 271 279 284 287 295 303 311 |
252 255 263 271 279 284 287 295 303 311 |
||
count of forbidden numbers below |
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 |
5 0 |
||
50 7 |
50 7 |
||
Line 414: | Line 455: | ||
50,000,000,000,000,000 8,333,333,333,333,325 |
50,000,000,000,000,000 8,333,333,333,333,325 |
||
500,000,000,000,000,000 83,333,333,333,333,323 |
500,000,000,000,000,000 83,333,333,333,333,323 |
||
Real time: |
Real time: 1.392 s User time: 1.343 s Sys. time: 0.042 s CPU share: 99.52 %</pre> |
||
</pre> |
|||
=={{header|Python}}== |
=={{header|Python}}== |
Revision as of 10:53, 29 April 2023
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
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
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
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
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
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
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 %
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[$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
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