Supporting files for the How Fast Is Crow? article · the programs in four languages, the harness and the raw timings
README.md
# Crow benchmark sources
The programs, harness and raw results behind the performance tables in
"Crow, a Programming Language" and "Crow vs. Other Languages" on aggnr.com.
Measured 2026-10-21 on an Apple M4 Max (16 cores, 64 GB), macOS 26.5,
with crow 0.7.0 and the changes after it (commit d3ee350), rustc 1.93.1,
go 1.23.0 and CPython 3.12.12. Every binary had run once before the timed
runs, since macOS checks a new executable at its first launch.
## Files
- `*.co` — the Crow programs. Sizes are written into the files
(`fannkuch[10]`, `n = 50000000` in `nbody.co`, `n = 1000000` in
`nbody_1m.co`, `n = 5500`, `n = 5000`, `n = 18`).
- `*.rs`, `*.go`, `*.py` — hand-written ports of the same programs. The
Rust and Go ports of spectral norm, k-nucleotide, reverse complement and
binary trees take a thread count as their last argument (16 here); the
Python ports of spectral norm and binary trees take a process count.
- `fasta.go` — makes the input for k-nucleotide and reverse complement:
`./fasta_go 5000000 1 > fasta.txt` (50.8 MB).
- `pidigits-Cargo.toml` — the Rust pidigits port needs the `num-bigint`
crate; build it as a cargo project with this manifest.
- `binarytrees-mimalloc-Cargo.toml` — the Rust binary-trees port linked
against mimalloc (prepend
`use mimalloc::MiMalloc; #[global_allocator] static GLOBAL: MiMalloc = MiMalloc;`
to the source).
- `bench5.py` — the harness: five runs per program and language (three for
the interpreter), median, min and max, and a check that every port's
output is byte-identical.
- `results5.json` — the raw timings this run produced.
## Building
```
crow build fannkuch.co # and the others; crow build = opt-level 3 + LTO
rustc -C opt-level=3 -C target-cpu=native -o fannkuch_rs fannkuch.rs
go build -o fannkuch_go fannkuch.go
python bench5.py # expects the binaries named as above
```
bench5.py
"""Five runs per program and language; median, min and max. Interpreted
Crow three runs. Writes results5.json and prints a table."""
import hashlib
import json
import os
import statistics
import subprocess
import time
PY = os.path.expanduser("~/.venvs/finch_ml/bin/python")
CROW = os.path.expanduser("~/.cargo/bin/crow")
def once(cmd, stdin=None):
f = open(stdin, "rb") if stdin else None
t0 = time.perf_counter()
r = subprocess.run(cmd, stdin=f, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
dt = time.perf_counter() - t0
if f:
f.close()
if r.returncode != 0:
raise RuntimeError(f"{cmd}: exit {r.returncode}: {r.stderr.decode(errors='replace')[:200]}")
return dt, hashlib.md5(r.stdout).hexdigest()
def runs(cmd, n=5, stdin=None):
times, digests = [], set()
for _ in range(n):
dt, d = once(cmd, stdin)
times.append(dt)
digests.add(d)
assert len(digests) == 1, f"{cmd}: output changed between runs"
return times, digests.pop()
CASES = [
("fannkuch", "n = 10", {
"crow": ["./fannkuch"], "rust": ["./fannkuch_rs", "10"], "go": ["./fannkuch_go", "10"],
"cpython": [PY, "fannkuch.py", "10"], "crow -i": [CROW, "run", "-i", "fannkuch.co"]}, None),
("n-body", "50 million steps", {
"crow": ["./nbody"], "rust": ["./nbody_rs", "50000000"], "go": ["./nbody_go", "50000000"]}, None),
("n-body", "1 million steps", {
"crow": ["./nbody_1m"], "rust": ["./nbody_rs", "1000000"], "go": ["./nbody_go", "1000000"],
"cpython": [PY, "nbody.py", "1000000"], "crow -i": [CROW, "run", "-i", "nbody_1m.co"]}, None),
("spectral norm", "n = 5,500", {
"crow": ["./spectralnorm"], "rust": ["./spectralnorm_rs", "5500", "16"], "go": ["./spectralnorm_go", "5500", "16"],
"cpython": [PY, "spectralnorm.py", "5500", "16"], "crow -i": [CROW, "run", "-i", "spectralnorm.co"]}, None),
("k-nucleotide", "5 million", {
"crow": ["./knucleotide"], "rust": ["./knucleotide_rs", "16"], "go": ["./knucleotide_go", "16"],
"cpython": [PY, "knucleotide.py"], "crow -i": [CROW, "run", "-i", "knucleotide.co"]}, "fasta.txt"),
("reverse complement", "5 million", {
"crow": ["./revcomp"], "rust": ["./revcomp_rs", "16"], "go": ["./revcomp_go", "16"],
"cpython": [PY, "revcomp.py"], "crow -i": [CROW, "run", "-i", "revcomp.co"]}, "fasta.txt"),
("pidigits", "5,000 digits", {
"crow": ["./pidigits"], "rust": ["./pidigits_crate/target/release/pidigits", "5000"], "go": ["./pidigits_go", "5000"],
"cpython": [PY, "pidigits.py", "5000"], "crow -i": [CROW, "run", "-i", "pidigits.co"]}, None),
("binary trees", "depth 18", {
"crow": ["./binarytrees"], "rust": ["./binarytrees_rs", "18", "16"], "rust (mimalloc)": ["./bt_mi/target/release/bt_mi", "18", "16"],
"go": ["./binarytrees_go", "18", "16"], "cpython": [PY, "binarytrees.py", "18", "16"],
"crow -i": [CROW, "run", "-i", "binarytrees.co"]}, None),
]
out = []
for prog, size, langs, stdin in CASES:
digests = {}
for lang, cmd in langs.items():
n = 3 if lang == "crow -i" else 5
try:
ts, d = runs(cmd, n, stdin)
except Exception as e: # noqa: BLE001
print(f"{prog:19s} {size:17s} {lang:16s} FAILED {e}", flush=True)
continue
digests[lang] = d
med = statistics.median(ts)
out.append({"program": prog, "size": size, "lang": lang, "runs": n, "median": med, "min": min(ts), "max": max(ts), "times": ts, "md5": d})
print(f"{prog:19s} {size:17s} {lang:16s} median {med:7.3f} min {min(ts):7.3f} max {max(ts):7.3f} spread {100*(max(ts)-min(ts))/med:4.0f}%", flush=True)
if len(set(digests.values())) != 1:
print(f" !! outputs differ for {prog}: {digests}", flush=True)
else:
print(f" outputs identical across {len(digests)} ports", flush=True)
json.dump(out, open("results5.json", "w"), indent=1)
print("done")
results5.json
[
{
"program": "fannkuch",
"size": "n = 10",
"lang": "crow",
"runs": 5,
"median": 0.12661862466484308,
"min": 0.12469450011849403,
"max": 0.12821362540125847,
"times": [
0.12719512544572353,
0.12643879186362028,
0.12469450011849403,
0.12821362540125847,
0.12661862466484308
],
"md5": "323202fa3c20601a3e135f4e04d8e1eb"
},
{
"program": "fannkuch",
"size": "n = 10",
"lang": "rust",
"runs": 5,
"median": 0.13213629182428122,
"min": 0.13150033447891474,
"max": 0.1324076671153307,
"times": [
0.13225845899432898,
0.13179741706699133,
0.13150033447891474,
0.1324076671153307,
0.13213629182428122
],
"md5": "323202fa3c20601a3e135f4e04d8e1eb"
},
{
"program": "fannkuch",
"size": "n = 10",
"lang": "go",
"runs": 5,
"median": 0.12332841660827398,
"min": 0.12322391662746668,
"max": 0.12442158348858356,
"times": [
0.12442158348858356,
0.12332841660827398,
0.12322391662746668,
0.1232340419664979,
0.1238984577357769
],
"md5": "323202fa3c20601a3e135f4e04d8e1eb"
},
{
"program": "fannkuch",
"size": "n = 10",
"lang": "cpython",
"runs": 5,
"median": 2.9719013329595327,
"min": 2.970015707425773,
"max": 2.981610666960478,
"times": [
2.9711843747645617,
2.970015707425773,
2.981610666960478,
2.9767775824293494,
2.9719013329595327
],
"md5": "323202fa3c20601a3e135f4e04d8e1eb"
},
{
"program": "fannkuch",
"size": "n = 10",
"lang": "crow -i",
"runs": 3,
"median": 1.306546333245933,
"min": 1.30296004191041,
"max": 1.3126017497852445,
"times": [
1.3126017497852445,
1.30296004191041,
1.306546333245933
],
"md5": "323202fa3c20601a3e135f4e04d8e1eb"
},
{
"program": "n-body",
"size": "50 million steps",
"lang": "crow",
"runs": 5,
"median": 4.2005200842395425,
"min": 4.195659499615431,
"max": 4.203349375165999,
"times": [
4.203349375165999,
4.202794540673494,
4.2005200842395425,
4.195659499615431,
4.196747042238712
],
"md5": "32c32132315472b32572eb1f52fc01a6"
},
{
"program": "n-body",
"size": "50 million steps",
"lang": "rust",
"runs": 5,
"median": 1.5316822919994593,
"min": 1.5277127092704177,
"max": 1.5345129994675517,
"times": [
1.5318213747814298,
1.5316822919994593,
1.5315814996138215,
1.5277127092704177,
1.5345129994675517
],
"md5": "32c32132315472b32572eb1f52fc01a6"
},
{
"program": "n-body",
"size": "50 million steps",
"lang": "go",
"runs": 5,
"median": 1.4819951253011823,
"min": 1.4794401666149497,
"max": 1.4821725832298398,
"times": [
1.4817267088219523,
1.4820147501304746,
1.4821725832298398,
1.4819951253011823,
1.4794401666149497
],
"md5": "32c32132315472b32572eb1f52fc01a6"
},
{
"program": "n-body",
"size": "1 million steps",
"lang": "crow",
"runs": 5,
"median": 0.08638154156506062,
"min": 0.08562891650944948,
"max": 0.08656199928373098,
"times": [
0.08656199928373098,
0.08562891650944948,
0.08638154156506062,
0.08636962529271841,
0.0865072924643755
],
"md5": "6f4826a164a3e707ddfedd4b5b6d38e2"
},
{
"program": "n-body",
"size": "1 million steps",
"lang": "rust",
"runs": 5,
"median": 0.03203687537461519,
"min": 0.03169116657227278,
"max": 0.03212766721844673,
"times": [
0.03169116657227278,
0.032065083272755146,
0.03203687537461519,
0.03212766721844673,
0.03185912501066923
],
"md5": "6f4826a164a3e707ddfedd4b5b6d38e2"
},
{
"program": "n-body",
"size": "1 million steps",
"lang": "go",
"runs": 5,
"median": 0.031459749676287174,
"min": 0.031059875153005123,
"max": 0.03165741730481386,
"times": [
0.031059875153005123,
0.031459749676287174,
0.03154745977371931,
0.03165741730481386,
0.031217041425406933
],
"md5": "6f4826a164a3e707ddfedd4b5b6d38e2"
},
{
"program": "n-body",
"size": "1 million steps",
"lang": "cpython",
"runs": 5,
"median": 2.3509047077968717,
"min": 2.3413364170119166,
"max": 2.3817698750644922,
"times": [
2.377620958723128,
2.3413364170119166,
2.3509047077968717,
2.3421438755467534,
2.3817698750644922
],
"md5": "6f4826a164a3e707ddfedd4b5b6d38e2"
},
{
"program": "n-body",
"size": "1 million steps",
"lang": "crow -i",
"runs": 3,
"median": 2.128798959776759,
"min": 2.087382416240871,
"max": 2.1333018336445093,
"times": [
2.128798959776759,
2.1333018336445093,
2.087382416240871
],
"md5": "6f4826a164a3e707ddfedd4b5b6d38e2"
},
{
"program": "spectral norm",
"size": "n = 5,500",
"lang": "crow",
"runs": 5,
"median": 0.08134200051426888,
"min": 0.0795697495341301,
"max": 0.0821498753502965,
"times": [
0.0821498753502965,
0.08134200051426888,
0.08207733277231455,
0.07960616704076529,
0.0795697495341301
],
"md5": "1584fbeab0a952f314fbf0fd7621885f"
},
{
"program": "spectral norm",
"size": "n = 5,500",
"lang": "rust",
"runs": 5,
"median": 0.07576483301818371,
"min": 0.07173583284020424,
"max": 0.07733591739088297,
"times": [
0.07173583284020424,
0.07332850061357021,
0.07579137478023767,
0.07576483301818371,
0.07733591739088297
],
"md5": "1584fbeab0a952f314fbf0fd7621885f"
},
{
"program": "spectral norm",
"size": "n = 5,500",
"lang": "go",
"runs": 5,
"median": 0.10663099959492683,
"min": 0.10592270828783512,
"max": 0.11032320838421583,
"times": [
0.10592270828783512,
0.1062610000371933,
0.10725591611117125,
0.10663099959492683,
0.11032320838421583
],
"md5": "1584fbeab0a952f314fbf0fd7621885f"
},
{
"program": "spectral norm",
"size": "n = 5,500",
"lang": "cpython",
"runs": 5,
"median": 8.844692875631154,
"min": 8.673282207921147,
"max": 8.948407417163253,
"times": [
8.948407417163253,
8.686181458644569,
8.844692875631154,
8.897352125495672,
8.673282207921147
],
"md5": "1584fbeab0a952f314fbf0fd7621885f"
},
{
"program": "spectral norm",
"size": "n = 5,500",
"lang": "crow -i",
"runs": 3,
"median": 6.570208791643381,
"min": 6.543285749852657,
"max": 6.594627249985933,
"times": [
6.594627249985933,
6.570208791643381,
6.543285749852657
],
"md5": "1584fbeab0a952f314fbf0fd7621885f"
},
{
"program": "k-nucleotide",
"size": "5 million",
"lang": "crow",
"runs": 5,
"median": 0.05431133322417736,
"min": 0.053443250246346,
"max": 0.05559858400374651,
"times": [
0.05559858400374651,
0.05431133322417736,
0.05362291540950537,
0.055286542512476444,
0.053443250246346
],
"md5": "0fe2fd08243d139606d5b50e1e16ed71"
},
{
"program": "k-nucleotide",
"size": "5 million",
"lang": "rust",
"runs": 5,
"median": 0.02782437577843666,
"min": 0.027199416421353817,
"max": 0.0286056250333786,
"times": [
0.027375875040888786,
0.027199416421353817,
0.02782437577843666,
0.0286056250333786,
0.028043375350534916
],
"md5": "0fe2fd08243d139606d5b50e1e16ed71"
},
{
"program": "k-nucleotide",
"size": "5 million",
"lang": "go",
"runs": 5,
"median": 0.1318800002336502,
"min": 0.13100854121148586,
"max": 0.1342447092756629,
"times": [
0.1342447092756629,
0.13100854121148586,
0.1318800002336502,
0.13119020778685808,
0.13189304154366255
],
"md5": "0fe2fd08243d139606d5b50e1e16ed71"
},
{
"program": "k-nucleotide",
"size": "5 million",
"lang": "cpython",
"runs": 5,
"median": 4.266343416646123,
"min": 4.169315208680928,
"max": 4.423589707352221,
"times": [
4.423589707352221,
4.169315208680928,
4.397985000163317,
4.266343416646123,
4.191954541951418
],
"md5": "0fe2fd08243d139606d5b50e1e16ed71"
},
{
"program": "k-nucleotide",
"size": "5 million",
"lang": "crow -i",
"runs": 3,
"median": 0.2608485845848918,
"min": 0.2605542913079262,
"max": 0.26258308440446854,
"times": [
0.2605542913079262,
0.2608485845848918,
0.26258308440446854
],
"md5": "0fe2fd08243d139606d5b50e1e16ed71"
},
{
"program": "reverse complement",
"size": "5 million",
"lang": "crow",
"runs": 5,
"median": 0.10683029145002365,
"min": 0.10193170793354511,
"max": 0.10796116665005684,
"times": [
0.10193170793354511,
0.10683029145002365,
0.10693079233169556,
0.10646695829927921,
0.10796116665005684
],
"md5": "1fb8dee22daca7e37c0728146f0922c3"
},
{
"program": "reverse complement",
"size": "5 million",
"lang": "rust",
"runs": 5,
"median": 0.06125212460756302,
"min": 0.06003108248114586,
"max": 0.0624286662787199,
"times": [
0.06003108248114586,
0.06125212460756302,
0.0624286662787199,
0.06237270962446928,
0.060746874660253525
],
"md5": "1fb8dee22daca7e37c0728146f0922c3"
},
{
"program": "reverse complement",
"size": "5 million",
"lang": "go",
"runs": 5,
"median": 0.0721376659348607,
"min": 0.07046554051339626,
"max": 0.07426141761243343,
"times": [
0.07050850056111813,
0.07271674927324057,
0.0721376659348607,
0.07426141761243343,
0.07046554051339626
],
"md5": "1fb8dee22daca7e37c0728146f0922c3"
},
{
"program": "reverse complement",
"size": "5 million",
"lang": "cpython",
"runs": 5,
"median": 0.23751424998044968,
"min": 0.23499695863574743,
"max": 0.24063954036682844,
"times": [
0.24056758359074593,
0.23751250002533197,
0.23499695863574743,
0.24063954036682844,
0.23751424998044968
],
"md5": "1fb8dee22daca7e37c0728146f0922c3"
},
{
"program": "reverse complement",
"size": "5 million",
"lang": "crow -i",
"runs": 3,
"median": 0.5537444585934281,
"min": 0.5533660007640719,
"max": 0.5548450835049152,
"times": [
0.5537444585934281,
0.5533660007640719,
0.5548450835049152
],
"md5": "1fb8dee22daca7e37c0728146f0922c3"
},
{
"program": "pidigits",
"size": "5,000 digits",
"lang": "crow",
"runs": 5,
"median": 0.30799200013279915,
"min": 0.3053967924788594,
"max": 0.3091033324599266,
"times": [
0.3091033324599266,
0.3053967924788594,
0.30762104224413633,
0.30887812469154596,
0.30799200013279915
],
"md5": "37df65c33d3422c721a1987269d50eb9"
},
{
"program": "pidigits",
"size": "5,000 digits",
"lang": "rust",
"runs": 5,
"median": 0.8134558331221342,
"min": 0.8106924993917346,
"max": 0.8214566251263022,
"times": [
0.8214566251263022,
0.8134558331221342,
0.8142017079517245,
0.8128723343834281,
0.8106924993917346
],
"md5": "37df65c33d3422c721a1987269d50eb9"
},
{
"program": "pidigits",
"size": "5,000 digits",
"lang": "go",
"runs": 5,
"median": 0.26464737486094236,
"min": 0.264099458232522,
"max": 0.26621820870786905,
"times": [
0.26588650047779083,
0.26464737486094236,
0.264099458232522,
0.2644134173169732,
0.26621820870786905
],
"md5": "37df65c33d3422c721a1987269d50eb9"
},
{
"program": "pidigits",
"size": "5,000 digits",
"lang": "cpython",
"runs": 5,
"median": 0.9903083741664886,
"min": 0.986634916625917,
"max": 0.9920188747346401,
"times": [
0.9920188747346401,
0.9889799160882831,
0.9903083741664886,
0.986634916625917,
0.9918315410614014
],
"md5": "37df65c33d3422c721a1987269d50eb9"
},
{
"program": "pidigits",
"size": "5,000 digits",
"lang": "crow -i",
"runs": 3,
"median": 0.3420087918639183,
"min": 0.34037216752767563,
"max": 0.3424863750115037,
"times": [
0.3420087918639183,
0.34037216752767563,
0.3424863750115037
],
"md5": "37df65c33d3422c721a1987269d50eb9"
},
{
"program": "binary trees",
"size": "depth 18",
"lang": "crow",
"runs": 5,
"median": 0.09356691595166922,
"min": 0.09192362520843744,
"max": 0.10606408305466175,
"times": [
0.09192362520843744,
0.09448933321982622,
0.10606408305466175,
0.09356691595166922,
0.09253508318215609
],
"md5": "04083ec3542623e5258db423ba0d7ee5"
},
{
"program": "binary trees",
"size": "depth 18",
"lang": "rust",
"runs": 5,
"median": 0.7552614593878388,
"min": 0.7305715000256896,
"max": 0.7952742492780089,
"times": [
0.7358365003019571,
0.7552614593878388,
0.7305715000256896,
0.7569287922233343,
0.7952742492780089
],
"md5": "04083ec3542623e5258db423ba0d7ee5"
},
{
"program": "binary trees",
"size": "depth 18",
"lang": "rust (mimalloc)",
"runs": 5,
"median": 0.1521114595234394,
"min": 0.15104408282786608,
"max": 0.15329466573894024,
"times": [
0.15324412565678358,
0.15329466573894024,
0.15168491657823324,
0.1521114595234394,
0.15104408282786608
],
"md5": "04083ec3542623e5258db423ba0d7ee5"
},
{
"program": "binary trees",
"size": "depth 18",
"lang": "go",
"runs": 5,
"median": 0.6415721671655774,
"min": 0.6383472923189402,
"max": 0.6574915423989296,
"times": [
0.6431570407003164,
0.6574915423989296,
0.6415721671655774,
0.6403929581865668,
0.6383472923189402
],
"md5": "04083ec3542623e5258db423ba0d7ee5"
},
{
"program": "binary trees",
"size": "depth 18",
"lang": "cpython",
"runs": 5,
"median": 0.6356999585404992,
"min": 0.6260153753682971,
"max": 0.6459652502089739,
"times": [
0.641351250000298,
0.6459652502089739,
0.6324360007420182,
0.6356999585404992,
0.6260153753682971
],
"md5": "04083ec3542623e5258db423ba0d7ee5"
},
{
"program": "binary trees",
"size": "depth 18",
"lang": "crow -i",
"runs": 3,
"median": 2.2552013341337442,
"min": 2.2539211669936776,
"max": 2.2572009582072496,
"times": [
2.2552013341337442,
2.2572009582072496,
2.2539211669936776
],
"md5": "04083ec3542623e5258db423ba0d7ee5"
}
]
Fannkuch-redux
fannkuch.co
--- fannkuch-redux: checksum and max flips over all permutations of 1..n
--- (benchmarks game), single-threaded
fn fannkuch[n]
perm1 = range[n]
perm = range[n]
count = range[n]
max_flips = 0
checksum = 0
perm_count = 0
r = n
while true
while r gt 1
count[r - 1] = r
r = r - 1
for i in range[n]
perm[i] = perm1[i]
flips = 0
k = perm[0]
while k gt 0
i = 0
j = k
while i lt j
t = perm[i]
perm[i] = perm[j]
perm[j] = t
i = i + 1
j = j - 1
flips = flips + 1
k = perm[0]
if flips gt max_flips
max_flips = flips
if perm_count % 2 eq 0
checksum = checksum + flips
else
checksum = checksum - flips
--- next permutation
while true
if r eq n
return [checksum, max_flips]
p0 = perm1[0]
for i in range[r]
perm1[i] = perm1[i + 1]
perm1[r] = p0
count[r] = count[r] - 1
if count[r] gt 0
break
r = r + 1
perm_count = perm_count + 1
res = fannkuch[10]
log['`res[0]`']
log['Pfannkuchen(10) = `res[1]`']
fannkuch.rs
fn fannkuch(n: usize) -> (i64, i64) {
let mut perm1: Vec<usize> = (0..n).collect();
let mut perm = vec![0usize; n];
let mut count = vec![0usize; n];
let (mut max_flips, mut checksum, mut perm_count) = (0i64, 0i64, 0i64);
let mut r = n;
loop {
while r > 1 {
count[r - 1] = r;
r -= 1;
}
perm.copy_from_slice(&perm1);
let mut flips = 0i64;
let mut k = perm[0];
while k > 0 {
perm[..=k].reverse();
flips += 1;
k = perm[0];
}
if flips > max_flips {
max_flips = flips;
}
checksum += if perm_count % 2 == 0 { flips } else { -flips };
loop {
if r == n {
return (checksum, max_flips);
}
let p0 = perm1[0];
for i in 0..r {
perm1[i] = perm1[i + 1];
}
perm1[r] = p0;
count[r] -= 1;
if count[r] > 0 {
break;
}
r += 1;
}
perm_count += 1;
}
}
fn main() {
let n: usize = std::env::args().nth(1).unwrap().parse().unwrap();
let (c, m) = fannkuch(n);
println!("{}\nPfannkuchen({}) = {}", c, n, m);
}
fannkuch.go
package main
import (
"fmt"
"os"
"strconv"
)
func fannkuch(n int) (int, int) {
perm1 := make([]int, n)
for i := range perm1 {
perm1[i] = i
}
perm := make([]int, n)
count := make([]int, n)
maxFlips, checksum, permCount := 0, 0, 0
r := n
for {
for r > 1 {
count[r-1] = r
r--
}
copy(perm, perm1)
flips := 0
for k := perm[0]; k > 0; k = perm[0] {
for i, j := 0, k; i < j; i, j = i+1, j-1 {
perm[i], perm[j] = perm[j], perm[i]
}
flips++
}
if flips > maxFlips {
maxFlips = flips
}
if permCount%2 == 0 {
checksum += flips
} else {
checksum -= flips
}
for {
if r == n {
return checksum, maxFlips
}
p0 := perm1[0]
for i := 0; i < r; i++ {
perm1[i] = perm1[i+1]
}
perm1[r] = p0
count[r]--
if count[r] > 0 {
break
}
r++
}
permCount++
}
}
func main() {
n, _ := strconv.Atoi(os.Args[1])
c, m := fannkuch(n)
fmt.Printf("%d\nPfannkuchen(%d) = %d\n", c, n, m)
}
fannkuch.py
import sys
def fannkuch(n):
perm1 = list(range(n))
count = [0] * n
max_flips = checksum = perm_count = 0
r = n
while True:
while r > 1:
count[r - 1] = r
r -= 1
perm = perm1[:]
flips = 0
k = perm[0]
while k:
perm[:k + 1] = perm[k::-1]
flips += 1
k = perm[0]
if flips > max_flips:
max_flips = flips
checksum += flips if perm_count % 2 == 0 else -flips
while True:
if r == n:
return checksum, max_flips
p0 = perm1[0]
perm1[:r] = perm1[1:r + 1]
perm1[r] = p0
count[r] -= 1
if count[r] > 0:
break
r += 1
perm_count += 1
n = int(sys.argv[1])
c, m = fannkuch(n)
print(f"{c}\nPfannkuchen({n}) = {m}")
N-body
nbody.co
--- n-body (benchmarks game): the Sun and four giant planets, advanced n
--- steps of 0.01; prints the energy before and after
pi = 3.141592653589793
solar_mass = 4.0 * pi * pi
days = 365.24
x = [0.0, 4.84143144246472090, 8.34336671824457987, 12.8943695621391310, 15.3796971148509165]
y = [0.0, -1.16032004402742839, 4.12479856412430479, -15.1111514016986312, -25.9193146099879641]
z = [0.0, -0.103622044471123109, -0.403523417114321381, -0.223307578892655734, 0.179258772950371181]
vx = [0.0, 0.00166007664274403694 * days, -0.00276742510726862411 * days, 0.00296460137564761618 * days, 0.00268067772490389322 * days]
vy = [0.0, 0.00769901118419740425 * days, 0.00499852801234917238 * days, 0.00237847173959480950 * days, 0.00162824170038242295 * days]
vz = [0.0, -0.0000690460016972063023 * days, 0.0000230417297573763929 * days, -0.0000296589568540237556 * days, -0.0000951592254519715870 * days]
m = [solar_mass, 0.000954791938424326609 * solar_mass, 0.000285885980666130812 * solar_mass, 0.0000436624404335156298 * solar_mass, 0.0000515138902046611451 * solar_mass]
fn energy[x, y, z, vx, vy, vz, m]
e = 0.0
for i in range[5]
e = e + 0.5 * m[i] * (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i])
j = i + 1
while j lt 5
dx = x[i] - x[j]
dy = y[i] - y[j]
dz = z[i] - z[j]
e = e - m[i] * m[j] / sqrt[dx * dx + dy * dy + dz * dz]
j = j + 1
return e
fn run[n, x, y, z, vx, vy, vz, m]
dt = 0.01
step = 0
while step lt n
i = 0
while i lt 5
j = i + 1
while j lt 5
dx = x[i] - x[j]
dy = y[i] - y[j]
dz = z[i] - z[j]
d2 = dx * dx + dy * dy + dz * dz
mag = dt / (d2 * sqrt[d2])
vx[i] = vx[i] - dx * m[j] * mag
vy[i] = vy[i] - dy * m[j] * mag
vz[i] = vz[i] - dz * m[j] * mag
vx[j] = vx[j] + dx * m[i] * mag
vy[j] = vy[j] + dy * m[i] * mag
vz[j] = vz[j] + dz * m[i] * mag
j = j + 1
i = i + 1
for k in range[5]
x[k] = x[k] + dt * vx[k]
y[k] = y[k] + dt * vy[k]
z[k] = z[k] + dt * vz[k]
step = step + 1
return [energy[x, y, z, vx, vy, vz, m]]
--- the Sun's velocity cancels the system's momentum
px = 0.0
py = 0.0
pz = 0.0
for k in range[5]
px = px + vx[k] * m[k]
py = py + vy[k] * m[k]
pz = pz + vz[k] * m[k]
vx[0] = -px / solar_mass
vy[0] = -py / solar_mass
vz[0] = -pz / solar_mass
n = 50000000
log[fmt[energy[x, y, z, vx, vy, vz, m], 9]]
after = run[n, x, y, z, vx, vy, vz, m]
log[fmt[after[0], 9]]
nbody_1m.co
--- n-body (benchmarks game): the Sun and four giant planets, advanced n
--- steps of 0.01; prints the energy before and after
pi = 3.141592653589793
solar_mass = 4.0 * pi * pi
days = 365.24
x = [0.0, 4.84143144246472090, 8.34336671824457987, 12.8943695621391310, 15.3796971148509165]
y = [0.0, -1.16032004402742839, 4.12479856412430479, -15.1111514016986312, -25.9193146099879641]
z = [0.0, -0.103622044471123109, -0.403523417114321381, -0.223307578892655734, 0.179258772950371181]
vx = [0.0, 0.00166007664274403694 * days, -0.00276742510726862411 * days, 0.00296460137564761618 * days, 0.00268067772490389322 * days]
vy = [0.0, 0.00769901118419740425 * days, 0.00499852801234917238 * days, 0.00237847173959480950 * days, 0.00162824170038242295 * days]
vz = [0.0, -0.0000690460016972063023 * days, 0.0000230417297573763929 * days, -0.0000296589568540237556 * days, -0.0000951592254519715870 * days]
m = [solar_mass, 0.000954791938424326609 * solar_mass, 0.000285885980666130812 * solar_mass, 0.0000436624404335156298 * solar_mass, 0.0000515138902046611451 * solar_mass]
fn energy[x, y, z, vx, vy, vz, m]
e = 0.0
for i in range[5]
e = e + 0.5 * m[i] * (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i])
j = i + 1
while j lt 5
dx = x[i] - x[j]
dy = y[i] - y[j]
dz = z[i] - z[j]
e = e - m[i] * m[j] / sqrt[dx * dx + dy * dy + dz * dz]
j = j + 1
return e
fn run[n, x, y, z, vx, vy, vz, m]
dt = 0.01
step = 0
while step lt n
i = 0
while i lt 5
j = i + 1
while j lt 5
dx = x[i] - x[j]
dy = y[i] - y[j]
dz = z[i] - z[j]
d2 = dx * dx + dy * dy + dz * dz
mag = dt / (d2 * sqrt[d2])
vx[i] = vx[i] - dx * m[j] * mag
vy[i] = vy[i] - dy * m[j] * mag
vz[i] = vz[i] - dz * m[j] * mag
vx[j] = vx[j] + dx * m[i] * mag
vy[j] = vy[j] + dy * m[i] * mag
vz[j] = vz[j] + dz * m[i] * mag
j = j + 1
i = i + 1
for k in range[5]
x[k] = x[k] + dt * vx[k]
y[k] = y[k] + dt * vy[k]
z[k] = z[k] + dt * vz[k]
step = step + 1
return [energy[x, y, z, vx, vy, vz, m]]
--- the Sun's velocity cancels the system's momentum
px = 0.0
py = 0.0
pz = 0.0
for k in range[5]
px = px + vx[k] * m[k]
py = py + vy[k] * m[k]
pz = pz + vz[k] * m[k]
vx[0] = -px / solar_mass
vy[0] = -py / solar_mass
vz[0] = -pz / solar_mass
n = 1000000
log[fmt[energy[x, y, z, vx, vy, vz, m], 9]]
after = run[n, x, y, z, vx, vy, vz, m]
log[fmt[after[0], 9]]
nbody.rs
const PI: f64 = 3.141592653589793;
const SOLAR_MASS: f64 = 4.0 * PI * PI;
const DAYS: f64 = 365.24;
#[derive(Clone, Copy)]
struct Body {
x: [f64; 3],
v: [f64; 3],
m: f64,
}
fn bodies() -> [Body; 5] {
[
Body { x: [0.0; 3], v: [0.0; 3], m: SOLAR_MASS },
Body {
x: [4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01],
v: [1.66007664274403694e-03 * DAYS, 7.69901118419740425e-03 * DAYS, -6.90460016972063023e-05 * DAYS],
m: 9.54791938424326609e-04 * SOLAR_MASS,
},
Body {
x: [8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01],
v: [-2.76742510726862411e-03 * DAYS, 4.99852801234917238e-03 * DAYS, 2.30417297573763929e-05 * DAYS],
m: 2.85885980666130812e-04 * SOLAR_MASS,
},
Body {
x: [1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01],
v: [2.96460137564761618e-03 * DAYS, 2.37847173959480950e-03 * DAYS, -2.96589568540237556e-05 * DAYS],
m: 4.36624404335156298e-05 * SOLAR_MASS,
},
Body {
x: [1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01],
v: [2.68067772490389322e-03 * DAYS, 1.62824170038242295e-03 * DAYS, -9.51592254519715870e-05 * DAYS],
m: 5.15138902046611451e-05 * SOLAR_MASS,
},
]
}
fn offset_momentum(b: &mut [Body; 5]) {
let mut p = [0.0f64; 3];
for body in b.iter() {
for k in 0..3 {
p[k] += body.v[k] * body.m;
}
}
for k in 0..3 {
b[0].v[k] = -p[k] / SOLAR_MASS;
}
}
fn energy(b: &[Body; 5]) -> f64 {
let mut e = 0.0;
for i in 0..5 {
e += 0.5 * b[i].m * (b[i].v[0] * b[i].v[0] + b[i].v[1] * b[i].v[1] + b[i].v[2] * b[i].v[2]);
for j in i + 1..5 {
let dx = b[i].x[0] - b[j].x[0];
let dy = b[i].x[1] - b[j].x[1];
let dz = b[i].x[2] - b[j].x[2];
e -= b[i].m * b[j].m / (dx * dx + dy * dy + dz * dz).sqrt();
}
}
e
}
fn advance(b: &mut [Body; 5], dt: f64) {
for i in 0..5 {
for j in i + 1..5 {
let dx = b[i].x[0] - b[j].x[0];
let dy = b[i].x[1] - b[j].x[1];
let dz = b[i].x[2] - b[j].x[2];
let d2 = dx * dx + dy * dy + dz * dz;
let mag = dt / (d2 * d2.sqrt());
let (mi, mj) = (b[i].m, b[j].m);
b[i].v[0] -= dx * mj * mag;
b[i].v[1] -= dy * mj * mag;
b[i].v[2] -= dz * mj * mag;
b[j].v[0] += dx * mi * mag;
b[j].v[1] += dy * mi * mag;
b[j].v[2] += dz * mi * mag;
}
}
for body in b.iter_mut() {
for k in 0..3 {
body.x[k] += dt * body.v[k];
}
}
}
fn main() {
let n: usize = std::env::args().nth(1).unwrap().parse().unwrap();
let mut b = bodies();
offset_momentum(&mut b);
println!("{:.9}", energy(&b));
for _ in 0..n {
advance(&mut b, 0.01);
}
println!("{:.9}", energy(&b));
}
nbody.go
package main
import (
"fmt"
"math"
"os"
"strconv"
)
const solarMass = 4 * math.Pi * math.Pi
const days = 365.24
type body struct{ x, y, z, vx, vy, vz, m float64 }
func bodies() []body {
return []body{
{0, 0, 0, 0, 0, 0, solarMass},
{4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01,
1.66007664274403694e-03 * days, 7.69901118419740425e-03 * days, -6.90460016972063023e-05 * days,
9.54791938424326609e-04 * solarMass},
{8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01,
-2.76742510726862411e-03 * days, 4.99852801234917238e-03 * days, 2.30417297573763929e-05 * days,
2.85885980666130812e-04 * solarMass},
{1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01,
2.96460137564761618e-03 * days, 2.37847173959480950e-03 * days, -2.96589568540237556e-05 * days,
4.36624404335156298e-05 * solarMass},
{1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01,
2.68067772490389322e-03 * days, 1.62824170038242295e-03 * days, -9.51592254519715870e-05 * days,
5.15138902046611451e-05 * solarMass},
}
}
func offsetMomentum(b []body) {
var px, py, pz float64
for _, o := range b {
px += o.vx * o.m
py += o.vy * o.m
pz += o.vz * o.m
}
b[0].vx, b[0].vy, b[0].vz = -px/solarMass, -py/solarMass, -pz/solarMass
}
func energy(b []body) float64 {
e := 0.0
for i := range b {
e += 0.5 * b[i].m * (b[i].vx*b[i].vx + b[i].vy*b[i].vy + b[i].vz*b[i].vz)
for j := i + 1; j < len(b); j++ {
dx, dy, dz := b[i].x-b[j].x, b[i].y-b[j].y, b[i].z-b[j].z
e -= b[i].m * b[j].m / math.Sqrt(dx*dx+dy*dy+dz*dz)
}
}
return e
}
func advance(b []body, dt float64) {
for i := range b {
for j := i + 1; j < len(b); j++ {
dx, dy, dz := b[i].x-b[j].x, b[i].y-b[j].y, b[i].z-b[j].z
d2 := dx*dx + dy*dy + dz*dz
mag := dt / (d2 * math.Sqrt(d2))
mi, mj := b[i].m, b[j].m
b[i].vx -= dx * mj * mag
b[i].vy -= dy * mj * mag
b[i].vz -= dz * mj * mag
b[j].vx += dx * mi * mag
b[j].vy += dy * mi * mag
b[j].vz += dz * mi * mag
}
}
for i := range b {
b[i].x += dt * b[i].vx
b[i].y += dt * b[i].vy
b[i].z += dt * b[i].vz
}
}
func main() {
n, _ := strconv.Atoi(os.Args[1])
b := bodies()
offsetMomentum(b)
fmt.Printf("%.9f\n", energy(b))
for i := 0; i < n; i++ {
advance(b, 0.01)
}
fmt.Printf("%.9f\n", energy(b))
}
nbody.py
import sys
from math import pi, sqrt
SOLAR_MASS = 4 * pi * pi
DAYS = 365.24
# x, y, z, vx, vy, vz, mass
BODIES = [
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, SOLAR_MASS],
[4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01,
1.66007664274403694e-03 * DAYS, 7.69901118419740425e-03 * DAYS, -6.90460016972063023e-05 * DAYS,
9.54791938424326609e-04 * SOLAR_MASS],
[8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01,
-2.76742510726862411e-03 * DAYS, 4.99852801234917238e-03 * DAYS, 2.30417297573763929e-05 * DAYS,
2.85885980666130812e-04 * SOLAR_MASS],
[1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01,
2.96460137564761618e-03 * DAYS, 2.37847173959480950e-03 * DAYS, -2.96589568540237556e-05 * DAYS,
4.36624404335156298e-05 * SOLAR_MASS],
[1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01,
2.68067772490389322e-03 * DAYS, 1.62824170038242295e-03 * DAYS, -9.51592254519715870e-05 * DAYS,
5.15138902046611451e-05 * SOLAR_MASS],
]
def offset_momentum(b):
px = sum(o[3] * o[6] for o in b)
py = sum(o[4] * o[6] for o in b)
pz = sum(o[5] * o[6] for o in b)
b[0][3], b[0][4], b[0][5] = -px / SOLAR_MASS, -py / SOLAR_MASS, -pz / SOLAR_MASS
def energy(b):
e = 0.0
for i in range(len(b)):
bi = b[i]
e += 0.5 * bi[6] * (bi[3] * bi[3] + bi[4] * bi[4] + bi[5] * bi[5])
for j in range(i + 1, len(b)):
bj = b[j]
dx, dy, dz = bi[0] - bj[0], bi[1] - bj[1], bi[2] - bj[2]
e -= bi[6] * bj[6] / sqrt(dx * dx + dy * dy + dz * dz)
return e
def run(n, dt=0.01):
pairs = [(BODIES[i], BODIES[j]) for i in range(5) for j in range(i + 1, 5)]
for _ in range(n):
for a, c in pairs:
dx, dy, dz = a[0] - c[0], a[1] - c[1], a[2] - c[2]
d2 = dx * dx + dy * dy + dz * dz
mag = dt / (d2 * sqrt(d2))
ma, mc = a[6], c[6]
a[3] -= dx * mc * mag
a[4] -= dy * mc * mag
a[5] -= dz * mc * mag
c[3] += dx * ma * mag
c[4] += dy * ma * mag
c[5] += dz * ma * mag
for o in BODIES:
o[0] += dt * o[3]
o[1] += dt * o[4]
o[2] += dt * o[5]
n = int(sys.argv[1])
offset_momentum(BODIES)
print("%.9f" % energy(BODIES))
run(n)
print("%.9f" % energy(BODIES))
Spectral norm
spectralnorm.co
--- spectral norm (benchmarks game): the largest singular value of the
--- infinite matrix A(i, j) = 1 / ((i+j)(i+j+1)/2 + i + 1), by power
--- iteration on A'A. Each product is a `for` over the output rows, so it
--- runs on every core; `cpu[1]` makes it single-threaded.
--- A(i, j)
fn a[i, j]
return 1.0 / float[(i + j) * (i + j + 1) / 2 + i + 1]
--- row i of A v
fn row_av[i, v, n]
s = 0.0
j = 0
while j lt n
s = s + a[i, j] * v[j]
j = j + 1
return s
--- row i of A' v
fn row_atv[i, v, n]
s = 0.0
j = 0
while j lt n
s = s + a[j, i] * v[j]
j = j + 1
return s
n = 5500
u = for i in range[n]
1.0
v = u
step = 0
while step lt 10
t = for i in range[n]
row_av[i, u, n]
v = for i in range[n]
row_atv[i, t, n]
t = for i in range[n]
row_av[i, v, n]
u = for i in range[n]
row_atv[i, t, n]
step = step + 1
vbv = 0.0
vv = 0.0
for i in range[n]
vbv = vbv + u[i] * v[i]
vv = vv + v[i] * v[i]
log[fmt[sqrt[vbv / vv], 9]]
spectralnorm.rs
// usage: spectralnorm N THREADS (THREADS = 1: single-threaded)
fn a(i: usize, j: usize) -> f64 {
1.0 / (((i + j) * (i + j + 1) / 2 + i + 1) as f64)
}
// out[i] = sum_j A(i,j) v[j] (transpose: A(j,i))
fn mul(out: &mut [f64], v: &[f64], transpose: bool, threads: usize) {
let n = v.len();
let row = |i: usize| -> f64 {
let mut s = 0.0;
for j in 0..n {
s += if transpose { a(j, i) } else { a(i, j) } * v[j];
}
s
};
if threads == 1 {
for i in 0..n {
out[i] = row(i);
}
} else {
let chunk = (n + threads - 1) / threads;
std::thread::scope(|s| {
for (c, o) in out.chunks_mut(chunk).enumerate() {
let row = &row;
s.spawn(move || {
for (k, oi) in o.iter_mut().enumerate() {
*oi = row(c * chunk + k);
}
});
}
});
}
}
fn main() {
let args: Vec<usize> = std::env::args().skip(1).map(|s| s.parse().unwrap()).collect();
let (n, threads) = (args[0], args[1]);
let mut u = vec![1.0; n];
let mut v = vec![0.0; n];
let mut t = vec![0.0; n];
for _ in 0..10 {
mul(&mut t, &u, false, threads);
mul(&mut v, &t, true, threads);
mul(&mut t, &v, false, threads);
mul(&mut u, &t, true, threads);
}
let (mut vbv, mut vv) = (0.0, 0.0);
for i in 0..n {
vbv += u[i] * v[i];
vv += v[i] * v[i];
}
println!("{:.9}", (vbv / vv).sqrt());
}
spectralnorm.go
package main
// usage: spectralnorm N THREADS (THREADS = 1: single-threaded)
import (
"fmt"
"math"
"os"
"strconv"
"sync"
)
func a(i, j int) float64 {
return 1.0 / float64((i+j)*(i+j+1)/2+i+1)
}
func mul(out, v []float64, transpose bool, threads int) {
n := len(v)
rows := func(lo, hi int) {
for i := lo; i < hi; i++ {
s := 0.0
for j := 0; j < n; j++ {
if transpose {
s += a(j, i) * v[j]
} else {
s += a(i, j) * v[j]
}
}
out[i] = s
}
}
if threads == 1 {
rows(0, n)
return
}
chunk := (n + threads - 1) / threads
var wg sync.WaitGroup
for lo := 0; lo < n; lo += chunk {
hi := min(lo+chunk, n)
wg.Add(1)
go func() { defer wg.Done(); rows(lo, hi) }()
}
wg.Wait()
}
func main() {
n, _ := strconv.Atoi(os.Args[1])
threads, _ := strconv.Atoi(os.Args[2])
u, v, t := make([]float64, n), make([]float64, n), make([]float64, n)
for i := range u {
u[i] = 1
}
for k := 0; k < 10; k++ {
mul(t, u, false, threads)
mul(v, t, true, threads)
mul(t, v, false, threads)
mul(u, t, true, threads)
}
var vbv, vv float64
for i := 0; i < n; i++ {
vbv += u[i] * v[i]
vv += v[i] * v[i]
}
fmt.Printf("%.9f\n", math.Sqrt(vbv/vv))
}
spectralnorm.py
# usage: spectralnorm.py N PROCS (PROCS = 1: single-process)
import sys
from math import sqrt
from multiprocessing import Pool
def rows(args):
lo, hi, v, transpose = args
n = len(v)
out = []
for i in range(lo, hi):
s = 0.0
if transpose:
for j in range(n):
s += v[j] / ((i + j) * (i + j + 1) // 2 + j + 1)
else:
for j in range(n):
s += v[j] / ((i + j) * (i + j + 1) // 2 + i + 1)
out.append(s)
return out
def mul(v, transpose, pool, procs):
n = len(v)
chunk = (n + procs - 1) // procs
jobs = [(lo, min(lo + chunk, n), v, transpose) for lo in range(0, n, chunk)]
parts = pool.map(rows, jobs) if pool else [rows(j) for j in jobs]
return [x for p in parts for x in p]
if __name__ == "__main__":
n, procs = int(sys.argv[1]), int(sys.argv[2])
pool = Pool(procs) if procs > 1 else None
u = [1.0] * n
for _ in range(10):
v = mul(mul(u, False, pool, procs), True, pool, procs)
u = mul(mul(v, False, pool, procs), True, pool, procs)
vbv = vv = 0.0
for i in range(n):
vbv += u[i] * v[i]
vv += v[i] * v[i]
print("%.9f" % sqrt(vbv / vv))
K-nucleotide
knucleotide.co
--- k-nucleotide (benchmarks game): the frequencies of the 1- and 2-mers of
--- section THREE of a fasta file, and the counts of five fixed k-mers.
---
--- The section's lines joined are one text, and every count is `occurs` of
--- it — a search at the speed of reading the text. A pattern that cannot
--- overlap itself (each of the five, and every pair of two different
--- letters) occurs as often as `occurs` counts it; a doubled letter XX
--- overlaps itself, so its count is the X's that have a next letter less
--- the X's followed by another letter. Each list of counts is a `for`: the
--- counts run on every core, and `cpu[1]` makes it single-threaded.
use io
--- the names by count, largest first, ties by name; one 'name percent' a line
fn report[names, counts, total]
n = len[names]
used = for i in range[n]
false
done = 0
while done lt n
best = -1
j = 0
while j lt n
if not used[j]
if best lt 0
best = j
else
if counts[j] gt counts[best]
best = j
j = j + 1
used[best] = true
log['`names[best]` `fmt[100.0 * float[counts[best]] / float[total], 3]`']
done = done + 1
log['']
return 0
rows = lines['fasta.txt']
first = 0
while not starts[rows[first], '>THREE']
first = first + 1
stop = len[rows]
if rows[stop - 1] eq ''
stop = stop - 1
seq = upper[join[slice[rows, first + 1, stop], '']]
n = len[seq]
base = ['A', 'C', 'G', 'T']
pairs = ['AA', 'AC', 'AG', 'AT', 'CA', 'CC', 'CG', 'CT', 'GA', 'GC', 'GG', 'GT', 'TA', 'TC', 'TG', 'TT']
pats = ['GGT', 'GGTA', 'GGTATT', 'GGTATTTTAATT', 'GGTATTTTAATTTATAGT']
ones = for b in base
occurs[seq, b]
apart = for p in pairs
occurs[seq, p]
found = for p in pats
occurs[seq, p]
--- a doubled letter: the X's with a next letter, less those followed by
--- another letter
last = seq[n - 1]
twos = for k in range[16]
bx = k / 4
by = k % 4
if bx ne by
apart[k]
else
others = 0
for o in range[4]
if o ne bx
others = others + apart[4 * bx + o]
trailing = 0
if last eq base[bx]
trailing = 1
ones[bx] - trailing - others
one = report[base, ones, n]
two = report[pairs, twos, n - 1]
for t in range[5]
log['`found[t]`\t`pats[t]`']
knucleotide.rs
// k-nucleotide (benchmarks game): read fasta on stdin, take section THREE and
// report the frequencies of 1- and 2-mers and the counts of five fixed k-mers.
// One pass of rolling 2-bit codes per chunk of the sequence; chunks overlap by
// 17 bases of warm-up so no window is lost, and are counted on THREADS workers.
// usage: knucleotide THREADS < fasta
use std::io::Read;
use std::sync::atomic::{AtomicUsize, Ordering};
const PATTERNS: [&str; 5] = ["GGT", "GGTA", "GGTATT", "GGTATTTTAATT", "GGTATTTTAATTTATAGT"];
const CHUNK: usize = 1 << 20;
fn code(s: &str) -> u64 {
s.bytes().fold(0, |c, b| c << 2 | b2c(b) as u64)
}
fn b2c(b: u8) -> u8 {
match b {
b'A' => 0,
b'C' => 1,
b'G' => 2,
_ => 3,
}
}
#[derive(Default, Clone)]
struct Counts {
c1: [u64; 4],
c2: [u64; 16],
pat: [u64; 5],
}
fn count_chunk(seq: &[u8], s: usize, e: usize, targets: &[(usize, u64); 5]) -> Counts {
let mut r = Counts::default();
let mut code = 0u64;
for j in s.saturating_sub(17)..e {
code = (code << 2 | seq[j] as u64) & ((1u64 << 36) - 1);
if j < s {
continue;
}
let have = j + 1;
r.c1[(code & 3) as usize] += 1;
if have >= 2 {
r.c2[(code & 15) as usize] += 1;
}
for (t, &(k, tc)) in targets.iter().enumerate() {
if have >= k && code & ((1u64 << (2 * k)) - 1) == tc {
r.pat[t] += 1;
}
}
}
r
}
fn report(names: &[String], counts: &[u64], total: u64) {
let mut idx: Vec<usize> = (0..names.len()).collect();
idx.sort_by(|&a, &b| counts[b].cmp(&counts[a]).then(names[a].cmp(&names[b])));
for i in idx {
println!("{} {:.3}", names[i], 100.0 * counts[i] as f64 / total as f64);
}
println!();
}
fn main() {
let threads: usize = std::env::args().nth(1).unwrap().parse().unwrap();
let mut input = Vec::new();
std::io::stdin().read_to_end(&mut input).unwrap();
let mut seq: Vec<u8> = Vec::new();
let mut in_three = false;
for line in input.split(|&b| b == b'\n') {
if line.starts_with(b">") {
in_three = line.starts_with(b">THREE");
} else if in_three {
seq.extend(line.iter().map(|&b| b2c(b.to_ascii_uppercase())));
}
}
let n = seq.len();
let targets: [(usize, u64); 5] = std::array::from_fn(|i| (PATTERNS[i].len(), code(PATTERNS[i])));
let nchunks = (n + CHUNK - 1) / CHUNK;
let next = AtomicUsize::new(0);
let total = std::thread::scope(|s| {
let hs: Vec<_> = (0..threads)
.map(|_| {
s.spawn(|| {
let mut acc = Counts::default();
loop {
let c = next.fetch_add(1, Ordering::Relaxed);
if c >= nchunks {
return acc;
}
let r = count_chunk(&seq, c * CHUNK, ((c + 1) * CHUNK).min(n), &targets);
for i in 0..4 { acc.c1[i] += r.c1[i]; }
for i in 0..16 { acc.c2[i] += r.c2[i]; }
for i in 0..5 { acc.pat[i] += r.pat[i]; }
}
})
})
.collect();
let mut t = Counts::default();
for h in hs {
let r = h.join().unwrap();
for i in 0..4 { t.c1[i] += r.c1[i]; }
for i in 0..16 { t.c2[i] += r.c2[i]; }
for i in 0..5 { t.pat[i] += r.pat[i]; }
}
t
});
let l1: Vec<String> = "ACGT".chars().map(|c| c.to_string()).collect();
let l2: Vec<String> = "ACGT".chars().flat_map(|a| "ACGT".chars().map(move |b| format!("{}{}", a, b))).collect();
report(&l1, &total.c1, n as u64);
report(&l2, &total.c2, n as u64 - 1);
for (i, p) in PATTERNS.iter().enumerate() {
println!("{}\t{}", total.pat[i], p);
}
}
knucleotide.go
package main
// k-nucleotide (benchmarks game): read fasta on stdin, take section THREE and
// report the frequencies of 1- and 2-mers and the counts of five fixed k-mers.
// One pass of rolling 2-bit codes per chunk of the sequence; chunks overlap by
// 17 bases of warm-up so no window is lost, and are counted on THREADS workers.
// usage: knucleotide THREADS < fasta
import (
"bytes"
"fmt"
"io"
"os"
"sort"
"strconv"
"sync"
"sync/atomic"
)
var patterns = []string{"GGT", "GGTA", "GGTATT", "GGTATTTTAATT", "GGTATTTTAATTTATAGT"}
const chunk = 1 << 20
func b2c(b byte) byte {
switch b {
case 'A', 'a':
return 0
case 'C', 'c':
return 1
case 'G', 'g':
return 2
}
return 3
}
type counts struct {
c1 [4]uint64
c2 [16]uint64
pat [5]uint64
}
func (a *counts) add(b *counts) {
for i := range a.c1 {
a.c1[i] += b.c1[i]
}
for i := range a.c2 {
a.c2[i] += b.c2[i]
}
for i := range a.pat {
a.pat[i] += b.pat[i]
}
}
type target struct {
k int
code uint64
}
func countChunk(seq []byte, s, e int, targets []target) counts {
var r counts
var code uint64
start := s - 17
if start < 0 {
start = 0
}
for j := start; j < e; j++ {
code = (code<<2 | uint64(seq[j])) & (1<<36 - 1)
if j < s {
continue
}
have := j + 1
r.c1[code&3]++
if have >= 2 {
r.c2[code&15]++
}
for t, tg := range targets {
if have >= tg.k && code&(1<<(2*uint(tg.k))-1) == tg.code {
r.pat[t]++
}
}
}
return r
}
func report(names []string, cnt []uint64, total uint64) {
idx := make([]int, len(names))
for i := range idx {
idx[i] = i
}
sort.Slice(idx, func(a, b int) bool {
if cnt[idx[a]] != cnt[idx[b]] {
return cnt[idx[a]] > cnt[idx[b]]
}
return names[idx[a]] < names[idx[b]]
})
for _, i := range idx {
fmt.Printf("%s %.3f\n", names[i], 100.0*float64(cnt[i])/float64(total))
}
fmt.Println()
}
func main() {
threads, _ := strconv.Atoi(os.Args[1])
input, _ := io.ReadAll(os.Stdin)
seq := make([]byte, 0, len(input))
inThree := false
for _, line := range bytes.Split(input, []byte{'\n'}) {
if len(line) > 0 && line[0] == '>' {
inThree = bytes.HasPrefix(line, []byte(">THREE"))
} else if inThree {
for _, b := range line {
seq = append(seq, b2c(b))
}
}
}
n := len(seq)
targets := make([]target, len(patterns))
for i, p := range patterns {
var c uint64
for j := 0; j < len(p); j++ {
c = c<<2 | uint64(b2c(p[j]))
}
targets[i] = target{len(p), c}
}
nchunks := (n + chunk - 1) / chunk
var next int64 = -1
res := make([]counts, threads)
var wg sync.WaitGroup
for w := 0; w < threads; w++ {
wg.Add(1)
go func(w int) {
defer wg.Done()
for {
c := int(atomic.AddInt64(&next, 1))
if c >= nchunks {
return
}
r := countChunk(seq, c*chunk, min((c+1)*chunk, n), targets)
res[w].add(&r)
}
}(w)
}
wg.Wait()
var total counts
for w := range res {
total.add(&res[w])
}
l1 := []string{"A", "C", "G", "T"}
var l2 []string
for _, a := range l1 {
for _, b := range l1 {
l2 = append(l2, a+b)
}
}
report(l1, total.c1[:], uint64(n))
report(l2, total.c2[:], uint64(n-1))
for i, p := range patterns {
fmt.Printf("%d\t%s\n", total.pat[i], p)
}
}
knucleotide.py
# k-nucleotide (benchmarks game): read fasta on stdin, take section THREE and
# report the frequencies of 1- and 2-mers and the counts of five fixed k-mers.
# usage: knucleotide.py < fasta
import sys
from collections import Counter
PATTERNS = ["GGT", "GGTA", "GGTATT", "GGTATTTTAATT", "GGTATTTTAATTTATAGT"]
def read_three():
parts, in_three = [], False
for line in sys.stdin:
if line.startswith(">"):
in_three = line.startswith(">THREE")
elif in_three:
parts.append(line.strip().upper())
return "".join(parts)
def frequencies(seq, k):
counts = Counter(seq[i:i + k] for i in range(len(seq) - k + 1))
total = len(seq) - k + 1
for name, c in sorted(counts.items(), key=lambda kv: (-kv[1], kv[0])):
print("%s %.3f" % (name, 100.0 * c / total))
print()
def count_overlapping(seq, pat):
n, i = 0, seq.find(pat)
while i >= 0:
n += 1
i = seq.find(pat, i + 1)
return n
seq = read_three()
frequencies(seq, 1)
frequencies(seq, 2)
for p in PATTERNS:
print("%d\t%s" % (count_overlapping(seq, p), p))
Reverse complement
revcomp.co
--- reverse-complement (benchmarks game): each sequence of a fasta file,
--- reverse-complemented (IUB codes) and written in lines of 60.
---
--- A section's lines joined are one text: `translate` maps each code to its
--- complement in one pass, `reverse` turns the text round, and the answer
--- is cut into lines of 60, a `for` over them (they run on every core, and
--- `cpu[1]` makes it single-threaded).
use io
fn revcomp[body]
seq = reverse[translate[join[body, ''], 'ACGTUMRWSYKVHDBNacgtumrwsykvhdbn', 'TGCAAKYWSRMBDHVNTGCAAKYWSRMBDHVN']]
n = len[seq]
outs = for k in range[(n + 59) / 60]
slice[seq, k * 60, min[n, k * 60 + 60]]
return join[outs, '\n']
rows = lines['fasta.txt']
stop = len[rows]
if rows[stop - 1] eq ''
stop = stop - 1
--- the section headers
heads = []
i = 0
while i lt stop
if starts[rows[i], '>']
heads = push[heads, i]
i = i + 1
heads = push[heads, stop]
s = 0
while s lt len[heads] - 1
log[rows[heads[s]]]
log[revcomp[slice[rows, heads[s] + 1, heads[s + 1]]]]
s = s + 1
revcomp.rs
// reverse-complement (benchmarks game): read fasta on stdin, write each
// sequence reverse-complemented (IUB codes) in lines of 60.
// Output lines are filled independently, THREADS at a time.
// usage: revcomp THREADS < fasta > output
use std::io::{Read, Write};
use std::sync::Mutex;
const WIDTH: usize = 60;
const CHUNK_LINES: usize = 4096;
fn table() -> [u8; 256] {
let mut t = [0u8; 256];
for i in 0..256 {
t[i] = i as u8;
}
for (a, b) in "ACGTUMRWSYKVHDBN".bytes().zip("TGCAAKYWSRMBDHVN".bytes()) {
t[a as usize] = b;
t[a.to_ascii_lowercase() as usize] = b;
}
t
}
fn rc_section(seq: &[u8], t: &[u8; 256], threads: usize) -> Vec<u8> {
let n = seq.len();
let lines = (n + WIDTH - 1) / WIDTH;
let mut out = vec![0u8; n + lines];
let it = Mutex::new(out.chunks_mut(CHUNK_LINES * (WIDTH + 1)).enumerate());
std::thread::scope(|s| {
for _ in 0..threads {
s.spawn(|| loop {
let item = it.lock().unwrap().next();
let Some((c, buf)) = item else { break };
let mut o = 0;
let mut line = c * CHUNK_LINES;
while o < buf.len() {
let w = WIDTH.min(buf.len() - o - 1);
for k in 0..w {
buf[o + k] = t[seq[n - 1 - (line * WIDTH + k)] as usize];
}
buf[o + w] = b'\n';
o += w + 1;
line += 1;
}
});
}
});
out
}
fn main() {
let threads: usize = std::env::args().nth(1).unwrap().parse().unwrap();
let mut input = Vec::new();
std::io::stdin().read_to_end(&mut input).unwrap();
let t = table();
let out = std::io::stdout();
let mut out = out.lock();
let mut header: Option<&[u8]> = None;
let mut seq: Vec<u8> = Vec::new();
let mut flush = |header: &Option<&[u8]>, seq: &mut Vec<u8>, out: &mut std::io::StdoutLock| {
if let Some(h) = header {
out.write_all(h).unwrap();
out.write_all(b"\n").unwrap();
out.write_all(&rc_section(seq, &t, threads)).unwrap();
}
seq.clear();
};
for line in input.split(|&b| b == b'\n') {
if line.first() == Some(&b'>') {
flush(&header, &mut seq, &mut out);
header = Some(line);
} else {
seq.extend_from_slice(line);
}
}
flush(&header, &mut seq, &mut out);
}
revcomp.go
package main
// reverse-complement (benchmarks game): read fasta on stdin, write each
// sequence reverse-complemented (IUB codes) in lines of 60.
// Output lines are filled independently, THREADS at a time.
// usage: revcomp THREADS < fasta > output
import (
"bufio"
"bytes"
"io"
"os"
"strconv"
"sync"
"sync/atomic"
)
const width = 60
const chunkLines = 4096
func table() [256]byte {
var t [256]byte
for i := range t {
t[i] = byte(i)
}
from, to := "ACGTUMRWSYKVHDBN", "TGCAAKYWSRMBDHVN"
for i := 0; i < len(from); i++ {
t[from[i]] = to[i]
t[from[i]+32] = to[i]
}
return t
}
func rcSection(seq []byte, t *[256]byte, threads int) []byte {
n := len(seq)
lines := (n + width - 1) / width
out := make([]byte, n+lines)
step := chunkLines * (width + 1)
nchunks := (len(out) + step - 1) / step
var next int64 = -1
var wg sync.WaitGroup
for w := 0; w < threads; w++ {
wg.Add(1)
go func() {
defer wg.Done()
for {
c := int(atomic.AddInt64(&next, 1))
if c >= nchunks {
return
}
buf := out[c*step:min((c+1)*step, len(out))]
line := c * chunkLines
for o := 0; o < len(buf); {
w := min(width, len(buf)-o-1)
for k := 0; k < w; k++ {
buf[o+k] = t[seq[n-1-(line*width+k)]]
}
buf[o+w] = '\n'
o += w + 1
line++
}
}
}()
}
wg.Wait()
return out
}
func main() {
threads, _ := strconv.Atoi(os.Args[1])
input, _ := io.ReadAll(os.Stdin)
t := table()
out := bufio.NewWriterSize(os.Stdout, 1<<20)
defer out.Flush()
var header []byte
var seq []byte
flush := func() {
if header != nil {
out.Write(header)
out.WriteByte('\n')
out.Write(rcSection(seq, &t, threads))
}
seq = seq[:0]
}
for _, line := range bytes.Split(input, []byte{'\n'}) {
if len(line) > 0 && line[0] == '>' {
flush()
header = line
} else {
seq = append(seq, line...)
}
}
flush()
}
revcomp.py
# reverse-complement (benchmarks game): read fasta on stdin, write each
# sequence reverse-complemented (IUB codes) in lines of 60.
# usage: revcomp.py < fasta > output
import sys
TABLE = bytes.maketrans(b"ACGTUMRWSYKVHDBNacgtumrwsykvhdbn", b"TGCAAKYWSRMBDHVNTGCAAKYWSRMBDHVN")
def emit(header, parts, out):
if header is None:
return
seq = b"".join(parts).translate(TABLE)[::-1]
out.write(header + b"\n")
out.write(b"\n".join(seq[i:i + 60] for i in range(0, len(seq), 60)) + b"\n")
out = sys.stdout.buffer
header, parts = None, []
for line in sys.stdin.buffer.read().split(b"\n"):
if line.startswith(b">"):
emit(header, parts, out)
header, parts = line, []
else:
parts.append(line)
emit(header, parts, out)
Fasta, the input of k-nucleotide and reverse complement
fasta.go
package main
// fasta (benchmarks game). The LCG is jumped ahead per chunk of lines, so the
// chunks are independent and the output is the same for any thread count.
// usage: fasta N THREADS (THREADS = 1: single-threaded) > output
import (
"bufio"
"os"
"strconv"
"sync"
"sync/atomic"
)
const im, ia, ic = 139968, 3877, 29573
const width = 60
const chunkLines = 1024
const alu = "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCCAGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"
type sym struct {
c byte
p float64
}
var iub = []sym{{'a', 0.27}, {'c', 0.12}, {'g', 0.12}, {'t', 0.27}, {'B', 0.02}, {'D', 0.02}, {'H', 0.02}, {'K', 0.02}, {'M', 0.02}, {'N', 0.02}, {'R', 0.02}, {'S', 0.02}, {'V', 0.02}, {'W', 0.02}, {'Y', 0.02}}
var homo = []sym{{'a', 0.3029549426680}, {'c', 0.1979883004921}, {'g', 0.1975473066391}, {'t', 0.3015094502008}}
func jump(s, k uint64) uint64 {
ra, rc := uint64(1), uint64(0)
a, c := uint64(ia), uint64(ic)
for ; k > 0; k >>= 1 {
if k&1 == 1 {
rc = (ra*c + rc) % im
ra = ra * a % im
}
c = (a*c + c) % im
a = a * a % im
}
return (ra*s + rc) % im
}
func runSection(header string, total, workers int, fill func(c int, buf []byte)) []byte {
lines := (total + width - 1) / width
out := make([]byte, total+lines)
step := chunkLines * (width + 1)
nchunks := (len(out) + step - 1) / step
var next int64 = -1
var wg sync.WaitGroup
for w := 0; w < workers; w++ {
wg.Add(1)
go func() {
defer wg.Done()
for {
c := int(atomic.AddInt64(&next, 1))
if c >= nchunks {
return
}
fill(c, out[c*step:min((c+1)*step, len(out))])
}
}()
}
wg.Wait()
return append([]byte(header+"\n"), out...)
}
func aluFill(c int, buf []byte) {
pos := c * chunkLines * width
for o := 0; o < len(buf); {
w := min(width, len(buf)-o-1)
for k := 0; k < w; k++ {
buf[o+k] = alu[(pos+k)%len(alu)]
}
buf[o+w] = '\n'
o += w + 1
pos += w
}
}
func randFill(seed uint64, table []sym, c int, buf []byte) {
cum := make([]float64, len(table))
acc := 0.0
for i, s := range table {
acc += s.p
cum[i] = acc
}
last := len(table) - 1
st := jump(seed, uint64(c*chunkLines*width))
for o := 0; o < len(buf); {
w := min(width, len(buf)-o-1)
for k := 0; k < w; k++ {
st = (st*ia + ic) % im
r := float64(st) / float64(im)
i := 0
for i < last && r >= cum[i] {
i++
}
buf[o+k] = table[i].c
}
buf[o+w] = '\n'
o += w + 1
}
}
func main() {
n, _ := strconv.Atoi(os.Args[1])
workers, _ := strconv.Atoi(os.Args[2])
seed2 := uint64(42)
seed3 := jump(seed2, uint64(3*n))
out := bufio.NewWriterSize(os.Stdout, 1<<20)
defer out.Flush()
out.Write(runSection(">ONE Homo sapiens alu", 2*n, workers, aluFill))
out.Write(runSection(">TWO IUB ambiguity codes", 3*n, workers, func(c int, b []byte) { randFill(seed2, iub, c, b) }))
out.Write(runSection(">THREE Homo sapiens frequency", 5*n, workers, func(c int, b []byte) { randFill(seed3, homo, c, b) }))
}
Pidigits
pidigits.co
--- pidigits (benchmarks game): Gibbons' unbounded spigot. An int is a whole
--- number of any size (spec, Numbers That Keep Every Digit), so the spigot
--- is its arithmetic, as Python's is: the compiled program holds its growing
--- numbers in the runtime's big ints and its counters in 64 bits.
---
--- With u = r + 3q (never negative) the spigot needs no signed numbers:
--- emit m: u <- 10(u - m t), m <- floor(u / t) of the new u, q <- 10q
--- else: m <- floor((q(k-1) + u l) / (t l)), u <- (u - q) l + 3qk,
--- q <- qk, t <- tl, l <- l+2, k <- k+1
--- and it emits m when q + u < (m+1) t. The sequential algorithm has no
--- loop to run on several cores.
n = 5000
q = 1
u = 3
t = 1
k = 1
l = 3
m = 3
count = 0
line = ''
while count lt n
if q + u lt (m + 1) * t
line = '`line``m`'
count = count + 1
if count % 10 eq 0
log['`line`\t:`count`']
line = ''
u = 10 * (u - m * t)
m = u / t
q = 10 * q
else
m = (q * (k - 1) + u * l) / (t * l)
u = (u - q) * l + 3 * q * k
q = q * k
t = t * l
l = l + 2
k = k + 1
if line ne ''
log['`fmt[line, width = 10, align = \'left\']`\t:`n`']
pidigits.rs
// pidigits (benchmarks game): Gibbons' unbounded spigot with num-bigint
// (pure Rust; the site's fast Rust programs use GMP).
// usage: pidigits N
use num_bigint::BigInt;
use num_traits::{One, Zero};
fn main() {
let n: usize = std::env::args().nth(1).unwrap().parse().unwrap();
let (mut q, mut r, mut t) = (BigInt::one(), BigInt::zero(), BigInt::one());
let (mut k, mut m, mut x) = (BigInt::one(), BigInt::from(3), BigInt::from(3));
let mut line = String::new();
let mut out = 0;
while out < n {
if (&q << 2) + &r - &t < &m * &t {
line.push_str(&m.to_string());
out += 1;
if out % 10 == 0 {
println!("{}\t:{}", line, out);
line.clear();
}
let nm = (BigInt::from(10) * (&q * 3 + &r)) / &t - &m * 10;
r = (&r - &m * &t) * 10;
q *= 10;
m = nm;
} else {
let tx = &t * &x;
let nm = (&q * (&k * 7 + 2) + &r * &x) / &tx;
r = (&q * 2 + &r) * &x;
q *= &k;
t = tx;
k += 1;
m = nm;
x += 2;
}
}
if !line.is_empty() {
println!("{:<10}\t:{}", line, n);
}
}
pidigits-Cargo.toml
[package]
name = "pidigits"
version = "0.1.0"
edition = "2021"
[dependencies]
num-bigint = "0.4"
num-traits = "0.2"
num-integer = "0.1"
[profile.release]
opt-level = 3
lto = "fat"
codegen-units = 1
pidigits.go
package main
// pidigits (benchmarks game): Gibbons' unbounded spigot with math/big.
// usage: pidigits N
import (
"fmt"
"math/big"
"os"
"strconv"
"strings"
)
func main() {
n, _ := strconv.Atoi(os.Args[1])
q, r, t := big.NewInt(1), big.NewInt(0), big.NewInt(1)
k, m, x := big.NewInt(1), big.NewInt(3), big.NewInt(3)
var tmp1, tmp2, tmp3 big.Int
ten, two, seven := big.NewInt(10), big.NewInt(2), big.NewInt(7)
var line strings.Builder
out := 0
for out < n {
// 4q + r - t < m t
tmp1.Lsh(q, 2)
tmp1.Add(&tmp1, r)
tmp1.Sub(&tmp1, t)
tmp2.Mul(m, t)
if tmp1.Cmp(&tmp2) < 0 {
line.WriteString(m.String())
out++
if out%10 == 0 {
fmt.Printf("%s\t:%d\n", line.String(), out)
line.Reset()
}
// q, r, m = 10q, 10(r - m t), 10(3q + r)/t - 10m
tmp1.Mul(q, big.NewInt(3))
tmp1.Add(&tmp1, r)
tmp1.Mul(&tmp1, ten)
tmp1.Quo(&tmp1, t)
tmp3.Mul(m, ten)
tmp1.Sub(&tmp1, &tmp3)
r.Sub(r, &tmp2)
r.Mul(r, ten)
q.Mul(q, ten)
m.Set(&tmp1)
} else {
// q, r, t, k, m, x = qk, (2q + r)x, tx, k+1, (q(7k+2) + rx)/(tx), x+2
tmp1.Mul(k, seven)
tmp1.Add(&tmp1, two)
tmp1.Mul(&tmp1, q)
tmp2.Mul(r, x)
tmp1.Add(&tmp1, &tmp2)
tmp3.Mul(t, x)
tmp1.Quo(&tmp1, &tmp3)
r.Add(r, new(big.Int).Lsh(q, 1))
r.Mul(r, x)
q.Mul(q, k)
t.Set(&tmp3)
k.Add(k, big.NewInt(1))
m.Set(&tmp1)
x.Add(x, two)
}
}
if line.Len() > 0 {
fmt.Printf("%-10s\t:%d\n", line.String(), n)
}
}
pidigits.py
# pidigits (benchmarks game): Gibbons' unbounded spigot with Python's ints.
# usage: pidigits.py N
import sys
def digits(n):
q, r, t, k, m, x = 1, 0, 1, 1, 3, 3
out = 0
while out < n:
if 4 * q + r - t < m * t:
yield m
out += 1
q, r, m = 10 * q, 10 * (r - m * t), (10 * (3 * q + r)) // t - 10 * m
else:
q, r, t, k, m, x = q * k, (2 * q + r) * x, t * x, k + 1, (q * (7 * k + 2) + r * x) // (t * x), x + 2
n = int(sys.argv[1])
line = ""
for i, d in enumerate(digits(n), 1):
line += str(d)
if i % 10 == 0:
print("%s\t:%d" % (line, i))
line = ""
if line:
print("%s\t:%d" % (line.ljust(10), n))
Binary trees
binarytrees.co
--- binary-trees (benchmarks game): build and check many trees of growing
--- depth. At each depth the trees are independent, so they are a `for`: they
--- run on every core, and `cpu[1]` makes it single-threaded.
--- a leaf is told by its own flag: a child is a Tree or 0, and a Tree is not
--- compared with 0
struct Tree [left, right, leaf]
fn make_tree[depth]
if depth eq 0
return Tree[0, 0, true]
return Tree[make_tree[depth - 1], make_tree[depth - 1], false]
fn check_tree[tree]
if tree.leaf
return 1
return 1 + check_tree[tree.left] + check_tree[tree.right]
--- the check of one tree of this depth
fn tree_check[depth]
return check_tree[make_tree[depth]]
n = 18
min_depth = 4
max_depth = max[n, min_depth + 2]
stretch = max_depth + 1
log['stretch tree of depth `stretch`\t check: `tree_check[stretch]`']
long_lived = make_tree[max_depth]
depth = min_depth
while depth lte max_depth
iterations = 2 ^ (max_depth - depth + min_depth)
checks = for j in range[iterations]
tree_check[depth]
total = 0
for c in checks
total = total + c
log['`iterations`\t trees of depth `depth`\t check: `total`']
depth = depth + 2
log['long lived tree of depth `max_depth`\t check: `check_tree[long_lived]`']
binarytrees.rs
// binary-trees (benchmarks game). At each depth the trees are built and
// checked THREADS at a time (1 = single-threaded).
// usage: binarytrees N THREADS
use std::sync::atomic::{AtomicU64, Ordering};
struct Node(Option<Box<Node>>, Option<Box<Node>>);
fn make(d: u32) -> Node {
if d == 0 {
Node(None, None)
} else {
Node(Some(Box::new(make(d - 1))), Some(Box::new(make(d - 1))))
}
}
fn check(n: &Node) -> u64 {
1 + n.0.as_ref().map_or(0, |c| check(c)) + n.1.as_ref().map_or(0, |c| check(c))
}
fn main() {
let a: Vec<usize> = std::env::args().skip(1).map(|s| s.parse().unwrap()).collect();
let (n, threads) = (a[0] as u32, a[1]);
let min_depth = 4u32;
let max_depth = n.max(min_depth + 2);
let stretch = max_depth + 1;
println!("stretch tree of depth {}\t check: {}", stretch, check(&make(stretch)));
let long_lived = make(max_depth);
let mut d = min_depth;
while d <= max_depth {
let iterations = 1u64 << (max_depth - d + min_depth);
let next = AtomicU64::new(0);
let total = AtomicU64::new(0);
std::thread::scope(|s| {
for _ in 0..threads {
s.spawn(|| {
let mut sum = 0;
// take iterations in blocks to keep the counter cold
loop {
let start = next.fetch_add(64, Ordering::Relaxed);
if start >= iterations {
break;
}
for _ in start..(start + 64).min(iterations) {
sum += check(&make(d));
}
}
total.fetch_add(sum, Ordering::Relaxed);
});
}
});
println!("{}\t trees of depth {}\t check: {}", iterations, d, total.load(Ordering::Relaxed));
d += 2;
}
println!("long lived tree of depth {}\t check: {}", max_depth, check(&long_lived));
}
binarytrees-mimalloc-Cargo.toml
[package]
name = "bt_mi"
version = "0.1.0"
edition = "2021"
[dependencies]
mimalloc = "0.1"
[profile.release]
opt-level = 3
lto = "fat"
codegen-units = 1
binarytrees.go
package main
// binary-trees (benchmarks game). At each depth the trees are built and
// checked THREADS at a time (1 = single-threaded).
// usage: binarytrees N THREADS
import (
"fmt"
"os"
"strconv"
"sync"
"sync/atomic"
)
type node struct{ l, r *node }
func make_(d int) *node {
if d == 0 {
return &node{}
}
return &node{make_(d - 1), make_(d - 1)}
}
func check(n *node) int64 {
if n.l == nil {
return 1
}
return 1 + check(n.l) + check(n.r)
}
func main() {
n, _ := strconv.Atoi(os.Args[1])
threads, _ := strconv.Atoi(os.Args[2])
minDepth := 4
maxDepth := max(n, minDepth+2)
stretch := maxDepth + 1
fmt.Printf("stretch tree of depth %d\t check: %d\n", stretch, check(make_(stretch)))
longLived := make_(maxDepth)
for d := minDepth; d <= maxDepth; d += 2 {
iterations := int64(1) << uint(maxDepth-d+minDepth)
var next, total int64
var wg sync.WaitGroup
for t := 0; t < threads; t++ {
wg.Add(1)
go func() {
defer wg.Done()
var sum int64
for {
start := atomic.AddInt64(&next, 64) - 64
if start >= iterations {
break
}
for i := start; i < min(start+64, iterations); i++ {
sum += check(make_(d))
}
}
atomic.AddInt64(&total, sum)
}()
}
wg.Wait()
fmt.Printf("%d\t trees of depth %d\t check: %d\n", iterations, d, total)
}
fmt.Printf("long lived tree of depth %d\t check: %d\n", maxDepth, check(longLived))
}
binarytrees.py
# binary-trees (benchmarks game). At each depth the trees are built and
# checked PROCS at a time (1 = single-process).
# usage: binarytrees.py N PROCS
import sys
from multiprocessing import Pool
def make(d):
if d == 0:
return (None, None)
return (make(d - 1), make(d - 1))
def check(t):
l, r = t
if l is None:
return 1
return 1 + check(l) + check(r)
def work(args):
d, count = args
return sum(check(make(d)) for _ in range(count))
if __name__ == "__main__":
n, procs = int(sys.argv[1]), int(sys.argv[2])
min_depth = 4
max_depth = max(n, min_depth + 2)
stretch = max_depth + 1
print("stretch tree of depth %d\t check: %d" % (stretch, check(make(stretch))))
long_lived = make(max_depth)
pool = Pool(procs) if procs > 1 else None
for d in range(min_depth, max_depth + 1, 2):
iterations = 1 << (max_depth - d + min_depth)
if pool:
blocks = max(procs * 4, 1)
sizes = [iterations // blocks + (1 if i < iterations % blocks else 0) for i in range(blocks)]
total = sum(pool.map(work, [(d, c) for c in sizes if c]))
else:
total = work((d, iterations))
print("%d\t trees of depth %d\t check: %d" % (iterations, d, total))
print("long lived tree of depth %d\t check: %d" % (max_depth, check(long_lived)))
