|
Array Math Functions, Index Lists & Slices, Solve(), Sort(), and More |
Array Math and Solve()
This page picks up where the Arrays page (Arrays) leaves
off. There, an array is a list you fill, grow and print; here it is
something you compute with as a whole — scale every
element in one statement, multiply matrices, and solve systems of
equations. Read the Arrays page first; everything below
assumes dim, fill, expandable arrays and
print of an array. Two companion pages continue from here:
outer() pairs every element of one array with every element
of another (outer(): Every Element of One Array Against Every Element of Another), and reduce() /
scan() fold an array down or run it up
(reduce() and scan(): Folding an Array Down to One Value, or Running It Up).
Each section below opens with the idea; click a Show me line for the code and the details, or Expand all.
1. Whole-Array Arithmetic
You can assign a whole array expression to an expandable array. Every element is computed at once, and the array takes the shape of the arrays on the right.
For an arithmetic expression, the target is a REAL array. An
expandable one (dim x(*)) is reshaped to fit
whatever the right side produces, so you never size it yourself. A
fixed array (dim x(5), dim g(2, 3))
works too — x = x * 2 is fine — but it is not
resized: its shape must already match the result (dimensions, bounds and
origins), or the catchable exception ARRAYSHAPE is raised.
So an expandable target adapts, a fixed target checks. String arrays take
only the move forms — copy, reshape(),
transpose(), filter() — shown in section 5.
On the right side, any mix of arrays, numbers and + - * / ^
works, including unary minus and parentheses. Every array in one
expression must have the same shape; if not, the catchable exception
ARRAYSHAPE is raised. A scalar part — the 1.12
above, or an expression such as 1 + tax_rate — is
computed once, not once per element.
You do not need to assign the result before printing it or passing
it to a function. print takes a whole-array expression
directly — print prices * 1.12,
print due + 5 — and prints the result exactly as it
prints an array. The same goes wherever a function expects an array:
stats$sum(prices * 1.12) and size(prices + 1)
work directly. The statistics functions do not care about shape,
either: stats$sum(x) on a 2 by 3 array adds all six cells.
2. Masks: Comparing a Whole Array
Six students sat an exam and the pass mark is 60. Two questions come up at once: how many passed, and what were the passing scores? Comparing the whole array against the mark answers both, element by element: 1 where the test holds, 0 where it does not. That mask of ones and zeros is something you can count, multiply by, and combine:
Show me: comparison masks in code
Summing the mask answers the first question. Multiplying by it keeps
the passing scores and zeroes the rest (section 4 has a cleaner way to
pull them out). All six comparisons work this way (= <> <
<= >= >), against a number or an array of the same
shape. and, or and not combine
masks.
The comparisons are exact. Sheerpower arithmetic is
exact decimal, so a comparison on decimal values comes out the way it
does on paper. In a floating-point language 0.1 * 3 is
0.30000000000000004 and the test below is false:
Counting with a mask is the idiom, and it is cheap.
"How many elements meet this condition?" is
stats$sum(array condition) — the comparison gives a
mask of 1s and 0s, and summing it counts the 1s. No loop, no counter
variable, one line that reads as the question. Throw a thousand dice
and ask how many came up over 3:
The same line counts ten million as easily, and beats the loop that would replace it:
It is cheap because a mask is not an array of numbers: every boolean
array — a ? name, and every comparison result —
is kept one bit a flag. A hundred million flags take
12.5 MB, a billion 125 MB, where an array of REALs would take
3.2 GB and 32 GB. Counting, and, or
and not run on the packed words, 64 flags a step: a billion
flags are counted in 0.11 s and combined with and in
0.3 s. A boolean element can only hold 0 or 1 — storing
anything else raises NOT0OR1 at the store — and
typeof$(flags?) shows *Packed*.
3. Shapes: reshape() and redim
Sometimes the numbers are right but the shape is wrong — you have a
flat list and want a grid, a 3 by 3 identity matrix, a repeating stripe, or
you want to flatten a grid back to a plain vector.
reshape(a, dims...) pours an array's elements, in storage
order, into any new shape; redim resizes an array you already
have (keeping the data); and reshape(a) with no sizes flattens
it back to one dimension.
Show me: reshape, recycle, redim and flattening
You need a 3 by 3 identity matrix, or a row of alternating ones and
zeros, or a small vector poured into a grid.
reshape(a, dims...) pours the elements of a,
in storage order, into a new shape. By default the first elements fill
it and any remaining cells are zero. With the option
recycle: true the source repeats until the shape is full:
The identity matrix works because the pattern is one element longer than a row. Recycled to nine elements and poured into rows of three, the 1 lands one place further along in each row:
Any size follows the same rule: a 1 and four zeros recycled into 4 by 4 is the 4 by 4 identity.
recycle: true is written the way every optional setting of
a built-in function is written in Sheerpower — by name, after the
positional arguments — so a size can never be mistaken for it.
Two tools now repeat things, and they read best in different places.
n of (...) in a fill repeats a literal
list you are typing; recycle: true repeats the contents of
an array you already have, into any shape, including a grid:
reshape() builds a new array. redim changes
the shape of an expandable array you already have — new bounds,
or a new number of dimensions — and keeps the data in storage
order (row by row). Array assignment does the same thing on its own.
The target takes the shape of the result: a two-dimensional operand
makes x two-dimensional, and a later one-dimensional
result makes it a plain expandable array again.
Flattening. To pour any array back into a plain vector,
call reshape with no sizes at all: reshape(a) is
every element of a, in storage order, as a one-dimensional
array — APL's ravel. It is the readable form of
reshape(a, size(a)) and works on numeric and string arrays
alike; it composes like any array expression
(stats$sum(reshape(grid)), print reshape(grid)).
For an expandable array you already hold, redim x(*) does the
same thing in place.
The shape can change while the program runs, so a reference with the
wrong number of subscripts — x(1) while
x is two-dimensional — is caught at runtime as the
catchable exception WRONGNUMDIMS. An array declared with a fixed shape
(dim grid(2, 3)) keeps its number of dimensions
for good: redim grid(4, 5) changes the bounds, but
redim grid(6) is a compile error, and so is
grid(6) — its subscript count is checked when the
program compiles.
4. Selecting and Replacing with a Mask: filter() and where()
A comparison on an array gives a mask of ones and zeros.
filter(values, mask) pairs the two arrays element by
element in storage order (row by row) and returns the values whose
mask element is 1. The result is always one-dimensional. The two
arrays need not have the same shape, but they must hold the same
number of elements — a mask built from the wrong array is a
mistake far more often than a plan, so a mismatch raises
ARRAYSHAPE naming both counts. When a short or long mask
is what you mean, say so: filter(values, mask, partial:
true) selects nothing past a short mask's end and ignores a long
mask's extra elements.
If you come from NumPy, this is the one rule to unlearn:
here parentheses and brackets mean different things
when what is inside them is an array. v(m) with an array
m is a MASK — ones and zeros choosing elements
(section 13 stores through it) — while v[y] with an
array y is an INDEX LIST, the elements at positions
y(1), y(2), ... in that order (section 11). In NumPy parentheses never
index at all, so v(m) reads as a call there; here it is the
mask, and a 2 in it is the loud NOT0OR1 exception rather than
a silent third element.
Show me: filter(), where() and picking one list by another
Picking one list by another. A fairground ride requires riders to be at least 120 cm tall. Hand it the riders' names and their heights, and a single comparison on the heights chooses names out of the names — that is filter's signature move, selecting from one list by a condition on a parallel one. The same masks feed the statistics, so the average height of those turned away is one more line:
filter() drops the elements you do not want, so its
result is shorter and always one-dimensional. Its companion,
where(cond, a, b), keeps the shape: it walks the mask
and, element by element, takes from a where the mask is 1 and
from b where it is 0 — the vectorized "if this, use that,
otherwise the other." Either value may be a single number or string (used
wherever it is chosen) or a whole array the same size as the mask.
Show me: where() — clamp, label, and choose between two lists
The mask (cond) sets the result's shape, so a 2-D mask gives a
2-D result. Any array branch must hold the same number of elements as the
mask (a mismatch raises ARRAYSHAPE); a single value broadcasts
to every element it is chosen for, and a and b
must agree in kind — both numbers or both strings. Where
filter answers "which elements pass?", where
answers "what should each element become?"
5. String Arrays Too
That last line selected from a string array with a numeric mask.
Moving elements needs no arithmetic, so the whole-array copy,
reshape(), transpose() and
filter() all work on string arrays, and they nest:
Show me: reshape, transpose, copy and filter on string arrays
Those are the whole-array MOVES a string array supports —
copy, reshape(), transpose() and
filter(). The one string OPERATOR is +,
concatenation element by element (section 20); every string
FUNCTION maps over a string array (section 10), brackets pick by
position (section 11) and split() / join$()
turn text into an array and back (section 14). A reshaped string array may be printed, measured
with size(), or handed to a routine that takes an array
parameter (by name, or through the pipe operator |> from
the Operators page, Mathematical and Logical Operators).
6. Matrix Multiply: matmul()
matmul(a, b) multiplies two matrices — rows from
a, columns from b, each cell an exact sum of
products. It is the building block of the linear algebra on this page: a
weighted total, or the normal equations behind a best-fit line. A plain
vector can stand in for a row or a column; two vectors are ambiguous, so
Sheerpower asks you to shape one.
Show me: the code, the vector rule, and summing along an axis
matmul(a, b) multiplies two matrices. With a
2 by 3 and b 3 by 2, z = matmul(a, b) makes
z 2 by 2: rows from a, columns from
b, each element an exact sum of products. Both operands must
be two-dimensional and the inner sizes must agree; if not, the catchable
exception ARRAYSHAPE names the two shapes. It is a function, so it nests
and mixes with the other operators like any array expression:
z = matmul(a, b) * 2.
A plain vector may stand in for a matrix. On the right of
matmul() it acts as a column, so a matrix times a vector
gives a vector. On the left it acts as a row. A vector of ones makes the
two cases easy to see — one gives the row sums of a,
the other its column sums:
The matmul(..., ones) form is a nice way to see how a vector
acts in matrix multiply, but it is not how you would normally sum along an axis.
For that, reach for reduce(): reduce(a, +) gives
the row sums (6, 15) directly, and
reduce(a, +, axis: 1) the column sums (5, 7, 9).
See reduce() and scan() for folding along any axis of any array.
No ones vector, no matmul() — axis: n
folds dimension n, so axis: 2 reduces across each
row and axis: 1 down each column. This is the idiomatic way to
sum (or max(), min(), ...) along an axis of any
array, at any rank.
The dot product of two vectors — the sum of their paired
products — is dot(a, b) (section 18). For
a = 1, 2, 3 and b = 4, 5, 6 that is
32; stats$sum(a * b) is the same number the long
way round. (matmul() is for matrices; two vectors are a job
for dot().)
Like every array function, matmul() works without an
assignment: print matmul(a, b) prints the product, and
stats$sum(matmul(a, b)) totals it.
7. Rows into Columns: transpose()
Data often arrives one way and is wanted the other. A shop sells
coffee beans and tea, and its sales come in as one row per month with
a column per product: in month 1 it sold 120 bags of beans and 80 of
tea, in month 2 150 and 95, in month 3 170 and 110. A report wants one
row per product with the months across. transpose() turns
the 3 by 2 into a 2 by 3:
Show me: transpose() in code
A vector transposes to a 1 by n matrix. transpose() works
on string arrays too (section 5), and it composes with the other
array functions — matmul(transpose(a), a) is the heart of
the last problem on this page.
8. Solving Equations: solve()
A surprising number of everyday questions have the shape "I know
the totals — what were the parts?" Each total is one
equation; the parts are the unknowns. When there are as many
independent totals as unknowns, solve() finds the parts
in one statement. Here are three such problems, and what each one
teaches.
solve() is exact up to 64 unknowns; past that it asks you to
choose, with one word, between the exact answer and double precision. The
note on size at the end of this section gives the measured figures and the
exact: option.
Show me: three worked problems (cafe prices, a materials mix, least squares)
Problem 1: what does the cafe charge?
A cafe sells only coffee and muffins, and the till records totals, not prices. On Monday it sold 30 coffees and 20 muffins for 130.00; on Tuesday, 25 coffees and 30 muffins for 145.00. What are the two prices?
Why solve(): two unknowns (the prices), two totals
(the days). Written out, Monday is 30c + 20m = 130 and
Tuesday is 25c + 30m = 145. Put the coefficients in a matrix
(one row per day, one column per unknown) and the totals in a vector,
and solve(sales, takings) returns the unknowns in column
order:
Coffee is 2.50 and a muffin 2.75 — printed exactly, not as 2.4999999999999999.
Why the answer is exact. solve() uses
fraction-free elimination in Sheerpower's exact decimal
arithmetic. The only place an answer can round is one final division
per unknown, so an answer that terminates comes out exact.
Why exact arithmetic, and not floating point. Most languages solve in IEEE double precision — about sixteen significant digits, rounded at every step. Sheerpower does the elimination in exact decimal and rounds only at the final division. That changes what limits the answer:
- Floating point loses digits to ill-conditioning. A double spends roughly as many of its sixteen digits as the system's condition number has digits — a measure of how close a system's facts are to saying the same thing, pushed up by near-duplicate rows (common in real pricing and allocation data). A well-behaved system loses two or three digits; one whose condition number nears 1016 can be left with none, and returns a wrong answer with no warning.
- Sheerpower does not. Rounding only once, the condition
number never eats into the answer. In testing,
solve()was exact for a dense system and still exact at a condition number of 10200. - The only real limit is the size of the numbers, not the conditioning — 54 whole digits and 16 decimal places kept exactly, 64 significant digits beyond that. For the sizes on this page that limit never comes into view.
A note on size. Because every step stays exact,
solve() does not lose accuracy as a system grows — what
grows instead is the time it takes:
- Accuracy holds. In testing it returned the exact answer for systems of several dozen unknowns built from four- and eight-digit numbers, and even for answers as large as eighteen digits.
- Speed is what you meet first, so past 64 unknowns
solve()makes you choose. The exact kernel's time multiplies roughly ten- to eighteen-fold each time the number of unknowns doubles — the whole-number intermediates grow as well as the operation count — and at a thousand unknowns it takes about two minutes. So up to 64 unknowns a call is exact, as above. Beyond that a plainsolve(a, b)does not quietly change what it computes: it raisesEXACTCHOICE, a catchable exception whose detail names the size and both options, and you add one of them —exact: truefor the exact answer and the wait, orexact: falsefor double precision: the matrix is copied to IEEE doubles, factored there (an LU with partial pivoting, as a numerical library would), corrected by one refinement step whose residual is computed in exact arithmetic, and the answers come back as each double's value to 17 significant digits. Onesolve()call, exact up to 64 andexact: falsebeyond, measured on a low-end laptop against a system whose true answer is known:
(Each time is one call, averaged over repeated calls.) The double kernel's time grows about with the cube of the unknowns, from a far lower floor: a thousand unknowns in about a sixth of a second, where the exact kernel takes about two minutes, and the answers within a few parts in 1015.unknowns kernel time largest error 16 exact 0.26 ms 0 32 exact 1.0 ms 0 48 exact 1.8 ms 0 64 exact 3.8 ms 0 80 double 1.8 ms 6e-14 100 double 1.9 ms 9e-16 200 double 5.6 ms 1.5e-15 500 double 40 ms 1.7e-15 1000 double 0.16 s 1.4e-15 - Why 64. The threshold is not about accuracy: the exact
kernel is exact at every size — past about fifteen unknowns it solves
the system modulo many primes at once, spread across the processor cores,
and rebuilds the exact answer from them (in testing, every unknown of
systems up to a thousand came back exact). What it decides is the size
past which a call must say which it wants, and the exact kernel's cost
around and past it is:
Two properties make a count of unknowns the right kind of line, whatever the number: it is deterministic — the same program behaves the same on every machine, where a line drawn on elapsed time would not — and it is data-independent: it depends only on the size of the system, never on its values. And the reason the choice is yours rather than automatic is the one honest cost of double precision: a price that should be 7.05 can come back as 7.0499999999999998, because the double kernel is accurate to about one part in 1015, not exact — a change of meaning a program should state, never receive in silence.unknowns exact kernel 48 1.7 ms 64 3.7 ms 100 16 ms 128 36 ms 200 0.18 s 300 0.88 s 500 6.5 s 1000 about 2 minutes - Choosing for yourself:
exact:.solve(a, b, exact: true)runs the exact elimination at any size — when the answer must be exact and you can wait — andsolve(a, b, exact: false)runs the double kernel at any size. The cafe again, both ways, and then a 100-unknown system whose true answer is known (the coefficients and the answer are two-decimal numbers):print solve(sales, takings) ! 2.5 2.75 -- exact, as always at this size print solve(sales, takings, exact: false) ! 2.5 2.75 -- doubles + refinement: clean here too when exception in z = solve(a, b) ! 100 unknowns, no choice made: it raises use print extext$ ! Too big for the exact default -- choose ! exact: true or exact: false print _string$ ! solve(): 100 unknowns is past the 64 the exact elimination end when ! is the default for -- add exact: true (the exact answer, ! whose time climbs steeply as systems grow) or exact: false ! (double precision, much faster on large systems) start timer z = solve(a, b, exact: false) ! double precision print _elapsed, z(1) ! 0 7.0499999999999998 start timer z = solve(a, b, exact: true) ! the exact kernel print _elapsed, z(1) ! 0.016 7.05The refinement step usually lands a small system on the exact double (2.5 and 2.75 are representable), soexact: falseon a handful of unknowns is hard to tell from exact; at a hundred unknowns the difference is visible in the last digits, andexact: truebuys them back for about 18 ms against the double kernel's 2.5 ms (one call of either is within a tick of the timer's 15.6 ms resolution, which is what the example prints).EXACTCHOICEisexceptiontype('exactchoice'), and like every runtime exception it leaves its explanation in_string$. In double precision aNaNorInfelement is refused (NUM_OUTOFRANGE, naming the element), and a matrix whose elimination leaves a pivot within about n × 10-13 of its largest entry is reportedSINGULAR— two rows nearly proportional, which the exact kernel would still separate. - The same rule for
determinant()andlstsq(). Both share the exact elimination's cost curve, so both take the 64-unknown line and theexact:option:determinant(a, exact: false)factors in doubles and multiplies the diagonal (a 1000 by 1000 determinant, near 101743, in 0.16 s — carried as a logarithm, since no double holds it, and returned as a scientific REAL), and a matrix it finds singular gives 0, so at every sizedeterminant(a)is 0 exactly whensolve(a, b)raisesSINGULAR.lstsq(a, b, exact: false)forms the normal equations in doubles too and refines with the exact residual: a 2000 by 200 fit in 74 ms (its exact form takes 0.61 s), 5000 by 500 in 0.44 s, both within 10-15 of the true coefficients. inverse()too.inverse(a)issolve(a, identity), so it takes the same 64 line and the same option.inverse(a, exact: false)factors in doubles once, solves every column of the identity, and refines them with the exact residual, which runs on the same 256-bit kernel asmatmul()across the processor cores: 500 by 500 in 0.5 s, 1000 by 1000 in 3 s, withmatmul(a, inverse(a))within 4 × 10-14 of the identity. Anysolve(a, b, exact: false)whose b has many columns gains the same way: a thousand right-hand sides at a thousand unknowns took 26 s before the residual moved to that kernel.matmul()stays exact, and got faster instead. A product of two-decimal matrices is exact by nature, so it was given no double kernel: when every element is an ordinary real, the row-and-column sums are accumulated exactly in 256-bit integers and rounded once at the end, the rows spread across the processor cores — 200 by 200 in 15 ms, 500 by 500 in 0.22 s, 1000 by 1000 in 1.7 s, every cell exact. One consequence worth knowing: because the sum is rounded once rather than product by product, a matrix whose products carry more than sixteen decimals can differ from the old answer by one unit in the sixteenth place, in the more exact direction; money data is unaffected.- In practice. Typical business systems run from a
handful to a few dozen unknowns, solved anywhere from a few times a day to
once per order or request — all of it exact and effectively instant.
Larger systems are a fraction of a second with
exact: false, and withexact: trueup to about three hundred unknowns; what still belongs to a numerical library is the territory the next section names — eigenvalues, factorizations, sparse matrices — not size. - The one place an answer still rounds is a value that never terminates — like one-third — which comes back to sixteen decimal places, exactly as any division would.
Problem 2: how many of each vehicle?
A delivery company knows it runs 8 vehicles, that together they carry 21 tons, and that the fleet costs 1650 a day. A van carries 1 ton and costs 100 a day; a truck 2 tons and 150; a lorry 5 tons and 400. How many vans, trucks and lorries are there?
Why solve(): three unknowns, three independent facts — a count, a capacity, a cost. Each fact is a row: what one van, one truck and one lorry contribute to it.
One van, five trucks, two lorries. The same shape serves many "blend" questions — ingredients meeting nutritional targets, products meeting material and cost totals, staffing meeting hours and budget — whenever each fact is a plain sum of the unknowns.
One limit to know: solve() solves the equations and
nothing more. It does not know that vehicles come in whole numbers or
that a count cannot be negative. This data happens to give 1, 5 and 2,
but change the daily cost to 1675 and it answers 1.75 vans, 4 trucks
and 2.25 lorries without complaint.
When the answer has to be whole, use diophantine()
instead (Problem 5 below). It takes the same two arrays, gives the
same 1, 5, 2 here, and refuses a fractional answer loudly:
When the answer must also be non-negative, or when there are more
unknowns than facts and one plan is to be preferred over the others,
that is the maximize / minimize block with
whole unknowns and = rows
(Resource Allocation with Maximize and Minimize). solve() itself enforces
none of this; it reports what the equations say.
Problem 3: which way are sales heading?
Six months of sales: 13, 13, 15, 19, 20, 22. They rise, but not in a straight line — no line passes through all six points. The question is the best straight line, the one that misses the points by the least (in the least-squares sense), and what it predicts for month 7.
Why transpose(): this time there are more facts than
unknowns. Each month is one equation (slope × month +
intercept = sales), so there are six equations and only two
unknowns, and solve() needs a square system — as
many equations as unknowns. It cannot take this one as it stands.
The remedy is to multiply both sides by the transpose of the
coefficient matrix. The transpose is 2 by 6, and 2 by 6 times 6 by 2
is 2 by 2 — square — and the solution of that smaller
system is the best-fit line. (Mathematicians call the smaller system the normal
equations.) Build the coefficient matrix with one row per month
(month, 1), and the whole method is one line:
solve(matmul(transpose(a), a), matmul(transpose(a), y)). It fits a
line through any set of points, or a plane if you add columns.
That normal-equations spelling — solve(matmul(transpose(a), a),
matmul(transpose(a), y)) — is common enough to have a shorthand:
lstsq(a, y) is exactly it, the best-fit answer in one call.
Give it an overdetermined system (more rows than unknowns) with a
full-rank a and it returns the least-squares solution in the
same exact arithmetic; a square a reduces to
solve(). If the columns are dependent, or there are fewer
rows than unknowns, aTa is singular and it raises
SINGULAR, just as solve() would.
Sales grow by 2 a month from a base of 10, and month 7 should bring
24. (matmul(transpose(design), sales) multiplies a matrix by
the plain sales vector — the vector-as-column rule
from the Matrix Multiply section.)
One thing to expect with real data: this data was chosen so the line comes out to exactly 2 and 10. Real best-fit lines almost never terminate — a slope such as 1.8571428571428571 is normal — and a long decimal is the arithmetic being honest, not wrong. The rule from Problem 1 still holds: an answer that terminates comes out exact, one that does not is rounded once, at the end.
A note for the numerically inclined. Forming
matmul(transpose(a), a) — the normal equations — is a
step floating-point libraries usually avoid, because it squares the
problem's condition number and can throw away precision. That warning
is about floating-point rounding, and it does not apply here:
Sheerpower's elimination is exact, so there is no accumulated rounding
to amplify. The normal-equations method is a sound, direct way to fit a
line or a plane in this arithmetic.
Problem 4: what is the best plan?
The first three problems had one right answer. The fourth kind of question has many workable answers and asks for the best one: a bakery makes croissants and muffins, earns 2.50 on a croissant and 2.75 on a muffin, and each morning has 100 kg of flour, 220 oven slots and 80 eggs. A croissant uses 0.5 kg of flour, one oven slot and 0.2 eggs; a muffin 0.4 kg, one slot and 0.5 eggs. How many of each should it bake?
Why maximize(): a profit per unit for each variable, a
limit per resource, and every quantity at least zero — that is a
linear program. Put the profits in a vector, the usage in a
matrix (one row per resource, one column per product) and the limits in
a vector, and maximize(profit, use, have) returns the plan
that earns the most while staying inside every limit:
Bake 100 croissants and 120 muffins for 580.00, and the two resources
that ran out — oven slots and eggs — are the ones worth
adding. minimize(cost, use, have) is the same call for a
cost to keep down. A "greater than" requirement is written with
negatives (x + y >= 4 is -x - y <= -4), and a plan that
must be met exactly as two rows, one each way.
Why the answer is exact. The method is the simplex,
which walks from corner to corner of the feasible region, and each step
is a Gauss-Jordan pivot — the same kind of step solve()
takes — so it runs in the same fraction-free exact arithmetic: no
feasibility tolerance, no "treat 10-9 as zero", and a
corner that terminates comes out exact, including fractional ones
(a plan of 1.25 and 1.25 prints as 1.25, not 1.2499999999999998).
When there is no plan. Two exceptions, both catchable
(see Exception Handling): INFEASIBLE when no quantities
satisfy every constraint (a "make at least 300" row against
220 oven slots), and UNBOUNDED when the profit can grow
without limit because a limit is missing — the detail says
"the objective can be made arbitrarily large -- a constraint is
missing".
Size. The same rule as solve(), with its
own line: past 128 variables plus constraints a plain call raises
EXACTCHOICE, and exact: true or
exact: false chooses. The line is higher than
solve()'s 64 because the exact simplex is about ten times
quicker than the exact elimination at the same size — 128 variables
plus constraints in a quarter of a second, 64 variables under 128
constraints in half a second, 100 under 200 in four — and
exact: false answers those in a tick. Those figures are
for money-like data. Because every row is first scaled to whole
numbers, data whose entries span many orders of magnitude (0.001 next
to 1000) makes the exact simplex several times slower at the same
size, while doubles do not care; there, exact: false is
the practical choice well below the line. In testing, the double
kernel landed on the same corner as the exact one every time, with the
plan good to about twelve digits at 750 variables plus constraints;
what it cannot promise is the corner itself on a degenerate problem,
where a tie decided against a tolerance can walk it one step short. Whole-unit
answers (you cannot bake 0.4 of a muffin) are the next problem.
Whole units, ranked goals, and the story problems are on
their own page: Resource Allocation with Maximize and Minimize — the woodshop in
whole units (whole: true), the block statement that writes a
plan as it is told, and the machine shop with two ranked goals. That page
compiles to the call above, so everything on this page about exactness and
size holds there too.
Problem 5: which whole numbers?
A last kind of question asks not for the best answer but for the whole-number ones. Stamps come in 3-cent and 5-cent values: how many of each make exactly 47 cents? A parcel line packs in boxes of 6 and 9: can it ship exactly 20? These are Diophantine equations — whole-number coefficients, whole-number unknowns — and they have a clean theory: a solution exists exactly when the greatest common divisor of the coefficients divides the total, and when one solution exists there is a whole family of them, spaced along a regular lattice.
Why diophantine(): diophantine(a, c) returns
one whole-number solution of a1x1 + ... +
anxn = c, and lattice(a) the
directions along which every other solution lies, one per row:
Every solution is how plus a whole multiple of each row of
trade: for two unknowns that is one direction, for
n unknowns n − 1 of them. The two-unknown
answer comes back in the textbook form, with the first unknown
reduced to the range 0 to a2 / g, which is why 47
cents gives 4 and 7 rather than one of its more distant relatives.
diophantine() raises the catchable NOSOLUTION
when the divisibility test fails, and its detail names the gcd and the
total, so the message is the explanation.
Systems. With a matrix of coefficients and a vector of totals the same two names solve several equations at once, through the Hermite normal form; the lattice then has as many rows as the system leaves free. The answer comes back small: the reduction alone lands on a solution some fifty digits out, so the lattice rows are shortened first (the LLL reduction) and the solution is walked back to the nearest lattice point — a system with a two-digit answer gets a two-digit answer, and the lattice rows are the short ones.
What this is not. A Diophantine solution can be
negative, as the second stamp solution above shows, and there is no
"fewest stamps" here. Non-negative or best whole-number
answers are the maximize / minimize blocks
(Resource Allocation with Maximize and Minimize) — there the totals are limits
and one plan is best; here the totals are exact and every answer
counts. What
diophantine() settles is whether a whole-number answer
exists at all and what every one of them looks like, exactly: the
coefficients and totals are whole numbers in the exact REAL (54 digits),
the work inside is done on 115-digit whole numbers, and it raises
NUMOVER rather than round if a step would need more. That
step is the back-substitution, which grows with the number of equations
and their digits: about eight equations of three-digit coefficients, or
five of six-digit ones, with any number of unknowns beyond that; a few
equations in many unknowns go much further. A square system
with one solution never meets that wall — it takes the exact
solve() path (64 unknowns, any digits) and asks only
whether that one solution is whole.
When there is no answer
If the facts are not independent — Tuesday was simply double
Monday — the system has no unique solution and
solve() raises the catchable exception
SINGULAR. Catch it where a bad data set is a real
possibility (exception handling has a page of its own:
Exception Handling):
Shape mistakes — a matrix that is not square, a totals vector
with the wrong number of rows — raise ARRAYSHAPE,
and the message names the shapes it saw.
The determinant: is my information enough?
Before you solve, a single number tells you whether an answer even exists. A cafe owner wants the price of a coffee, a tea and a muffin, and has three group receipts — each showing what was bought and the total paid:
- 2 coffees, 1 tea, 1 muffin — paid 10.75
- 1 coffee, 3 teas — paid 10.75
- 1 tea, 2 muffins — paid 8.75
Three unknowns and three receipts — but are the receipts
enough? Only if the three orders are genuinely different.
determinant() of the who-bought-what grid answers that in one
number: any value that is not zero means the orders are independent enough
to recover all three prices; a value of zero means they are not,
and no arithmetic can pull the prices apart.
Now suppose the third receipt were just "two of everything on the
first" — it repeats what you already knew and adds nothing new.
The grid's determinant collapses to zero, warning you before
you trust a result that these receipts cannot pin the prices down (it is the
same thing solve() would report as SINGULAR):
The determinant is computed by the same exact arithmetic as
solve(), so a grid of whole numbers or terminating decimals
gives an exact answer — and it follows the same size rule: past 64
unknowns it asks for exact: true or exact: false
(see the note on size above). Where solve() finds the
prices, determinant() tells you first whether they can be
found at all.
The inverse: work out a mix once, use it every week
A coffee roaster sells three blends, each a fixed mix of three beans: the
house blend is 40% Brazil, 30% Colombia and 30% Ethiopia, the espresso
10%, 70% and 20%, the light roast 30%, 10% and 60%. Every week new stock
arrives and the question is the same: how many kilos of each blend use
this stock up exactly? That is a solve() every week —
or one inverse(), worked out once, that turns any
week's stock into blends with a single matmul():
Every entry of recipe is exact: the mix's determinant is
0.1 (determinant(mix) says so), and dividing by it leaves
terminating decimals, so the blends come out as whole kilos with nothing
left over. A negative answer is the arithmetic saying that no amount of
the three recipes uses that stock up — there is too much Brazil for
them — so the roaster blends what they can and carries the rest.
inverse() follows solve()'s rules, because
it is solve(a, identity): exact by default up to 64 by 64,
exact: true or exact: false beyond (past that
line with no option it raises EXACTCHOICE), and a mix with no
inverse — a light roast that is just half house and half espresso,
so it adds nothing new — raises SINGULAR, with
_string$ naming the column: inverse(): a zero pivot at
column 3. For a single system, solve(a, b) is still
the tool (less work, and one rounding); inverse() is for the
matrix itself — the same one met by many right-hand sides over time,
as here.
What the array toolkit does, and where it stops. For linear algebra Sheerpower gives you the everyday tools, all in exact arithmetic:
solve(a, b [, exact:])— solve a square system (exact to 64 unknowns;exact: trueorexact: falsebeyond)lstsq(a, b [, exact:])— least-squares fit of an overdetermined system (the normal equations, in one call; the same size rule)matmul(a, b)— multiply two matrices (always exact; a 256-bit fast path for ordinary reals)transpose(a)— flip rows and columnsdeterminant(a [, exact:])— the determinant (zero = no unique solution; the same size rule)inverse(a [, exact:])— the inverse matrix, for when you need the matrix itself rather than one answer (the same size rule)maximize(c, a, b [, exact:] [, whole:])/minimize(c, a, b [, exact:] [, whole:])— the best plan: the x >= 0 that optimizes c · x under a x <= b (linear programming, the exact simplex; the same size rule);whole: truefor the best plan in whole units (branch and bound over the same simplex)maximize name = expr / [then minimize name = expr] / where names are whole | decimal / subject to / label: expr <= expr ... / end maximize(andminimize) — the same as a block statement with named unknowns, whole or decimal per name, ranked goals, one labelled constraint per line; the answer lands in the names (its own page: Resource Allocation with Maximize and Minimize)diophantine(a, c)/lattice(a)— whole-number solutions of a · x = c (or of a system a x = c): one solution, and the directions every other one lies along; exact, withNOSOLUTIONwhen the gcd does not divide cnorm(a [, kind:])— the length (Euclidean),manhattanormaxnorm of an array
Beyond these — eigenvalues, singular values (SVD), matrix
factorizations (QR, Cholesky), or very large sparse matrices
— you are into specialist numerical territory whose answers are
usually irrational and computed in floating point, so they belong in a
dedicated numerical library rather than here. And when the job is one
system, still reach for solve(a, b) rather than
matmul(inverse(a), b): it does less work and rounds once.
Vector length and distance: norm()
norm(a) is the size of an array as a single number — its
length. The default is the Euclidean (2-)norm,
sqrt(sum of squares): for a vector that is the straight-line
length from the origin, and the most common use is a distance
— norm(a - b) is the straight-line distance between two
vectors (nearest-neighbour, clustering, how far a prediction misses). It
reduces an array of any shape to one scalar, working entrywise over
every element in storage order, so a matrix or a 3-D array works too (on a
matrix the Euclidean norm is the Frobenius norm).
Two other norms are chosen with kind: — the option word
is a bare keyword:
Which depot is nearest? A parcel has to be collected, and
three depots could send a van. The parcel and each depot are points on a
map (kilometres east, kilometres north). Send the depot whose straight-line
distance to the parcel is smallest — that is a norm of
the difference, one per depot, and the smallest wins:
Did the part pass inspection? A machined part is measured
on five dimensions against its target, and it passes only if no single
dimension is off by more than the 0.20 mm tolerance. "The
worst single error" is exactly the max norm of the error
vector — one line decides it, and the manhattan norm on
the same errors gives the total drift for a process-quality log:
All three are entrywise — they treat the whole array
as one flat list of numbers, at any rank. (The specialist operator
matrix norms — the largest singular value, the max column sum —
stay out of scope with SVD.) The result is exact: the Euclidean norm rounds
once, at the final square root. A string array raises ARRAYSHAPE.
A whole batch of vectors at once: axis:.
A single call reduces the whole array to one number; add axis: n
and it reduces just one dimension, giving a 1‑D vector — the
norm of each row (or column). This is the batch case: N points or
feature‑vectors stored as an N by D matrix, and you want all N lengths (or,
with norm(a - b, axis:), all N distances) in one line. Like
reduce(), axis: n folds dimension n,
so on a matrix axis: 2 is per‑row and axis: 1 is
per‑column:
Each result is a 1‑D vector (one norm per lane), so it flows into anything:
stats$max(norm(points, axis: 2)) is the farthest point,
filter(names$, norm(offsets, axis: 2) < 5) keeps the near
ones. Without axis:, norm(m) is still the single
whole‑array number.
9. Any Function, Over a Whole Array
Every built-in math function that takes a number also takes an array,
and works on every element at once: abs(v),
round(v, 2), sqrt(v), mod(v, 3),
clamp(v, 0, 100), int(v), sin(g).
The result is a REAL array of the same shape (a 2 by 3 array stays 2 by
3). The other arguments are single values, evaluated once — or
arrays of the identical shape, paired element by element:
max(a, b), min(a, b), mod(a, b),
clamp(v, lo_arr, hi_arr). A number with an array
broadcasts: max(v, 0) clips every element at zero.
rnd(v) draws once per element.
Show me: math functions mapping over an array
One-argument max(v) is still the largest element, a
number. x = max(v, 10) with a plain numeric
x is a compile error — that form gives an array now,
so assign it to an array or use stats$max(v). Two arrays
of different counts raise ARRAYSHAPE; a function's own
exception (sqrt() of a negative) is the ordinary catchable
one.
The array need not come first: max(7, v) and
round(3.14159, decimals) map just as well — whichever
argument is the first array sets the shape, and the single values are used
for every element.
Your own routines map the same way: a function-form routine whose
first with parameter is a value, given an array there, runs
once per element — tax = calc_tax(amount = prices).
See “Mapping a routine over an array” on the Routines page,
Routine Parameters in Detail.
10. String Functions Over String Arrays
The same rule for strings: a string function given a string ARRAY as
its first argument maps element by element. A string result gives a
string array (ucase$(names$), trim$(padded$),
mid$(names$, 1, 2), replace$(names$, "a=A"));
a numeric result gives a REAL array (len(names$),
pos(names$, "b"), val(strs$)); and the other
way round, str$(v) or chr$(codes) turn a
numeric array into strings. Extra arguments are single values, or
arrays of the identical shape. Named options work inside a mapped call:
contains(names$, "b", exact: true).
Show me: string functions mapping, mixed types, and pipes
The arguments need not be of one type. A function whose parameters mix
strings and numbers — left$(string$, length),
mid$(string$, start, length),
repeat$(string$, count), lpad$(string$, width)
— takes an array for any of them, and the arrays are paired
element by element: the first word with the first length, the second with
the second. A single value in any position — the first included
— is used for every element: left$(words$, lens) cuts
each word to its own length, and left$("hello", lens) cuts
the one word to each length in turn. The first array among the
arguments gives the result its shape, so any arrays that are paired must
have that identical shape; a mismatch raises ARRAYSHAPE,
naming both.
And an array goes through the pipe like any value: |>
feeds it as the first argument, so a chain of maps reads left to right,
with the paired arrays and single values as the remaining arguments.
11. Picking by Position: x$[y]
Parentheses on an array hold a 0/1 mask (section 4); BRACKETS hold an
index list. x$[y] is the elements of x$ at
the positions listed in y, in y's order and
shape — duplicates allowed, any count, a 2 by 2 index array
gives a 2 by 2 result, and the index array may be REAL or INT. It
works on numeric and string arrays alike. On the left of
= it scatters — with one value it writes that value
into every listed position; with an ARRAY on the right it pairs the two,
writing w(k) into position y(k). That array form
is the write-by-index twin of the read: gather reads by index, scatter
writes by index. Because the source is taken as a snapshot,
x[p] = x is a correct in-place reorder — move every
element to a new spot in one line.
Show me: gather, scatter, and slices (v[a:b])
The counts must agree: if y and w have
different lengths you get a catchable ARRAYSHAPE. An index
that is not a whole number raises
NUM_OUTOFRANGE; one outside 1..size raises
SUBOUTBND — both name the position. A string array
takes a string source, a numeric one a numeric source. x$(m)
keeps its mask meaning and x$(3) is still one element; the
array-on-the-right form is only for the brackets.
A range of positions: v[a:b]
An index list picks scattered positions; a slice picks a
contiguous run. Put a colon in the brackets and you get the elements from
a to b — both ends included, the same rule
as a string slice s$[a:b]. * stands for the last
index, so v[5:*] runs to the end, v[:3] starts at
the first, and v[*-2:*] is the last three. A third number is a
step, and a negative step walks backwards:
A slice reads as a value — a copy of the run — so it
composes like any array: stats$sum(v[1:3]),
filter(v[2:*], v[2:*] > 30). On a grid, an index subscript
picks a single row or column and a range keeps that dimension, so
g[2, *] is row 2 (one-dimensional) and g[1:2, 2:3]
is a 2 by 2 tile. On the LEFT of = a slice writes through —
one value fills the run, an array pairs into it (counts must agree, else
ARRAYSHAPE):
Brackets and a colon are what make a slice: v[y] (an array,
no colon) is the index list above, v[3] is one element, and
v[a:b] is the range. Parentheses v(a:b) are not a
slice — slicing is a bracket feature. The seq() function
from the Arrays page is the movable cousin: v[seq(1, 3)] is
v[1:3], but seq() can also count backwards, skip,
or be reused as its own array.
12. Sorting: sort() and sortindex()
sortindex(x) does not sort x; it returns the
positions that WOULD sort it — a REAL array of 1-based
subscripts, the same count. Put that through the brackets and you have
the sorted array; put another array through it and you have sorted
one array by another, which is the reason it exists. The sort is
STABLE: equal keys keep their original order, also under
descending: true. Numbers compare exactly; strings are
case-blind unless nocase: false.
Two defaults, stated together so neither surprises you:
the array functions sort(), sortindex(),
unique() and isin() compare strings
case-blind by default (nocase: false for exact),
while the sort by clause inside a collect or
extract is case-sensitive by default
(sort nocase by for case-blind). The clause came first and
client programs depend on its order, so its default cannot change; the
array functions took the case-blind default of dim x(n)
sorted arrays, which is what most name lists want. When the two
meet in one program, write the option out on both sides.
Show me: sortindex(), sort(), and ranking
When it is the sorted array itself you want, sort(x) gives
it — the same options, the same stable order:
sort(x) is exactly x[sortindex(x)]. Three ways
to use it. x = sort(x) sorts x in place
— any array, a fixed, integer or two-dimensional one too, its shape
kept. y = sort(x) fills y, an expandable array
of the same kind (strings sort to a string array, numbers to a REAL
array). And inside an expression it is simply a temporary array:
print sort(x), stats$sum(sort(x)),
ucase$(sort(names$)). It is quick — a million numbers
sort in well under a tenth of a second — because plain numbers are
sorted by their digits (a radix sort) rather than by comparing pairs.
A sort where a single value is expected (total = sort(v))
is a compile error — an array has nowhere to go.
The other direction: from an order to a ranking.
sortindex reads by position —
names$[sortindex(ages)] pulls the names out in age order. The
reverse question, what place did each finish in?, writes
by position, and that is exactly the scatter from section 11
(x[y] = w). Four runners A, B, C, D with their finish times:
sortindex gives the finishing order, and scattering the places
1, 2, 3, 4 back into those slots labels each runner in their own roster
position — a leaderboard in two lines:
Reading by index and writing by index are the two halves of one idea:
sortindex tells you the order, and a scatter turns that order
into a rank for each item — the kind of standing you see on any
leaderboard.
13. Storing Through a Mask: vals(bools) = value
A mask on the left of = selects the elements to store
into: every element whose mask element is 1 receives the value, in one
statement. The mask is usually a comparison on the array itself, which
makes the clamp a one-liner. The mask must hold only 0s and 1s and the
same number of elements as the array — an index list mistaken
for a mask is loud at its first 2 (NOT0OR1) instead of
silently selecting elements 1, 2, 3; a wrong count is
ARRAYSHAPE. The read side, y = v(m), is
filter(v, m).
Show me: masked store, the clamp idiom, and a floor rule
A real rule, applied in one line. A shop takes 2 off every
price for a weekend sale — but nothing may sell below its 1.50 floor.
Subtract the discount across the whole list (the array arithmetic from
earlier), then a single masked assignment lifts just the items that fell too
far and leaves the rest untouched. _integer tells you how many
it changed:
That is what a mask is for: apply a rule to the exceptions — a floor, a cap, a correction — reaching only the elements that meet the condition, with no loop and no element-by-element test.
14. Text to Array and Back: split() and join$()
split(text$ [, delim$]) cuts a string into a string array
of its pieces; join$(array [, delim$]) writes an array
back as one string. The delimiter defaults to a comma for both and may
be any length. The pieces are verbatim — no trimming, so
trim$(split(line$)) is the cleaned form (a string map).
An empty text gives an empty array; leading, trailing or adjacent
delimiters give empty pieces; an EMPTY delimiter raises
BADFORMAT — a blank sep$ is an
unassigned variable far more often than a plan. join$()
writes numbers the way str$() does and takes an array
expression as readily as a name.
Show me: split() and join$() in code
The target of words$ = split(...) must be an expandable
string array; a plain string on the left raises
ARRAYRESULT. That statement fills the array directly
— no temporary copy — so a million-field line splits in a
fraction of a second.
15. Reshaping in Place
When the array you reshape is the array you assign to,
x = reshape(x, 1000, 10) does not build a copy: the array
changes shape where it stands, keeping its data in storage order, and a
growing recycle: true repeats the elements in place. On a
ten-million-element array that is the difference between half a second
and nothing measurable.
Part Two: Across Arrays
Sections 1 to 15 work on one array at a time. The sections from here on work across arrays — shapes that stretch to meet each other, one statistic per row or column, tables built from rows, and the set operations — and each leans on a section above. They are short; take them as needed.
16. Shapes That Stretch: Broadcasting
Section 1 said two arrays in one expression must have the same shape. There is one relaxation, and it is the useful one: a dimension of size 1 stretches to match the other array. Shapes pair up by position from the last dimension, so a 1 by 4 row against a 3 by 4 table repeats the row down every row of the table, and a 3 by 1 column repeats across every column. That is the whole rule (it is NumPy's, borrowed rather than invented), and it turns two of the most common table chores — "apply this per-column adjustment" and "scale each row by its own factor" — into one line each.
Show me: a row down a table, a column across it, and centering every column
The classic use is centering: subtract each column's mean from that
column. The means come out as a vector (section 17), and
reshape(colmean, 1, 4) lays it out as a row so it broadcasts
down the table:
A bare vector against a table is refused on purpose. NumPy reads it as a row, half of everyone expects a column, and a rule can be loosened later but never tightened — so the exception says which reshape you meant:
Broadcasting works in comparisons and masks too, and in the extra
arguments of the mapped functions of section 9:
max(sales, row) clips every column by its own floor,
clamp(sales, col, row + 10) takes a column of lows and a row
of highs at once. A single number still broadcasts everywhere, as it
always did, and a fixed target must still match the result's shape.
17. One Statistic per Row or Column: axis: and keep:
stats$sum(sales) is one number for the whole table. Add
axis: n and you get one statistic per lane along
that dimension — on a table, axis: 1 collapses the rows
and gives one value per column, axis: 2 collapses the columns
and gives one per row. Every statistic that takes an array takes the
option: sum, mean, min,
max, median, stddev,
var, percentile. Each lane runs through the
same kernel the plain call uses, so every exactness rule holds per lane.
Show me: column totals, row totals, and keep: true for centering and scaling
The result of an axis reduction is a vector, one dimension fewer than the
source. keep: true keeps the folded dimension, as size 1, so
the result is a 1 by 4 row or a 3 by 1 column that
broadcasts (section 16) straight back against the source. Centering every
column, or scaling every row to sum to 1, is then one expression with no
reshape():
reduce() and norm() take the same two options,
so a fold you write yourself behaves the same way:
An axis call is an array, so it lives on the right of an array assignment
or inside an array expression; in a scalar spot
(total = stats$sum(g, axis: 1)) the compiler says so. An axis
of 0 or past the array's rank raises ARRAYSHAPE.
18. The Dot Product: dot()
dot(a, b) multiplies two arrays element by element and adds
the products — one exact number. It replaces the older spelling
stats$sum(a * b) (same value, one intermediate array more), and it is the
calculation behind a total bill, a weighted score, and the cosine
similarity of two vectors.
Show me: a bill, a squared length, and cosine similarity
The two arrays pair in storage order and must hold the same number of
elements (else ARRAYSHAPE naming both); a 2-D array is fine,
so are slices and array expressions. matmul(v, w) still
refuses two vectors — its message points here.
19. Keeping Whole Rows: filter(g, mask, axis: n)
Section 4's filter() picks elements and always answers
with a flat list. On a table you usually want whole rows: the
orders over a limit, the readings from one sensor. axis: 1
does that — the mask has one element per row, and every chosen row
is kept with all its columns, so the result is still a table.
axis: 2 keeps chosen columns.
Show me: rows by a column test, columns by a mask, and the flat pick for contrast
The mask is usually built from one column, g[*, 2] > 4
(section 11's slices), and it must have exactly one element per
index along the axis — a wrong count raises ARRAYSHAPE
naming both, partial: true tolerates it. String tables
filter the same way.
20. String Arrays: + and the Comparisons
String arrays take + and the six comparisons.
+ concatenates element by element (section 5 lists the
moves, section 10 the functions). A single
string is joined to every element, two string arrays pair up, and numbers
join in through str$(). The shapes broadcast like the
numbers do (a 1 by c row of suffixes across a grid), and the
result feeds the string maps of section 10.
Show me: full names, titles, and a label with a number in it
A numeric array on the other side of the + is a compile error
that says str$(); shapes that do not broadcast raise
ARRAYSHAPE.
Comparisons make masks, on strings as on numbers.
names$ = "bob" is a 0/1 mask — one element per name,
the SCALAR rules per element: = and <>
byte-exact, < and the rest in the order the scalar
< uses (capitals before lower case). The other side may
be one string or a string array of the same (or a broadcast) shape.
That is the string-column filter the mask idioms of section 4 need, and
the masked store of section 13 works through it too. For a
case-blind test use same(), which maps over an array
like every string function: same(names$, "bob") is 1 for
bob and Bob alike.
Show me: filter a string column, count matches, replace a value, case-blind
A string array compared with a numeric array is a compile error
(str$() the numbers, or val() the strings);
shapes that do not broadcast raise ARRAYSHAPE.
21. Where Is the Largest? stats$argmax() and stats$argmin()
stats$max(scores) tells you the best score;
stats$argmax(scores) tells you which one — its
position, a 1-based index in storage order. Put that through the
brackets of a parallel array and you have the name that owns the score,
which is the question people actually ask. Ties go to the first
occurrence. It is a
statistics function, so it takes axis: n (the position within
each row or column) and collected: true on a cluster field
like the rest of the family.
Show me: the top scorer, the cheapest supplier, and per-column winners
On a table the plain call gives the flat storage-order position;
axis: 1 gives the winning row of every column,
axis: 2 the winning column of every row:
An empty source raises NOSTATISTIC, as stats$max
does.
22. The Distinct Values: unique()
unique(v) is v with the duplicates removed, in
first-occurrence order — the order collect's
unique uses for its representative row, and the order you
would get by hand. sort(unique(v)) when you want it sorted;
size(unique(v)) when you only want to know how many
distinct values there are. Strings compare case-blind by default
(nocase: false for exact), numbers exactly.
Show me: distinct values, the spellings that survive, and in place
The three shapes sort() has apply: x = unique(x)
compacts the array in place (an expandable array shrinks; a fixed one
must keep its count, else ARRAYSHAPE), y =
unique(x) fills y directly, and inside an expression
the result is a hidden temporary — unique(split(line$))
is the distinct fields of a CSV line. A 2-D source is read in storage
order and gives a vector.
23. Building a Table from Rows: stack()
Section 3 pours one array into a shape. stack(a, b, ...)
puts several arrays together: with no option each operand goes
under the previous one as rows — a vector counts as a row,
so three vectors make a 3-row table — and axis: 2 puts
them side by side. Two to sixteen operands, names or array
expressions, all strings or all numbers. The extent that is not being
stacked must agree (the column count when stacking rows, the row count
side by side), or ARRAYSHAPE names both.
Show me: rows into a table, a column onto a table, and vectors end to end
Vectors alone on axis: 2 concatenate into one vector, which
is what stack(v, w, axis: 2) can only mean; an empty operand
contributes nothing. With unique() it gives the union of two lists:
unique(stack(a, b, axis: 2)) (end to end, then the
distinct values; lists of different lengths cannot stack as rows).
24. Membership and the Set Operations: isin()
isin(values, set) is a mask: 1 where the element of
values occurs anywhere in set, 0 elsewhere.
That one function, fed to the mask idioms of section 4, gives every set
operation without a name for each: the elements of a also in
b is filter(a, isin(a, b)); those not in
b is filter(a, not isin(a, b)); the union is
unique(stack(a, b, axis: 2)); the count of members is
stats$sum(isin(a, b)). Strings compare case-blind by default
(nocase: false for exact), numbers exactly; the mask has
values' shape.
Show me: intersection, difference, union, and a stop-word list
Under the hood the two arrays are sorted together once, so a big list
against a big set costs a sort, not a scan per element. Duplicates on
either side are fine. For ONE value against a list use
match() or a cluster's findrow, as before.
Every Exception on This Page
Each is catchable with when exception in:
ARRAYSHAPE when two arrays in one expression do not fit
together or broadcast (sections 1, 6, 7, 16 to 19 and 23),
WRONGNUMDIMS when a subscript count does not match a
reshaped array (section 3), SINGULAR when a system has no
unique answer (section 8), NOT0OR1 for a mask element that
is not 0 or 1 (sections 4 and 13), SUBOUTBND /
NUM_OUTOFRANGE for a bad index in brackets (section 11),
BADFORMAT for an empty split() delimiter and
ARRAYRESULT for an array function assigned to a plain
variable (section 14), and NOSTATISTIC for a statistic
— stats$argmax() included — of an empty source
(section 21). None of these produce a NaN or an infinity: those are
values you put in (_nan, _inf, a CSV cell) and
the fine print below says how they behave once they are there.
- Arithmetic:
due = prices - discount
Compute a whole array at once; an expandable target takes the shape of the arrays on the right, a fixed target must already match it (else ARRAYSHAPE), andprint due + 5prints an array expression without assigning it. - Masks:
passes = scores >= 60
Comparisons andand/or/notwork per element, giving 1s and 0s you can count (stats$sum) or multiply by — exact decimal, sod * 3 = 0.3holds for 0.1. - Shapes:
reshape(v, 3, 3, recycle: true)/redim x(2, 3)
Pour an array into a new shape, or reshape an expandable array in place keeping its data; recycling the source gives identity matrices and stripes in one call.reshape(a)with no sizes flattens to a vector;redim x(*)does it in place. - Selecting:
filter(scores, scores >= 60)
The values whose mask element is 1, as a one-dimensional array. - String arrays:
grid$ = reshape(fruits$, 2, 2)
Copy,reshape(),transpose()andfilter()work on string arrays; the one operator is+(item 20). - Matrix multiply:
z = matmul(a, b)
Rows froma, columns fromb, exact sums; a vector acts as a column on the right and a row on the left. - Rows into columns:
transpose(a)
Turn a per-month table into a per-product one; the building block of the best-fit line. - Solving:
solve(a, b)/lstsq(a, b)/diophantine(a, b)/inverse(a)
When you know the totals and want the parts: prices from takings, counts from capacities.solve()answers exactly, and raisesSINGULARwhen the facts do not pin the answer down; it does not know that a count must be whole, so when it must be, calldiophantine(a, b)on the same two arrays: the same answer when it is whole,NOSOLUTIONwhen it is not.lstsq(a, b)is the best-fit (least-squares) answer for more facts than unknowns, in one call: the named form ofsolve(matmul(transpose(a), a), matmul(transpose(a), y)).inverse(a)works a matrix out once, for when the same one meets many right-hand sides. - Whole-number solutions:
diophantine(a, c)/lattice(a)
One whole-number solution of a · x = c (one equation or a system), and the directions every other solution lies along;NOSOLUTIONnames the gcd that does not divide the total. Which stamps make 47 cents, which boxes ship exactly 20. - Best plan:
maximize(c, a, b, whole: true)/minimize(c, a, b)
The plan that makes a goal largest or smallest within limits, every quantity zero or more;whole: truefor counted things. Written as a block with named unknowns on its own page (Resource Allocation with Maximize and Minimize). - Functions over arrays:
abs(v),round(v, 2),max(a, b),ucase$(names$),len(names$)
Any math or string function maps over a whole array, keeping its shape; extra arguments are one value or an identical-shape array. - Index lists and slices:
x$[y],v[y] = w,v[2:4]
Brackets pick elements by position, in the index array's order and shape; on the left they scatter — one value, or an array paired by position (sox[p] = xreorders in place). A colon makes a slice — a contiguous run (v[5:*],g[2, *]), readable and writable. - Sorting:
names$[sortindex(ages)]
The permutation that sorts an array, stable, case-blind by default; one array sorted by another. - Storing through a mask:
prices(prices < 0) = 0
A 0/1 mask on the left selects the elements to assign; a 2 in the mask is NOT0OR1, a wrong count ARRAYSHAPE. - Text and arrays:
words$ = split(line$),join$(v, " | ")
Pieces verbatim, comma by default, an empty delimiter raises; the inverse writes numbers as str$() does. - Broadcasting:
sales + row,sales * col
A dimension of size 1 stretches to match; shapes pair from the last dimension. A bare vector against a table is refused — sayreshape(v, 1, n)(row) orreshape(v, n, 1)(column). - Axis reductions:
stats$sum(g, axis: 1),keep: true
One statistic per column (axis 1) or row (axis 2);keep: truekeeps the folded dimension as size 1 so the result broadcasts back (g - stats$mean(g, axis: 1, keep: true)centers every column). Also onreduce()andnorm(). - Dot product:
dot(qty, price)
The exact sum of the products, one number;dot(v, v)is the squared length. - Whole rows:
filter(g, g[*, 2] > 4, axis: 1)
Keep the rows (axis 1) or columns (axis 2) a mask chooses, table shape kept. - String + and comparisons:
first$ + " " + last$,filter(names$, names$ = "bob")
Element-wise concatenation, and the six comparisons as 0/1 masks (byte-exact;same(names$, "bob")for case-blind); numbers join throughstr$(). - Where is the largest:
names$[stats$argmax(scores)]
The position of the maximum / minimum, first on ties; per row or column withaxis:. - Distinct values:
unique(v)
Duplicates dropped, first-occurrence order;sort(unique(v))sorted,size(unique(v))the count. - Stacking:
stack(r1, r2, r3)
Arrays as rows (a vector is a row) or side by side withaxis: 2. - Membership:
filter(a, isin(a, b))
A 0/1 mask of which values occur in a set; with filter() and unique(stack(a, b, axis: 2)) it is intersection, difference and union.
Array Math Fine Print
The exact rules behind the sections above, for when a program is near an edge. Every statement here was run, not reasoned.
-
NaN and infinity are values, never surprises.
_nanand_inf(and-_inf) exist,val("NaN")and a CSV cell readingInfproduce them, and once present they propagate by the IEEE rules (_inf - _infis NaN,1 / _infis 0, a NaN in a sum makes the sum NaN). What REAL does NOT do is manufacture them from ordinary numbers: dividing by zero, the square root of a negative and an overflow past the printed range still RAISE (DIVBY0,NUM_OUTOFRANGE,NUMOVER), so a bad value never appears silently. A NaN compares FALSE to everything, itself included, except through<>;isnan(v),isinf(v)andisfinite(v)are the masks that find them, and sorting puts NaN last. Everystats$function propagates a NaN, so the skip is explicit:stats$mean(filter(v, isfinite(v))). They print asNaN,Inf,-Inf; json$ writesnull. For rows with holes the cluster idiom still reads best:collect ... exclude isnan(c->f).div0(a, b)gives 0 for a division by zero when 0 is the answer you want, anddiv0inf(a, b)gives the IEEE one (Inf,-Inf, or NaN for0 / 0) when a hole or an infinity is — it maps over arrays, sodiv0inf(this_year - last_year, last_year)is a growth table whose zero bases come out asInfand NaN, ready forfilter(g, isfinite(g)). -
Same shape means the same number of dimensions
AND the same bounds. A 6-element vector and a 2 by 3 grid are
different shapes, and so are 2 by 3 and 3 by 2, even though each
pair holds six elements:
g + wraises ARRAYSHAPE in both cases. A one-dimensional expandable operand's shape is its appended elements. -
Vectors in
matmul(). A one-dimensional operand acts as a column on the right (result: a vector with one element per row of the matrix) and as a row on the left (one element per column). The result is one-dimensional in both cases.matmul(vector, vector)raises ARRAYSHAPE — the dot product isdot(a, b)(section 18), and the message says so. A 2D matrix with one column is not a vector:matmul(a, m)withm2 by 1 gives a 2 by 1 result, not a vector. -
What shape comes back.
filter(): always one-dimensional.reshape(a, n): one- dimensional.transpose()of an r by c matrix: c by r; of a vector: 1 by n (two-dimensional); of that: n by 1, still two-dimensional.solve(a, b): a vector whenbis a vector, n by k whenbis n by k. An array expression's result has the operands' shape. -
What
dim x(*)permits. An expandable array takes ANY shape from an assignment — the bounds and the number of dimensions — and stays expandable. After a one-dimensional result,x(*) = v(andfill x(*) with ...) appends after the last element. After a multi-dimensional result every slot is live, sox(*) = vraises SUBOUTBND (andfill x(*)FILLOVER) until the array is one-dimensional again (x = reshape(x, size(x))orredim x(*)). A fixed-shape array (dim g(2, 3)) never changes its number of dimensions:redim g(4, 5)is legal,redim g(6)andg(6)are compile errors. Routine parameters declaredvalues(*)accept an array of any shape, by reference — see the Routines page. -
filter()lengths. Values and mask pair element by element in storage order (row by row), whatever their shapes, and must hold the same number of elements: a mismatch raises ARRAYSHAPE naming both counts, like every other shape mistake on this page.partial: trueopts into the tolerant pairing (a short mask selects nothing past its end, a long mask's extras are ignored). A mask element selects on 1 only (a comparison mask holds only 1s and 0s). -
Maps keep the shape and always give REAL.
abs(g)of a 2 by 3 array is 2 by 3;len(names$)is a REAL array even though each length is whole;pos(),contains()and the other INT-result functions come back REAL too. A second ARRAY argument must match the first in bounds AND dimension count (for an expandable 1D array, the appended count) — not merely the element count — else ARRAYSHAPE names both; there is nopartial:form for maps. -
Brackets and parentheses never collide. With an
ARRAY subscript,
x(m)is a mask (0/1, same count) andx[y]an index list (positions, any count, any shape); with ONE NUMBER, both are the element. A scatter stores one scalar, or an ARRAY paired with the index list by position (x[y] = wwritesw(k)intoy(k), sox[p] = xreorders in place); repeated positions just store again. An INT array anywhere in an array expression is read through a REAL copy. -
sortindex()versussorted.dim x(100) sortedkeeps ONE array in order as you write it;sortindex()answers the question "in what order?" without moving anything, so several arrays can follow one key. Its string default is case-blind (nocase: true, the sorted-array rule) — note thesort byCLAUSE in a collect is case-sensitive unless you writesort nocase by. A 2D source gives storage-order positions in a 1D result. -
Mixed-type maps: the first argument decides. A
built-in maps when its FIRST argument is an array (a string array for
the string functions, a numeric one for the math functions); every
later argument may be a single value or an array of the identical
shape, of whatever type that parameter takes —
left$(words$, lens)pairs strings with numbers,round(v, decimals)numbers with numbers,pos(words$, needles$)strings with strings. A single value in the first position with an array later maps too (left$("hello", lens),max(7, v)): the first array argument, wherever it sits, sets the shape and the single values broadcast. Arrays pipe:words$ |> left$(lens)isleft$(words$, lens). -
sort()issortindex()applied. Same order, same options, same stability.x = sort(x)moves the elements insidex's own block, so the shape is kept and a fixed, INT or 2D array (storage order) sorts in place too;y = sort(x)fillsydirectly, and the kinds must agree (a string source needs a string target; an INT source may fill a REAL one); a 2D source sorted into ANOTHER array gives a 1D result of every element. Anywhere else the result is a temporary of the source's kind. The engine sorts plain numbers with a radix sort on their digits (no comparisons at all), falls back to a stable merge when a value is WIDE or SCI, and skips the work when the array is already in order. -
split()is a snapshot. The pieces are copies; changing the text afterwards does not change them (aVIEW ... pieceis the live alternative for one field).split(s$)where a string function's array argument is expected works everywhere an array expression works —ucase$(split(s$)),size(split(s$)),print split(s$), insidef$.join$(g, "/")of a 2D array walks storage order;join$of an empty array is"";join$(w$, "")concatenates (an empty delimiter is fine on this side). Not the same as the olderjoin(), which concatenates its arguments INTO a target and returns the length. -
SINGULAR means no UNIQUE solution. It is raised
both when infinitely many answers fit (Tuesday's sales AND
takings exactly double Monday's — the facts are consistent
but say the same thing twice) and when no answer fits at all
(double the sales but not double the takings).
solve()cannot tell the two apart; both are a zero pivot. Wrong sizes (a non-squarea, abwith the wrong number of rows) are ARRAYSHAPE, not SINGULAR.
All told, Sheerpower gives you 17 core array operations, made open-ended by array-mapping over the whole function library (built-in and user-defined), with slicing/indexing and the fill/print/structure machinery around them. The list below shows every one, all in one place.
(Show/Hide the Full List of Array Features)The Full List of Array Features
The 17 core operations
Everything here returns (or reduces) an array, and every one is exact.
matmul(a, b)— matrix multiply. R×K by K×C gives R×C, with exact sums. A vector stands for a column on the right and a row on the left, somatmul(m, v)andmatmul(v, m)give vectors; two vectors raiseARRAYSHAPE(say which is the row: for the dot product usedot()).dot(a, b)— the dot product, a single exact number: the sum ofa(i) * b(i)over two arrays of the same count (dot(v, v)is the squared length).reshape(a, d1 [, d2 ...] [, recycle: true])— pour a's elements into a new shape;reshape(a)with no sizes flattens to a 1-D vector.transpose(a)— swap rows and columns (r×c becomes c×r).solve(a, b)— solve the square linear systemmatmul(a, x) = bexactly.lstsq(a, b)— the least-squares best fit of an overdetermined (full-rank) system.determinant(a)— the determinant of a square matrix; 0 means there is no unique solution.inverse(a [, exact:])— the inverse of a square matrix:matmul(a, inverse(a))is the identity; the same exact rules and size line assolve().norm(a [, kind: manhattan|max] [, axis: n] [, keep: true])— vector length (the Frobenius norm of a matrix);kindgives the Manhattan (sum of |x|) or max norm;axisgives a vector of per-row / per-column norms.
Deliberately not here: eigenvalues, SVD, QR and Cholesky, which are iterative floating-point results with no exact answer: reach for a numerical library for those.outer(a, b, op)— every element of a paired with every element of b throughop(the times / addition table);opis a bare operator, a built-in, or your own routine.reduce(a, op [, axis: n] [, keep: true])— fold an array down to a scalar, or one dimension smaller along an axis: on a gridreduce(g, +)is the row sums (the last axis),reduce(g, +, axis: 1)the column sums; the same with*,max(),min()or your own routine.keep: truekeeps the folded dimension as size 1, so the result broadcasts straight back against the source.stats$sum(a, axis: n [, keep: true])— andstats$mean,min,max,median,stddev,var,percentile: one statistic per row or column (axis: 1= per column,axis: 2= per row), the same exact answer the whole-array call gives each lane.g - stats$mean(g, axis: 1, keep: true)centres every column in one line. Withoutaxis:astats$*()call takes the whole array as one population.scan(a, op [, axis: n])— the running (cumulative) fold, in a's shape: running totals per row withaxis.filter(values, mask [, axis: n] [, partial: true])— the values whose mask element is 1, as a 1-D array (storage order); withaxis: nthe mask picks whole rows or columns:filter(g, g[*, 1] > 100, axis: 1)keeps the rows whose first value exceeds 100, all their columns, andaxis: 2keeps columns.sort(a [, descending: true] [, nocase: false])— the sorted array.sortindex(a [, descending: true] [, nocase: false])— the permutation vector that sorts a (sort one array by another:names$[sortindex(ages)]).unique(a [, nocase: false])— the distinct elements in first-occurrence order (sort(unique(a))sorted;size(unique(a))the count).stats$argmax(a [, axis: n])/stats$argmin(a)— the POSITION of the maximum / minimum (first on ties);names$[stats$argmax(scores)].stack(a, b, ... [, axis: 2])— the arrays one under the other as rows (a vector is a row), or side by side.isin(a, set [, nocase: false])— a 0/1 mask of which elements of a occur in set:filter(a, isin(a, b))is the intersection,filter(a, not isin(a, b))the difference,unique(stack(a, b, axis: 2))the union.split(text$ [, delim$])— a string into a string array.join$(a [, delim$])— an array into one string.seq(lo, hi [, step])— the numbers lo..hi as a REAL array (APL's iota).
Array-mapping — what makes it open-ended
You are never limited to the 17. Any function applies element-wise to a whole array, so the entire function library — and anything you write yourself — works array-wide.
- Math over arrays:
abs(v),round(v, 2),sqrt(v),sin(v),mod(v, 3),clamp(v, lo, hi),int(v),exp(v),log10(v), ... — any built-in whose first argument is a REAL, giving a same-shape REAL array. - Strings over arrays:
ucase$(names$),trim$(v$),mid$(v$, 1, 2),replace$(v$, "a=A"),len(v$),pos(v$, "x"),val(v$),str$(v),chr$(codes), ... — string built-ins over a string array (a string→number function gives a REAL array). The operators work element by element too:+concatenates (names$ + "!",first$ + " " + last$,first$ + " (" + str$(ages) + ")") and the six comparisons give a 0/1 mask with the scalar rules (names$ = "bob",names$ < "m";same(names$, "bob")for case-blind) — sofilter(names$, names$ = "bob")filters a string column andnames$(names$ = "old") = "new"replaces a value.sort(),unique(),isin()andfilter()take string arrays; join the whole array into one string withjoin$(v$). - Your own routines:
calc_tax(amount = prices)— a function-form routine whose first parameter is a scalar maps element-wise; pipes too:prices |> calc_tax(). - Extra arguments are single values (used for every
element) or same-shape arrays (paired element by element):
max(a, b),clamp(v, lo, hi).
Slicing & indexing
v[a:b]— a slice (inclusive both ends);v[a:*]to the end,v[:b]from the start,v[a:b:step]every step-th (v[1:*:2]every other; a negative step walks down,v[*:1:-1]reverses);*is the last index and*-ncounts back from it, sov[*-1:*]is the last two andv[*-2:*]the last three (there is nov[-1]: a negative subscript is a real index ondim a(-10 to 30)). On a gridg[2, *]is a row,g[*, 3]a column,g[1:2, 2:3]a tile.x[y]— gather: the elements at the positions in index array y, shaped like y.x[y] = w/x[y] = value— scatter: write w's elements (or one value) to those positions.v(mask)— masked read (the same asfilter);vals(mask) = value— masked store (a 1 selects each element to change).
Element-wise operators
- Arithmetic:
+ - * / ^and unary-— whole arrays, element by element. A single value broadcasts (v * 2,g - stats$mean(g)), and so do shapes, by one rule: shapes pair by position from the last dimension, and a dimension of 1 stretches. So a 1×c row repeats down the rows of an r×c matrix and an r×1 column repeats across the columns — per-column centering isg - reshape(colmeans, 1, c), a per-row offsetg + reshape(rowoffs, r, 1). A bare vector against a matrix is refused (say which way it lies:reshape(v, 1, n)for a row,reshape(v, n, 1)for a column); any other mismatch raisesARRAYSHAPEnaming both shapes. - Comparisons give 0/1 masks:
= <> < <= >= >—passes = scores >= 60. A mask is a boolean array, kept one bit a flag (a REAL or integer array of 0s and 1s is accepted as a mask too); a mask element that is not 0 or 1 raisesNOT0OR1, so an index list mistaken for a mask is loud.stats$sum(mask)counts the 1s — the way to count the elements that meet a condition (Array Math, section 2). - Logic on masks:
and,or,not.
Statistics
The whole stats$*() family (48 aggregates —
stats$sum, stats$mean,
stats$median, stats$stddev,
stats$percentile, stats$pcorr,
stats$min, stats$max, ...) works on any
array.
Declaring, filling & shaping
dim x(n)/dim g(r, c)/dim a(lo to hi)— fixed arrays (1-origin; custom or negative bounds).dim x(*)— expandable;x(*) = vappends,x(*)reads the last, and(*)is the append cursor of any array.redim x(...)— resize (keeps data) or reshape an expandable array.fill x with ...— load in one statement: broadcast, a list,n of v,seq(a, b), positionedx(pos), or appendx(*).print x/print array x: list, index, csv, quoted— print an array (quotedwraps each string element in quotes, revealing empty / whitespace-bearing elements).dim x(n) sorted— a self-maintaining sorted array.freeze x/permafreeze x— make a whole array immutable.routine r with values(*), returning sq(*)— pass arrays by reference to routines.x = y— whole-array copy;typeof$(x)— an array's shape and state.lbound(x [, dim])/ubound(x [, dim])— the low and high index of a dimension (dim grid(2 to 5, 0 to 3):lbound(grid, 2)is 0,ubound(grid, 1)is 5); with no dimension number they answer for the first dimension, whatever the shape. For an expandable arrayubound(x)is the last index appended, the same assize(x)when it starts at 1 (0 while nothing has been appended, so a loop fromlbound(x)toubound(x)runs zero times).maxsize(x)— how many elements the array's current allocation holds before it must grow (a fixed array: its size; an expandable one: at leastsize(x)).finditem(names$, 'fig' [, nth [, method]])— the index of an element in a STRING array, 0 when absent; case-blind unlessmethodis 1;nthpicks the nth match ('pear', 'apple', 'fig', 'Apple':finditem(x$, 'apple')is 2,finditem(x$, 'apple', 2)is 4,finditem(x$, 'Apple', 1, 1)is 4). A binary search over a sorted index the array keeps, so it stays fast on a large array;_integeris the match's rank in sorted order. For a numeric array, or a whole mask of hits, useisin()/where().
Array Math Takeaways
- A REAL array takes a whole array expression: x = prices * 1.12, x = x + 1, x = a - b (exact REAL, element-wise; every array in one expression needs the same shape, else catchable ARRAYSHAPE). An EXPANDABLE target is reshaped to the result's shape; a FIXED target (dim x(5), dim g(2, 3)) works too but is not resized -- its shape must already match the result, else ARRAYSHAPE. dim x(*) is expandable AND shape-flexible: the result's shape (even 2D or 3D) becomes x's, and redim x(2, 3) / redim x(*) reshape it keeping the data in storage order (a fixed-shape array can be resized, not reshaped). A wrong subscript count on a reshaped array is a runtime, catchable WRONGNUMDIMS. reshape(a, dims...) pours an array into a new shape (zero-padded; recycle: true repeats the source: identity matrices, stripes); reshape(a) with no sizes FLATTENS to a 1-D vector (ravel), the readable form of reshape(a, size(a)), numeric or string. Comparisons (= <> < <= >= >) and and/or/not work element-wise, each element 1 or 0 -- a MASK: passes = scores >= 60, stats$sum(scores >= 60) counts them, (mask) * values keeps the matching values; a fixed-shape array is checked when the program compiles. matmul(a, b) is MATRIX multiply (2 x 3 by 3 x 2 gives 2 x 2; exact sums; ARRAYSHAPE if the inner sizes differ; a vector acts as a column on the right, a row on the left). print z + 1 prints a whole-array expression without assigning it first.
- solve(a, b) returns the z with matmul(a, z) = b -- the parts behind known totals (prices from takings, vehicle counts from capacity and cost). b may be a vector or n x k for k systems at once. Fraction-free elimination in exact decimal: a terminating answer is EXACT (2.5, not 2.4999999999999999). transpose(a) swaps rows and columns; a best-fit line is solve(matmul(transpose(a), a), matmul(transpose(a), y)) -- or, named, lstsq(a, y): full-rank overdetermined least squares in one call, the same exact arithmetic, a square a reducing to solve. A singular / rank-deficient / underdetermined system raises catchable SINGULAR; shape mistakes raise ARRAYSHAPE.
- filter(values, mask) returns a one-dimensional array of the values whose mask element is 1 -- filter(scores, scores >= 60). The values and the mask must have the SAME count -- unequal counts raise ARRAYSHAPE (pass partial: true for the tolerant pairing, where a short mask selects nothing past its end). String arrays take the three shape moves -- x$ = y$, reshape(), transpose(), nested freely -- filter() with a numeric mask, and + (element-wise concatenation, section 20). Every array function works in print, in a function argument, in a pipe, and as a routine argument.
- Any function maps over a whole array: abs(v), round(v, 2), max(a, b) (identical shapes, else ARRAYSHAPE), ucase$(names$), len(names$) (a REAL array), str$(v) -- same shape, always REAL for numbers. Extra arguments may be arrays too, of any type the parameter takes, paired element by element (left$(words$, lens), round(v, decimals)); any argument may be the array -- a single value anywhere broadcasts (left$("hello", lens), max(7, v)). Arrays pipe: words$ |> left$(lens), v |> abs() |> round(1). Brackets are an INDEX LIST: x$[y] picks by position in y's order and shape, v[y] = 9 scatters; parentheses stay the 0/1 mask, and vals(bools) = 9 STORES through a mask (a 2 in the mask is NOT0OR1). sortindex(x [, descending: true] [, nocase: false]) is the STABLE permutation that sorts x -- names$[sortindex(ages)] sorts one array by another; sort(x [, the same options]) is the sorted array itself: x = sort(x) sorts in place (a fixed, INT or 2D array too, shape kept), y = sort(x) fills y, elsewhere it is a temporary of the source's kind. words$ = split(line$ [, delim$]) fills a string array with the verbatim pieces (default ",", empty delimiter raises BADFORMAT, empty text = empty array, direct into the target); join$(array [, delim$]) is the inverse (numbers as str$ writes them). x = reshape(x, ...) reshapes in place.
