cvec ← larg (fn ##.nats) rarg ⍝ Natural number arithmetic.
Nats applies its operand function between natural (non-negative, whole) number
arguments. The operand may be one of:
+ sum
- difference
× product
÷ quotient
* exponent
⌊ min
⌈ max
| residue
∨ greatest common divisor
∧ least common multiple
<≤=≥>≠ relational
Arguments [larg] and [rarg] are typically vectors of characters '0'-'9' but may
be numeric for (natural) numbers with no more than ⎕PP digits.
NB: [nats] deals only in single (scalar) numbers, rather than whole arrays.
The result is a character vector of decimal digits.
1234567890 ×nats '12345678901234567890'
15241578751714678875019052100
[nats] is similar to operator →big← except that it:
- deals only in non-negative numbers,
- accepts only dyadic operand functions,
- is quicker for large numbers.
Operand (-) takes the "elementary school" approach of returning 0 if the minuend
is smaller than the subtrahend (larg<rarg).
2 -nats 3 ⍝ larg < rarg
0
Operand (÷) returns the floor of the real quotient:
100 ÷nats 13 ⍝ floor of real quotient.
7
13 |nats 100 ⍝ residue (remainder).
9
9 +nats 7 ×nats 13 ⍝ dividend = remainder + quotient × divisor
100
Gianluigi Quario comments that it is difficult to enter large-magnitude numbers:
456 +nats 123E20 ⍝ raw E-notation won't work.
Bad number: 1.23E22
456+nats 1.23E22
∧
456 +nats '12300000000000000000000' ⍝ this works but is ungainly.
12300000000000000000456
456+nats 123×nats 10*nats 20 ⍝ OK but inconvenient.
12300000000000000000456
Gianluigi suggests this function:
pfmt←{⎕IO ⎕PP←0 16 ⍝ plain format of a numeric scalar.
rep←⍕⍵
⍵≡rep:⍵
~'E'∊rep:rep
exp←⍎(1+rep⍳'E')↓rep
0>exp:'_ '~⍨(⎕PP-exp)⍕⍵
sig←'¯'=⊃rep
man←⎕PP↑'.'~⍨(sig↓(rep⍳'E')↑rep),⎕PP⍴'0' ⍝ mantissa
(sig⍴'¯'),(1+exp)↑man,exp⍴'0'
}
Then:
456 +nats pfmt 123E20
12300000000000000000456
456 +nats pfmt 1.23E22
12300000000000000000456
Technical notes:
Performance benefits from representing the numbers internally in a fairly large
base (radix). However, if too large a base is chosen, intermediate results run
the risk of exceeding the 53-bit precision of IEEE floating point arithmetic and
so losing accuracy. A base of 10*6 seems optimal. A power of 10 base is chosen
to simplify conversion to and from internal format.
Product
-------
Roger Hui points out that however small a base we choose, the product of suffic-
iently long vectors of digits will cause 53-bit overflow. To avoid this loss of
precision, inner function "mul" tests the lengths of its arguments and, if nec-
essary, splits them into sub-products of shorter vectors. If ⍺ and ⍵ are vectors
of radix-rx digits, and + and × are vector sum and product respectively, then:
¯ ¯
⍺ × ⍵ ←→ ((⍺ × s↑⍵),s↓⍵⍴0) + ⍺ × s↓⍵
¯ ¯ ¯ ¯
For example, suppose a 7-digit product were enough to cause overflow. We could
avoid the overflow by splitting the product into the sum of a 3-digit and a 4-
digit product, like so:
1 2 × 3 4 5 6 7 8 9 ←→ ((1 2 × 3 4 5),0 0 0 0) + 1 2 × 6 7 8 9
¯ └─────7─────┘ ¯ └─3─┘ ¯ ¯ └──4──┘
As a test case, Roger observes that the square of natural number ⍵/'9' is
(⍵ 2 ⍵ 2-1)/'9801':
×nats⍨'9'
81
×nats⍨'99'
9801
×nats⍨'999'
998001
×nats⍨'99999999'
9999999800000001
With nat's radix of 10*6, the square of 999...9 is sufficent to require split-
ting: └60000┘
(×nats⍨6e4/'9')≡(6e4 2 6e4 2-1)/'9801' ⍝ (takes around a minute).
1
See →xtimes← for a faster multi-digit product.
Exponent
--------
For ⍺*⍵, we avoid calculating ⍵ products ⍺ × ⍺ × ···, by using the recurrence
relation: └─── ⍵ ───┘
⍺*0 = 1
⍺*1 = ⍺
⍺*⍵ = (⍺*2|⍵) × {⍵×⍵} ⍺*⌊⍵÷2
which needs only O(2⍟⍵) products.
(muse:
Notice in the examples below that */5/2 is a number with 19,729 decimal dig-
its. An amusing way to distinguish mathematicians from computer programmers
at a social gathering (perhaps a cocktail party or wedding reception) is to
sing this little song (to the tune of "The Grand old Duke of York"):
Two,
Two-to-the-two,
Two-to-the-two-to-the-two,
Two-to-the-two-to-the-two-to-the-two,
Two-to-the-two-to-the-two-to-the-two-to-the-two,
...
You will find that after six or seven lines, the mathematicians begin to
look bored, while the computer programmers turn pale and start to faint.
See also: staples
)
Quotient
--------
We use the traditional "high-school" long division method (invented by Charles
Babbage, see: eval.dws/notes.Long_division). A minor refinement is that, as the
radix is large, it pays to do a binary, rather than linear search for the next
"digit" of the result. With a radix of 10*6 (around 2*20), we can find the digit
in about 20 comparisons when using (0, radix-1) as the lower and upper search
bounds.
But we can do a little better. We take the two most significant digits of the
dividend and divisor ((p0 p1 ···)(q0 q1 ···)) and decode each pair using the
radix. This gives a pair of numbers (pp qq), each in the range 0,radix*2, which
is still well within our required 53-bit precision limit.
To find a lower bound, we assume that all digits following p0 and p1 are 0 and
that all those following q0 and q1 are (radix-1). This gives a lower bound of
(⌊pp÷qq+1).
Likewise, an upper bound assumes all digits following p0 and p1 are (radix-1)
and all those following q0 and q1 are 0, giving a upper bound of (⌈(pp+1)÷qq).
It makes things a little easier if we soften this upper bound very slightly by
observing that:
( ⌈⍵) = ⌊1+⍵ ⍝ for non-integer ⍵
(1+⌈⍵) = ⌊1+⍵ ⍝ for integer ⍵
so we can reasonably take an upper bound of: (⌊1+(pp+1)÷qq).
In other words, the lower-upper bound pair is (⌊pp÷qq+1) (⌊1+(pp+1)÷qq).
Then, using expression transformation to simplify:
(⌊pp÷qq+1) (⌊1+(pp+1)÷qq) ⍝ (⌊⍺)(⌊⍵) = ⌊⍺ ⍵
_¯ ¯
= ⌊(pp÷qq+1) (1+(pp+1)÷qq) ⍝ (0+⍺)(1+⍵) = 0 1+⍺ ⍵
____ ¯¯
= ⌊0 1+(pp÷qq+1) ((pp+1)÷qq) ⍝ (⍺÷⍵)(∊÷⍴) = (⍺ ∊)÷(⍵ ⍴)
¯ ¯
= ⌊0 1+(pp(pp+1))÷(qq+1)qq ⍝ ⍺(⍺+1) = ⍺+0 1, if ⍺ scalar.
¯¯¯¯¯¯¯¯ ¯¯¯¯¯¯¯¯
= ⌊0 1+(pp+0 1) ÷ qq+1 0 ⍝ ⍺÷⍵ = ↑÷/⍺ ⍵
¯
= ⌊0 1+↑÷/(pp+0 1) (qq+1 0) ⍝ (⍺+⍵)(∊+⍴) = ⍺ ⍵+∊ ⍴
¯ ¯
= ⌊0 1+↑÷/pp qq+(0 1)(1 0) ⍝ thus (pp qq) can use a single name ppqq.
Examples:
⎕d ×nats ⎕d ⍝ 123456789 * 2
15241578750190521
2*nats 1024 ⍝ 2*1024
17976931348623159077293051907890247336179769789423065727343008115773267580550096
31327084773224075360211201138798713933576587897688144166224928474306394741
24377767893424865485276302219601246094119453082952085005768838150682342462
88147391311054082723716335051068458629823994724593847971630483535632962422
4137216
⍴↑*nats/5/2 ⍝ */5/2 has 19,729 decimal digits.
19729
⎕d *nats 8 ⍝ 123456789 * 8
53965948844821664748141453212125737955899777414752273389058576481
↑÷nats/⎕d ⎕d *nats¨ 101 100 ⍝ test accuracy
123456789
0 1{⍵=0:⊃⍺ ⋄ (1↓⍺,+nats/⍺)∇ ⍵-1}1000 ⍝ 1,000th Fibonacci number (JRC).
43466557686937456435688527675040625802564660517371780402481729089536555417949051
89040387984007925516929592259308032263477520968962323987332247116164299644
0906533187938298969649928516003704476137795166849228875
factorial ← ↑∘(×nats/)∘⍳ ⍝ !⍵
factorial 1000 ⍝ !1000
40238726007709377354370243392300398571937486421071463254379991042993851239862902
05920442084869694048004799886101971960586316668729948085589013238296699445
90997424504087073759918823627727188732519779505950995276120874975462497043
60141827809464649629105639388743788648733711918104582578364784997701247663
28898359557354325131853239584630755574091142624174743493475534286465766116
67797396668820291207379143853719588249808126867838374559731746136085379534
52422158659320192809087829730843139284440328123155861103697680135730421616
87476096758713483120254785893207671691324484262361314125087802080002616831
51027341827977704784635868170164365024153691398281264810213092761244896359
92870511496497541990934222156683257208082133318611681155361583654698404670
89756029009505376164758477284218896796462449451607653534081989013854424879
84959953319101723355556602139450399736280750137837615307127761926849034352
62520001588853514733161170210396817592151090778801939317811419454525722386
55414610628921879602238389714760885062768629671466746975629112340824392081
60153780889893964518263243671616762179168909779911903754031274622289988005
19544441428201218736174599264295658174662830295557029902432415318161721046
58320367869061172601587835207515162842255402651704833042261439742869330616
90897968482590125458327168226458066526769958652682272807075781391858178889
65220816434834482599326604336766017699961283186078838615027946595513115655
20360939881806121385586003014356945272242063446317974605946825731037900840
24432438465657245014402821885252470935190620929023136493273497565513958720
55965422874977401141334696271542284586237738753823048386568897646192738381
49001407673104466402598994902222217659043399018860185665264850617997023561
93897017860040811889729918311021171229845901641921068884387121855646124960
79872290851929681937238864261483965738229112312502418664935314397013742853
19266498753372189406942814341185201580141233448280150513996942901534830776
44569099073152433278288269864602789864321139083506217095002597389863554277
19674282224875758676575234422020757363056949882508796892816275384886339690
99598262809561214509948717012445164612603790293091208890869420285106401821
54399457156805941872748998094254742173582401063677404595741785160829230135
35808184009699637252423056085590370062427124341690900415369010593398383577
79394109700277534720000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000
binomial ← {(factorial ⍵)÷nats(factorial ⍺)×nats factorial ⍵-⍺} ⍝ ⍺!⍵
1000 binomial 2000
20481516269894897143351625029808250443964248879813970338203826376717481862020837
55828932994182610206201464766319998023692415481798004524792018047549769261
57856301289663432064714851152395251651227768588611539546256147907378668464
15444453361761377007385567381458963007130651045595951447988874620636871851
45518285511731662762536637730846829322553890497438594814317550307837964443
70810085163724827462791417016619883764840843541430817785947037746565188475
51468074969467492380303310181872329800966856745856025254991011811352535346
58887941966653674904511306110096311906270342502293155911108976733963991149
120
⍝ Handy prefixes for the first few 2*10×⍵:
⍝
prefs ← 'kilo' 'mega' 'giga' 'tera' 'peta' 'exa' 'zetta' 'yotta' '...'
⍉↑prefs{⍺ ⍵}2*nats¨ 10 20 to 200
kilo 1024
mega 1048576
giga 1073741824
tera 1099511627776
peta 1125899906842624
exa 1152921504606846976
zetta 1180591620717411303424
yotta 1208925819614629174706176
... 1237940039285380274899124224
1267650600228229401496703205376
1298074214633706907132624082305024
1329227995784915872903807060280344576
1361129467683753853853498429727072845824
1393796574908163946345982392040522594123776
1427247692705959881058285969449495136382746624
1461501637330902918203684832716283019655932542976
1496577676626844588240573268701473812127674924007424
1532495540865888858358347027150309183618739122183602176
1569275433846670190958947355801916604025588861116008628224
1606938044258990275541962092341162602522202993782792835301376
See also: numbers xtimes xpower big gcd rats bt adic
Back to: contents
Back to: Workspaces