rational arithmetic
Forum rules
This forum is for discussing APL-related issues. If you think that the subject is off-topic, then the Chat forum is probably a better place for your thoughts !
This forum is for discussing APL-related issues. If you think that the subject is off-topic, then the Chat forum is probably a better place for your thoughts !
1 post
• Page 1 of 1
rational arithmetic
by Roger|Dyalog on Wed Jul 08, 2020 4:17 am
The note describes a facility for extended-precision rational arithmetic. It is expected that changes to the text will be made as improvements are devised and bugs (if any☺) are found. In fact, the original design was changed, informed by the analysis in A Design Exercise. See also the APL Chat Forum post Rational Arithmetic Implementation Notes.
The facility provides a representation of rational numbers, a monadic operator Q for working with rational numbers, and a function qi (“rational input”) to specify rational numbers conveniently. Some examples:
(Misalignments in the display are due to defects in the APL Chat Forum software.)
After I implemented the current facility, I discovered the operator “rats” http://dfns.dyalog.com/n_rats.htm in the dfns workspace. The current facility is more general and more capable than in “rats”.
Representation
A rational number is a scalar of an enclosed 3-element vector, whose items are:
Because a rational number so represented is a scalar, “structural” (⍴, ↑, etc.) and “selection” (⌿, [;], ⌷, etc.) functions just work on arrays of rational numbers.
For example:
The Function qi
The function qi (“rational input”) provides a convenient way to specify rational numbers. It applies independently to items of the argument; each item can be one of several types:
Examples:
The Operator Q
f Q derives a function which works on rational numbers and produces a rational result except as noted otherwise.
Scalar Functions
The boolean functions < ≤ = ≠ ≥ > return a boolean result, an ordinary 0 or 1 instead of a rational 0 or 1. Likewise ×⍵ returns an ordinary ¯1, 0, or 1.
Non-Scalar Functions
⊤Q ⍵ applies to rational scalar positive integer ⍵ and produces its binary representation. (⊤ does not have a monadic definition. One possibility is ⊤⍵ ↔ 2⊤⍵, that is, bits or binary representation. (Similarly, ⊥⍵ would be 2⊥⍵, binary value.) This is in fact the definition used in A Formal Description of SYSTEM/360 [Falkoff, Iverson, and Sussenguth 1964, Table 1]. It was actually implemented in APL\360 for a time but was later removed due to space limitations [Hui 1992]. For example, ⊤Q qi 123 ↔ 1 1 1 1 0 1 1.
⍺ ⊤Q ⍵ applies to rational scalar positive integer ⍵ and produces its base-⍺ representation, for an ordinary integer ⍺. For example, 5 ⊤Q qi 1234567 ↔ 3 0 4 0 0 1 2 3 2.
⍕Q ⍵ formats rational array ⍵, with an “r” separating the numerator and the denominator (if not 1).
⍺ ⍕Q ⍵ formats rational array ⍵ in fixed-point (if ⍺≥0) or exponential (if ⍺<0) format. The left argument ⍺ must be a single integer: ⍺≥0 specifies the number of digits after the decimal point; ⍺<0 specifies the number of decimal digits in exponential notation.
⌹Q is analoguous to ⌹ for rational arguments. As usual, a DOMAIN ERROR is signalled if ⍵ is non-singular (but not as usual, if the error is signalled ⍵ really is singular).
Applications
Computing with rational numbers opens new vistas.
Continued Fractions
How about that: our friends the Fibonacci numbers.
Square Root via Newton’s Method
If W is a number whose square root is sought and ⍵ is an approximation of the root, Newton’s method improves the approximation by computing the average of ⍵, and W divided by ⍵. By using rational arithmetic, a sequence of rational approximations obtains.
The following code displays the rational approximations to 40 decimal places, and the differences between 2 and the squares of the approximations. The differences illustrate that Newton’s method doubles the number of correct digits with each iteration.
Matrix Inverse
Hilbert’s matrix is known to be invertible, but is notorious for being nearly singular.
Inverting the 8-by-8 Hilbert matrix as 64-bit floating point numbers:
The Code
Tests
These tests provide a prima facie case that the code actually works. qtest ignores its argument and returns a result of 1 if all tests passed.
⍕Q has a separate set of tests because implementing the functionality was ... interesting.
created:_2020-07-07
updated: 2020-08-24
The facility provides a representation of rational numbers, a monadic operator Q for working with rational numbers, and a function qi (“rational input”) to specify rational numbers conveniently. Some examples:
qi 3.5
┌───────┐
│┌─┬─┬─┐│
││7│2│1││
│└─┴─┴─┘│
└───────┘
(qi 3.5) ×Q (qi ¯120)
┌────────────┐
│┌─────┬─┬──┐│
││4 2 0│1│¯1││
│└─────┴─┴──┘│
└────────────┘
3.5ׯ120
¯420
!Q qi 40
┌─────────────────────────────────────────────────────────────────────────────────────────────────────┐
│┌───────────────────────────────────────────────────────────────────────────────────────────────┬─┬─┐│
││8 1 5 9 1 5 2 8 3 2 4 7 8 9 7 7 3 4 3 4 5 6 1 1 2 6 9 5 9 6 1 1 5 8 9 4 2 7 2 0 0 0 0 0 0 0 0 0│1│1││
│└───────────────────────────────────────────────────────────────────────────────────────────────┴─┴─┘│
└─────────────────────────────────────────────────────────────────────────────────────────────────────┘
⍕Q !Q qi 40
815915283247897734345611269596115894272000000000
!40
8.15915E47
+Q\ qi 3 1 4 1 ¯5.9
┌───────┬───────┬───────┬───────┬───────────┐
│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌───┬───┬─┐│
││3│1│1│││4│1│1│││8│1│1│││9│1│1│││3 1│1 0│1││
│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└───┴───┴─┘│
└───────┴───────┴───────┴───────┴───────────┘
+\ 3 1 4 1 ¯5.9
3 4 8 9 3.1
(Misalignments in the display are due to defects in the APL Chat Forum software.)
After I implemented the current facility, I discovered the operator “rats” http://dfns.dyalog.com/n_rats.htm in the dfns workspace. The current facility is more general and more capable than in “rats”.
Representation
A rational number is a scalar of an enclosed 3-element vector, whose items are:
- a vector of the decimal digits of the numerator
- a vector of the decimal digits of the denominator
- the sign, ¯1, 0, or 1
Because a rational number so represented is a scalar, “structural” (⍴, ↑, etc.) and “selection” (⌿, [;], ⌷, etc.) functions just work on arrays of rational numbers.
For example:
Hilbert←{÷1+∘.+⍨⍳⍵}
Hilbert 5 ⍝ an array of floating point numbers
1 0.5 0.333333 0.25 0.2
0.5 0.333333 0.25 0.2 0.166667
0.333333 0.25 0.2 0.166667 0.142857
0.25 0.2 0.166667 0.142857 0.125
0.2 0.166667 0.142857 0.125 0.111111
qi Hilbert 5 ⍝ an array of rational numbers
┌───────┬───────┬───────┬───────┬───────┐
│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│
││1│1│1│││1│2│1│││1│3│1│││1│4│1│││1│5│1││
│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│
├───────┼───────┼───────┼───────┼───────┤
│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│
││1│2│1│││1│3│1│││1│4│1│││1│5│1│││1│6│1││
│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│
├───────┼───────┼───────┼───────┼───────┤
│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│
││1│3│1│││1│4│1│││1│5│1│││1│6│1│││1│7│1││
│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│
├───────┼───────┼───────┼───────┼───────┤
│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│
││1│4│1│││1│5│1│││1│6│1│││1│7│1│││1│8│1││
│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│
├───────┼───────┼───────┼───────┼───────┤
│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│
││1│5│1│││1│6│1│││1│7│1│││1│8│1│││1│9│1││
│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│
└───────┴───────┴───────┴───────┴───────┘
⍕Q qi Hilbert 5 ⍝ rational formatted output
1 1r2 1r3 1r4 1r5
1r2 1r3 1r4 1r5 1r6
1r3 1r4 1r5 1r6 1r7
1r4 1r5 1r6 1r7 1r8
1r5 1r6 1r7 1r8 1r9 The Function qi
The function qi (“rational input”) provides a convenient way to specify rational numbers. It applies independently to items of the argument; each item can be one of several types:
- an ordinary number per APL, e.g. 1234, ¯2.75, 1.234e¯10
- the decimal digits of the numerator, e.g. 1 2 3 4 5 6 7 8 (the denominator is assumed to be 1)
- a character vector of the number per APL, but without the usual limits, e.g. '1.2345e6789'. This is a way to specify exactly a large or tiny rational number.
Examples:
qi ¯2.75
┌──────────┐
│┌───┬─┬──┐│
││1 1│4│¯1││
│└───┴─┴──┘│
└──────────┘
⍕Q qi ¯2.75
¯11r4
qi ⊂1 2 3 4 5 6 7 8
┌─────────────────────┐
│┌───────────────┬─┬─┐│
││1 2 3 4 5 6 7 8│1│1││
│└───────────────┴─┴─┘│
└─────────────────────┘
⍕Q qi ⊂1 2 3 4 5 6 7 8
12345678
qi ○1 ⍝ ⎕ct=1e¯14
┌───────────────────────────────┐
│┌─────────────┬─────────────┬─┐│
││5 4 1 9 3 5 1│1 7 2 5 0 3 3│1││
│└─────────────┴─────────────┴─┘│
└───────────────────────────────┘
18 ⍕ 5419351÷1725033
3.141592653589815___
0 qi ○1 ⍝ ⎕ct=0
┌─────────────────────────────────────┐
│┌─────────────────┬───────────────┬─┐│
││2 4 5 8 5 0 9 2 2│7 8 2 5 6 7 7 9│1││
│└─────────────────┴───────────────┴─┘│
└─────────────────────────────────────┘
18 ⍕ 245850922÷78256779
3.141592653589793___
18 ⍕ ○1
3.141592653589793___
t←qi ⊂'1.2345e6789'
⍴¨¨t ⍝ ≢ decimal digits in the numerator and denominator
┌─────────┐
│┌────┬─┬┐│
││6790│1│││
│└────┴─┴┘│
└─────────┘
¯10 ⍕Q t ⍝ exponential format
1.234500000e6789
The Operator Q
f Q derives a function which works on rational numbers and produces a rational result except as noted otherwise.
Scalar Functions
- monadic and dyadic: + - × ÷ ⌊ ⌈ | !
- dyadic only: * ∨ ∧ < ≤ = ≠ ≥ >
- monadic only: ?
The boolean functions < ≤ = ≠ ≥ > return a boolean result, an ordinary 0 or 1 instead of a rational 0 or 1. Likewise ×⍵ returns an ordinary ¯1, 0, or 1.
Non-Scalar Functions
⊤Q ⍵ applies to rational scalar positive integer ⍵ and produces its binary representation. (⊤ does not have a monadic definition. One possibility is ⊤⍵ ↔ 2⊤⍵, that is, bits or binary representation. (Similarly, ⊥⍵ would be 2⊥⍵, binary value.) This is in fact the definition used in A Formal Description of SYSTEM/360 [Falkoff, Iverson, and Sussenguth 1964, Table 1]. It was actually implemented in APL\360 for a time but was later removed due to space limitations [Hui 1992]. For example, ⊤Q qi 123 ↔ 1 1 1 1 0 1 1.
⍺ ⊤Q ⍵ applies to rational scalar positive integer ⍵ and produces its base-⍺ representation, for an ordinary integer ⍺. For example, 5 ⊤Q qi 1234567 ↔ 3 0 4 0 0 1 2 3 2.
⍕Q ⍵ formats rational array ⍵, with an “r” separating the numerator and the denominator (if not 1).
⍺ ⍕Q ⍵ formats rational array ⍵ in fixed-point (if ⍺≥0) or exponential (if ⍺<0) format. The left argument ⍺ must be a single integer: ⍺≥0 specifies the number of digits after the decimal point; ⍺<0 specifies the number of decimal digits in exponential notation.
⌹Q is analoguous to ⌹ for rational arguments. As usual, a DOMAIN ERROR is signalled if ⍵ is non-singular (but not as usual, if the error is signalled ⍵ really is singular).
Applications
Computing with rational numbers opens new vistas.
Continued Fractions
+∘÷\ 20⍴1
1 2 1.5 1.66667 1.6 1.625 1.61538 1.61905 1.61765 1.61818 1.61798
1.61806 1.61803 1.61804 1.61803 1.61803 1.61803 1.61803 1.61803 1.61803
⍕Q (+Q)∘(÷Q)\ 20 ⍴ qi 1
1 2 3r2 5r3 8r5 13r8 21r13 34r21 55r34 89r55 144r89 233r144 377r233
610r377 987r610 1597r987 2584r1597 4181r2584 6765r4181 10946r6765
How about that: our friends the Fibonacci numbers.
30 ⍕Q ⍪ (+Q)∘(÷Q)\ qi 3 7 15 1 292 1 1 1 2 1 3 1 14 ⍝ OEIS A001203
3.000000000000000000000000000000
3.142857142857142857142857142857
3.141509433962264150943396226415
3.141592920353982300884955752212
3.141592653011902604072261494774
3.141592653921421044708715941593
3.141592653467436705520454785349
3.141592653618936623397500301411
3.141592653581077771204419306582
3.141592653591403978482542414219
3.141592653589389171543687321707
3.141592653589815383241943777307
3.141592653589792659375626945712
Square Root via Newton’s Method
If W is a number whose square root is sought and ⍵ is an approximation of the root, Newton’s method improves the approximation by computing the average of ⍵, and W divided by ⍵. By using rational arithmetic, a sequence of rational approximations obtains.
⍝
sqrt←{ ⍝ apply Newton’s method ⍺ times to find sqrt ⍵
q2←qi 2
W←⍵
f←{q2÷Q⍨⍵+Q W÷Q ⍵}
⍺ {1≥⍺:⍵ ⋄ (⍺-1)∇ ⍵,f ¯1↑⍵} qi 1
}
t←8 sqrt qi 2
t
┌───────┬───────┬───────────┬───────────────┬───────────────────────────┬──
│┌─┬─┬─┐│┌─┬─┬─┐│┌───┬───┬─┐│┌─────┬─────┬─┐│┌───────────┬───────────┬─┐│
││1│1│1│││3│2│1│││1 7│1 2│1│││5 7 7│4 0 8│1│││6 6 5 8 5 7│4 7 0 8 3 2│1││...
│└─┴─┴─┘│└─┴─┴─┘│└───┴───┴─┘│└─────┴─────┴─┘│└───────────┴───────────┴─┘│
└───────┴───────┴───────────┴───────────────┴───────────────────────────┴──
⍕Q t
1 3r2 17r12 577r408 665857r470832 886731088897r627013566048
1572584048032918633353217r1111984844349868137938112 ...
The following code displays the rational approximations to 40 decimal places, and the differences between 2 and the squares of the approximations. The differences illustrate that Newton’s method doubles the number of correct digits with each iteration.
(40 ⍕Q ⍪t), ' ', ¯6 ⍕Q ⍪ (qi 2) -Q t ×Q t
1.0000000000000000000000000000000000000000 1.00000e0
1.5000000000000000000000000000000000000000 ¯2.50000e¯1
1.4166666666666666666666666666666666666667 ¯6.94444e¯3
1.4142156862745098039215686274509803921569 ¯6.00730e¯6
1.4142135623746899106262955788901349101166 ¯4.51095e¯12
1.4142135623730950488016896235025302436150 ¯2.54358e¯24
1.4142135623730950488016887242096980785697 ¯8.08728e¯49
1.4142135623730950488016887242096980785697 ¯8.17550e¯98
Matrix Inverse
Hilbert’s matrix is known to be invertible, but is notorious for being nearly singular.
A← qi Hilbert 8
⍕Q A
1 1r2 1r3 1r4 1r5 1r6 1r7 1r8
1r2 1r3 1r4 1r5 1r6 1r7 1r8 1r9
1r3 1r4 1r5 1r6 1r7 1r8 1r9 1r10
1r4 1r5 1r6 1r7 1r8 1r9 1r10 1r11
1r5 1r6 1r7 1r8 1r9 1r10 1r11 1r12
1r6 1r7 1r8 1r9 1r10 1r11 1r12 1r13
1r7 1r8 1r9 1r10 1r11 1r12 1r13 1r14
1r8 1r9 1r10 1r11 1r12 1r13 1r14 1r15
B← ⌹Q A
⍕Q A (+Q).(×Q) B
1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0
0 0 1 0 0 0 0 0
0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0 0 1 0 0
0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 1
Inverting the 8-by-8 Hilbert matrix as 64-bit floating point numbers:
a←Hilbert 8
b←⌹a
a +.× b
1.00000E0 5.82077E¯11 0.00000E0 ¯7.45058E¯9 1.49012E¯8 ¯1.49012E¯8 ¯1.49012E¯8 0.00000E0
¯2.72848E¯12 1.00000E0 9.31323E¯10 ¯2.60770E¯8 0.00000E0 1.49012E¯8 5.96046E¯8 ¯7.45058E¯9
¯3.63798E¯12 5.82077E¯11 1.00000E0 7.45058E¯9 ¯6.70552E¯8 2.98023E¯8 2.23517E¯8 ¯1.11759E¯8
0.00000E0 8.73115E¯11 1.39698E¯9 1.00000E0 ¯5.21541E¯8 1.49012E¯8 2.23517E¯8 ¯1.11759E¯8
0.00000E0 ¯8.73115E¯11 1.39698E¯9 3.72529E¯9 1.00000E0 2.98023E¯8 2.98023E¯8 ¯1.11759E¯8
1.36424E¯12 1.74623E¯10 ¯2.79397E¯9 1.86265E¯9 ¯1.49012E¯8 1.00000E0 2.98023E¯8 5.58794E¯9
¯2.72848E¯12 1.74623E¯10 ¯1.39698E¯9 1.11759E¯8 ¯3.72529E¯8 3.72529E¯8 1.00000E0 ¯3.72529E¯9
2.27374E¯12 2.91038E¯11 4.65661E¯10 ¯3.72529E¯9 0.00000E0 ¯3.72529E¯8 3.72529E¯8 1.00000E0
The Code
⍝
)copy dfns cfract nats
assert←{⍺←'assertion failure' ⋄ 0∊⍵:⍺ ⎕signal 8 ⋄ shy←0}
Err←{0::⎕io⊃⎕dm ⋄ 0=⎕nc'⍺':⍺⍺ ⍵ ⋄ ⍺ ⍺⍺ ⍵}
qi←{ ⍝ rational input
⍺←⎕CT ⋄ ⎕CT←⍺
qi0←{
d←10|⎕DR ⍵
r←{0=d:0 ⋄ (1E¯15≤|⍵)∧(|⍵)≤1000000000000000}⍵
(d∊1 3)∧1=≢⍴⍵:(⊂,⍵),(⊂,1),⍵≢,0
d∊1 3:(⊂{(1⌈≢⍵)↑⍵}10⊥⍣¯1|⍵),(⊂,1),×⍵
{r∧×d:⍵≡⌊⍵ ⋄ 0}⍵:(⊂10⊥⍣¯1|⍵),(⊂,1),×⍵
r∧5=d:⊃(+Q)∘(÷Q)/∇¨⍺ cfract ⍵
5=d:∇{⎕PP←18 ⋄ ⍕⍵}⍵
assert 0=d: ⍝ from here on ⍵ must be char
assert ⍵∊⎕D,'¯.Ee':
k←⌊/(,⍵)⍳'Ee'
x←(⊂⎕D⍳'.¯'~⍨k↑⍵),(⊂(1⌈k-(,⍵)⍳'.')↑1),¯1*'¯'=⊃⍵
x←x+Q(,¨0 1),0 ⍝ +Q zero to reduce to lowest terms
k=≢⍵:x
e←⍎(1+k)↓⍵
x×Q 1,⍨(0>e)⌽(⊂1,(|e)⍴0),⊂,1
}
⍺ qi0¨⍵
}
Q←{
icarry←{((∧/b)-⍨+/∧\b←0=t)↓t←{1↓+⌿1 0⌽0,0 10⊤⍵}⍣≡0,⍵}
iplus←{icarry(m↑⍺)+m↑⍵⊣m←-(≢⍺)+≢⍵}
iminus←{icarry ⍺-(-≢⍺)↑⍵}
icmp←{0≠s←×(≢⍺)-≢⍵:s ⋄ (s,0)[1⍳⍨0≠s←×⍺-⍵]}
itimes←{
roots←{×\1,1↓(1⌈⍵÷2)⍴¯1*2÷⍵}
cube←{⍵⍴⍨2⍴⍨2⍟≢⍵}
extend←{(2*⌈2⍟¯1+(≢⍺)+≢⍵)↑¨⍺ ⍵}
floop←{(⊣/⍺)∇⍣(×m)⊢(+⌿⍵),[m-0.5]⍺×[⍳m←≢⍴⍺]-⌿⍵}
FFT←{,(cube roots≢⍵)floop cube ⍵}
iFFT←{(≢⍵)÷⍨,(cube+roots≢⍵)floop cube ⍵}
rconvolve←{(¯1+(≢⍺)+≢⍵)↑iFFT⊃×/FFT¨⍺ extend ⍵}
icarry 0,⌊0.5+9○⍺ rconvolve ⍵
}
idiv←{⎕D⍳⎕D[⍺]÷nats ⎕D[⍵]}
iresidue←{⍵ iminus ⍺ itimes ⍵ idiv ⍺}
igcd←{⍵≡,0:⍺ ⋄ ⍵ ∇ ⍵ iresidue ⍺}
irand←{
p←((⊃⍵)+0∨.≠1↓⍵),(¯1+≢⍵)⍴10
{0≥⊃0~⍨⍵-z←?p:∇ ⍵ ⋄ (⍳0∧.=z),(∨\0≠z)⌿z}⍵
}
qnat←{assert(¯1≠⊃⌽⍵)∧(,1)≡1⊃⍵: ⋄ ⍵} ⍝ must be a natural number (≥0)
qcancel←{(⊂,1)≡g←igcd/2↑⍵:⍵ ⋄ ((2↑⍵)idiv¨g),2↓⍵}
qcmp←{
(an ad as)←⍺
(wn wd ws)←⍵
(as≠ws)∨0=as:(as×0≠as)-ws×0=as
as×(an itimes wd)icmp(ad itimes wn)
}
qplus←{
(an ad as)←⍺
(wn wd ws)←⍵
g←ad igcd wd
zd←ad itimes⊢r←wd idiv g
ax←an itimes r
wx←wn itimes ad idiv g
as=ws:qcancel(⊂ax iplus wx),(⊂zd),ws
s←ax icmp wx
0=s:(,¨0 1),0
1=s:qcancel(⊂ax iminus wx),(⊂zd),as
qcancel(⊂wx iminus ax),(⊂zd),ws
}
qnegate←{(2↑⍵),-2↓⍵}
qminus←{⍺ qplus qnegate ⍵}
qtimes←{
(an ad as)←⍺
(wn wd ws)←⍵
qcancel(⊂an itimes wn),(⊂ad itimes wd),as×ws
}
qsignum←{⊃⌽⍵}
qrecip←{assert 0≠⊃⌽⍵: ⋄ ⍵[1 0 2]}
qdiv←{⍺ qtimes qrecip ⍵}
qfloor←{(,1)≡1⊃⍵:⍵ ⋄ (⊂1 iplus⍣(¯1=⊃⌽⍵)⊢(⊃⍵)idiv 1⊃⍵),(⊂,1),¯1↑⍵}
qmin←{⍺⊣⍣(0>⍺ qcmp ⍵)⊢⍵}
qceiling←{qnegate qfloor qnegate ⍵}
qmax←{⍺⊣⍣(0<⍺ qcmp ⍵)⊢⍵}
qabs←{(2↑⍵),|⊃⌽⍵}
qresidue←{0=⊃⌽⍺:⍵ ⋄ ⍵ qminus ⍺ qtimes qfloor ⍵ qdiv ⍺}
qfact←{((⊂,1),1),⍨⊂{(,0)≡⍵:,1 ⋄ ⍵ itimes ∇ ⍵ iminus,1}⊃qnat ⍵}
qbin←{
a←⊃qnat ⍺
w←⊃qnat ⍵
0≤c←a icmp w:(⊂,~c),(⊂,1),~c
d←w iminus a
m←10⊥a⊣⍣(1>a icmp d)⊢d
n←⊃itimes/(0⌷⍵)iminus¨↓⍉10⊥⍣¯1⍳m
d←⊃itimes/↓⍉10⊥⍣¯1⊢1+⍳m
(⊂n idiv d),(⊂,1),1
}
qpow←{
'nonce error'assert(,1)≡1⊃⍵:
0=⊃⌽⍺:(0=⊃⌽qnat ⍵)⊃(⊂⍺),⊂(,¨1 1),1
(2↑⍺)≡,¨1 1:(,¨1 1),¯1*(¯1=⊃⌽⍺)∧2|⊃⌽⊃⍵
'limit error'assert 15≥≢⊃⍵:
a←qrecip⍣(¯1=⊃⌽⍵)⊢⍺
n←10⊥⊃⍵
a{⍺←(,¨1 1),1 ⋄ 0=⍵:⍺ ⋄ (⍺⍺ qtimes⍣(2|⍵)⊢⍺)((⍺⍺ qtimes ⍺⍺)∇∇)⌊⍵÷2}n
}
qgcd←{(qabs ⍺){0=⊃⌽⍵:⍺ ⋄ ⍵ ∇ ⍵ qresidue ⍺}qabs ⍵}
qlcm←{⍺ qtimes ⍵ qdiv ⍺ qgcd ⍵}
qlt←{0>⍺ qcmp ⍵}
qle←{0≥⍺ qcmp ⍵}
qeq←{⍺≡⍵}
qne←{⍺≢⍵}
qge←{0≤⍺ qcmp ⍵}
qgt←{0<⍺ qcmp ⍵}
qrand←{assert((,1)≡1⊃⍵)∧0<⊃⌽⍵: ⋄ (⊂z),(⊂,1),0≢z←irand⊃⍵}
f←(⍺⍺{aa←⍺⍺ ⋄ ⊃⎕CR'aa'}⍵),1+0≠⎕NC'⍺'
e←¯3=≡⍵
⍝ scalar functions
e<'+' 2≡f:⍺ qplus ⍵
e∧'+' 2≡f:⍺ qplus¨⍵
e<'-' 1≡f:qnegate ⍵
e∧'-' 1≡f:qnegate¨⍵
e<'-' 2≡f:⍺ qminus ⍵
e∧'-' 2≡f:⍺ qminus¨⍵
e<'×' 1≡f:qsignum ⍵
e∧'×' 1≡f:qsignum¨⍵
e<'×' 2≡f:⍺ qtimes ⍵
e∧'×' 2≡f:⍺ qtimes¨⍵
e<'÷' 1≡f:qrecip ⍵
e∧'÷' 1≡f:qrecip¨⍵
e<'÷' 2≡f:⍺ qdiv ⍵
e∧'÷' 2≡f:⍺ qdiv¨⍵
e<'⌊' 1≡f:qfloor ⍵
e∧'⌊' 1≡f:qfloor¨⍵
e<'⌊' 2≡f:⍺ qmin ⍵
e∧'⌊' 2≡f:⍺ qmin¨⍵
e<'⌈' 1≡f:qceiling ⍵
e∧'⌈' 1≡f:qceiling¨⍵
e<'⌈' 2≡f:⍺ qmax ⍵
e∧'⌈' 2≡f:⍺ qmax¨⍵
e<'|' 1≡f:qabs ⍵
e∧'|' 1≡f:qabs¨⍵
e<'|' 2≡f:⍺ qresidue ⍵
e∧'|' 2≡f:⍺ qresidue¨⍵
e<'!' 1≡f:qfact ⍵
e∧'!' 1≡f:qfact¨⍵
e<'!' 2≡f:⍺ qbin ⍵
e∧'!' 2≡f:⍺ qbin¨⍵
e<'*' 2≡f:⍺ qpow ⍵
e∧'*' 2≡f:⍺ qpow¨⍵
e<'∨' 2≡f:⍺ qgcd ⍵
e∧'∨' 2≡f:⍺ qgcd¨⍵
e<'∧' 2≡f:⍺ qlcm ⍵
e∧'∧' 2≡f:⍺ qlcm¨⍵
e<'<' 2≡f:⍺ qlt ⍵
e∧'<' 2≡f:⍺ qlt¨⍵
e<'≤' 2≡f:⍺ qle ⍵
e∧'≤' 2≡f:⍺ qle¨⍵
e<'=' 2≡f:⍺ qeq ⍵
e∧'=' 2≡f:⍺ qeq¨⍵
e<'≠' 2≡f:⍺ qne ⍵
e∧'≠' 2≡f:⍺ qne¨⍵
e<'≥' 2≡f:⍺ qge ⍵
e∧'≥' 2≡f:⍺ qge¨⍵
e<'>' 2≡f:⍺ qgt ⍵
e∧'>' 2≡f:⍺ qgt¨⍵
e<'?' 1≡f:qrand ⍵
e∧'?' 1≡f:qrand¨⍵
qbits←{
assert 0=≢⍴⍵:
ihalve←{(1≥⊃⍵)↓(⌊⍵÷2)+¯1↓0,5×2|⍵}
⍬{⍵≡⍬:⍺ ⋄ (⍺,⍨2|⊃⌽⍵)∇(ihalve ⍵)}⊃qnat⊃⍵
}
qbaserep←{
assert 0=≢⍴⍵:
assert(0=≢⍴⍺)∧(1<⍺)∧(⍺<1000000000)∧⍺≡⌊⍺:
2=⍺:qbits ⍵
n←⊃qnat⊃⍵
1=k←10⍟⍺:n
k=⌊k:10⊥⍉((⌈(≢n)÷k),k)⍴(-k×⌈(≢n)÷k)↑n
B←10⊥⍣¯1⊢⍺
⍬{⍵≡,0:⍺ ⋄ (⍺,⍨10⊥⍵ iminus B itimes q)∇ q←⍵ idiv B}⊃qnat⊃⍵
}
qfmt1a←{((¯1=⊃⌽⍵)⍴'¯'),⎕D[0⊃⍵],((,1)≢1⊃⍵)/'r',⎕D[1⊃⍵]}
qfmt2da←{
0=⊃⌽⍵:'0',((×⍺)⍴'.'),⍺⍴'0'
a←⊃qfloor((⊂(⊃⍵),⍺⍴0),(1⌷⍵),1)qplus(,¨1 2),1
('¯'⍴⍨(a≢,0)∧0>⊃⌽⍵),⎕D[(⍳⍺≥≢a),(-⍺)↓a],((×⍺)⍴'.'),⎕D[(-⍺)↑a]
}
qfmt2ea←{
0=⊃⌽⍵:'0.',((¯1-⍺)⍴'0'),'e0'
a←|⍺
(p q)←2↑⍵
k←⌈1+a+(≢q)-≢p
c←(1+a)↑d←((k+≢p)↑p)idiv q
c←⎕D[(a↑c)iplus(-a)↑5≤c[a]]
e←(¯1+≢d)-k
((0>⊃⌽⍵)⍴'¯'),(⊃c),((¯1≠⍺)⍴'.'),(1↓c),'e',⍕e
}
ac←{ ⍝ align columns
s←(×/¯1↓s),¯1↑s←(-2⌈≢⍴⍵)↑1 1,⍴⍵
m←-(~(⊃⌽s)↑1)+⌈⌿s⍴≢¨⍵
0=⎕NC'⍺':m{∊⍺↑¨⍵}⍤1⊢⍵
n←⌈⌿s⍴≢¨⍺
∊⍤1⊢(m↑¨⍤1⊢⍵),¨(n↑¨⍤1⊢⍺)
}
qfmt1←{(326=⎕DR ⍵)⍲∧/,3=≢¨⍵:qfmt1a ⍵ ⋄ ac qfmt1a¨⍵}
qfmt2d←{(326=⎕DR ⍵)⍲∧/,3=≢¨⍵:⍺ qfmt2da ⍵ ⋄ ac ⍺ qfmt2da¨⍵}
qfmt2e←{(326=⎕DR ⍵)⍲∧/,3=≢¨⍵:⍺ qfmt2ea ⍵ ⋄ (e↓¨t)ac(e←t⍳¨'e')↑¨t←⍺ qfmt2ea¨⍵}
qfmt2←{assert(0=≢⍴⍺)∧⍺≡⌊⍺: ⋄ 0≤⍺:⍺ qfmt2d ⍵ ⋄ ⍺ qfmt2e ⍵}
qGauss←{
⍺=n←≢⍵:⍵
'singular'assert n>i←⍺+(0=⊃∘⌽¨⍵[⍺↓⍳n;⍺])⍳0:
w←⍵ ⋄ w[i,⍺;]←w[⍺,i;] ⋄ w[⍺;]←w[⍺;]÷Q¨w[⍺;⍺]
j←(⍳n)~⍺
w[j;]←w[j;]-Q¨w[j;⍺]∘.(×Q)w[⍺;]
(1+⍺)∇ w
}
x←(+Q).(×Q)
qmatinv←{assert 2=≢⍴⍵: ⋄ ≠/⍴⍵:(∇(⍉⍵)x ⍵)x⍉⍵ ⋄ (0,n)↓0 qGauss ⍵,qi∘.=⍨⍳n←≢⍵}
qmatdiv←{assert 2=≢⍴⍵: ⋄ ≠/⍴⍵:(⌹Q(⍉⍵)x ⍵)x(⍉⍵)x ⍺ ⋄ (⌹Q ⍵)x ⍺}
⍝ non-scalar functions
'⊤' 1≡f:qbits ⍵
'⊤' 2≡f:⍺ qbaserep ⍵
'⍕' 1≡f:qfmt1 ⍵
'⍕' 2≡f:⍺ qfmt2 ⍵
'⌹' 1≡f:qmatinv ⍵
'⌹' 2≡f:⍺ qmatdiv ⍵
'not implemented'assert 0:
}
Tests
These tests provide a prima facie case that the code actually works. qtest ignores its argument and returns a result of 1 if all tests passed.
qtest ⍬
1
⍝
qtest←{
assert qtestqi ⍬:
assert qtesta ⍬:
assert qtestsf ⍬:
assert qtestfmt ⍬:
1
}
qtestqi←{
assert(⊂(,¨0 1),0)≡qi 0:
assert(⊂(,¨1 1),1)≡qi 1:
assert(⊂(,¨1 1),¯1)≡qi ¯1:
assert(qi 0)≡qi</1 1:
assert(qi 1)≡qi</0 1:
assert(↓(10⊥⍣¯1¨x),⍤0 1⊢(⊂,1),1)≡qi⊢x←3 7 59:
assert(↓(10⊥⍣¯1¨x×100),⍤0 1⊢(⊂,1),1)≡qi x×100:
assert(↓(10⊥⍣¯1¨x×100000),⍤0 1⊢(⊂,1),1)≡qi x×100000:
assert(⊂(⊂,5),(⊂,2),1)≡qi 2.5:
assert(⊂(⊂,5),(⊂,2),¯1)≡qi ¯2.5:
assert(qi x)≡qi⍤0⊢x←4÷⍨¯50+?3 4⍴100:
assert(⊂(⊂x),(⊂,1),1)≡qi⊂x←1,?40⍴2:
assert(⊂(⊂x),(⊂,1),1)≡qi⊂x←2,?40⍴10:
assert(qi 123)≡qi⊂'123':
assert(qi 12.3 ¯456.78)≡qi'12.3' '¯456.78':
assert(⊂(⊂1,118⍴0),(⊂,1),1)≡qi⊂'1e118':
assert(⊂(⊂,1),(⊂1,118⍴0),¯1)≡qi⊂'¯1e¯118':
assert(⊂(⊂2 4 6 9),(⊂,2 0),¯1)≡qi⊂'¯1.2345e2':
assert(qi 0 0)≡qi'0e15' '0e¯20':
1
}
qtesta←{
assert ¯1≡×Q qi ¯12345:
assert 0≡×Q qi 0:
assert 1≡×Q qi 59.3:
assert(⌊Q qi 56.4)≡qi 56:
assert(⌊Q qi ¯56.4)≡qi ¯57:
assert(⌊Q qi 12345)≡qi 12345:
assert(⌊Q qi ¯12345)≡qi ¯12345:
assert(⌊Q qi 0)≡qi 0:
assert(⌈Q qi 56.4)≡qi 57:
assert(⌈Q qi ¯56.4)≡qi ¯56:
assert(⌈Q qi 12345)≡qi 12345:
assert(⌈Q qi ¯12345)≡qi ¯12345:
assert(⌈Q qi 0)≡qi 0:
assert(|Q qi 4.5)≡qi 4.5:
assert(|Q qi ¯4.5)≡qi 4.5:
assert(|Q qi 0)≡qi 0:
assert((qi x)|Q qi y)≡qi(x←2/1.3 ¯1.3)|y←4⍴2.8 ¯2.8:
assert((qi 0)|Q qi x)≡qi⊢x←¯4.5 0 12.4:
assert(!Q qi 0)≡qi 1:
assert(!Q qi 1)≡qi 1:
assert(!Q qi 12)≡qi!12:
assert(⍕Q!Q qi 39)≡'20397882081197443358640281739902897356800000000':
assert((qi 10)!Q qi 19)≡qi 10!19:
x←'1647278652451762678788128833110870712983038446517480945400'
assert(⍕Q(qi 120)!Q qi 200)≡x:
assert((qi 4÷3)*Q qi⍳10)≡qi(4÷3)*⍳10:
assert((qi 4÷3)*Q qi-⍳10)≡qi(3÷4)*⍳10:
assert((qi 0)*Q qi 5)≡qi 0:
assert((qi ¯1)*Q qi x)≡qi ¯1*x←?4 5⍴100:
assert((qi 1)*Q qi x)≡qi(⍴x)⍴1:
x←qi⊢y←¯3.47 ¯2 0 2 3.47
assert(∘.(<Q)⍨x)≡∘.<⍨y:
assert(∘.(≤Q)⍨x)≡∘.≤⍨y:
assert(∘.(=Q)⍨x)≡∘.=⍨y:
assert(∘.(≠Q)⍨x)≡∘.≠⍨y:
assert(∘.(≥Q)⍨x)≡∘.≥⍨y:
assert(∘.(>Q)⍨x)≡∘.>⍨y:
assert 1=x>Q?Q 100⍴x←qi 13:
assert(⊤Q qi⊂x)≡2⊥⍣¯1⊢10⊥x←?10⍴10:
assert(⊤Q(qi 2)*Q qi 99)≡100↑1:
assert(⊤Q qi 0)≡,0:
assert(⊤Q qi 1)≡,1:
x←(+Q).(×Q)
assert(⍕∘.=⍨⍳5)≡⍕Q(⌹Q a)x⊢a←qi?5 5⍴50:
assert(b⌹Q a)≡(⌹Q a)x⊢b←qi?5⍴50:
assert(a x b⌹Q a)≡b:
assert(⌹Q a)≡(⌹Q(⍉a)x a)x⍉a←qi?10 5⍴50:
assert(b⌹Q a)≡(⌹Q(⍉a)x a)x(⍉a)x⊢b←qi?10⍴50:
assert qtestfmt ⍬:
assert'assertion failure'≡(qi 1 2)÷Q Err qi 0:
assert'assertion failure'≡!Q Err qi ¯1:
assert'assertion failure'≡!Q Err qi 1.2:
1
}
qtestsf←{ ⍝ test scalar functions
do←{
f←⍎⍺⍺
s←⍺⍺∊'!*'
p←qi⊢x←1+?7 3[s]⍴10 3[s]
q←qi⊢y←1+?7 3[s]⍴10 3[s]
assert(qi⍣⍵⊢x f y)≡p f Q q:
assert(qi⍣⍵⊢(0⌷x)f y)≡(0⌷p)f Q q:
assert(qi⍣⍵⊢x f 0⌷y)≡p f Q 0⌷q:
assert(qi⍣⍵⊢x∘.f y)≡p∘.(f Q)q:
~⍵:1
assert(qi⍣⍵⊢f/x)≡f Q/p:
assert(qi⍣⍵⊢f\x)≡f Q\p:
1
}
assert{⍵ do 0}¨'<≤=≠≥>':
assert{⍵ do 1}¨'+-×÷⌊⌈|!*∨∧':
x←4÷⍨¯20+?4 5⍴50
y←4÷⍨¯20+?5 7⍴60
assert((qi x)(+Q).(×Q)qi y)≡qi x+.×y:
1
}
qtestfmt←{
assert(⍕Q qi(⍳5)÷2)≡'0 1r2 1 3r2 2':
assert(⍕Q qi-÷1+⍳5)≡'¯1 ¯1r2 ¯1r3 ¯1r4 ¯1r5':
x←3 2⍴qi'123' '0.6789' '¯1.0625' '¯0.5' '2' '¯0.125e¯1'
assert(⍕Q x)≡∊⍤1⊢¯6 ¯11↑¨⍤1⍕Q¨x:
f←{('E',⎕AV)[('e',⎕AV)⍳¯6⍕Q ⍵]}
g←{1↓¯6⍕⍵}
assert{assert(f qi ⍵)≡g ⍵: ⋄ 1}¨⍳200:
assert{assert(f qi ⍵)≡g ⍵: ⋄ 1}¨÷1+⍳200:
assert(x←f qi 3.14159)≡y←g 3.14159:
assert(¯6⍕Q qi⊂'1.2345e6789')≡'1.23450e6789':
assert(¯6⍕Q qi⊂'¯1.2345e6789')≡'¯1.23450e6789':
assert(¯6⍕Q qi⊂'1.2345e¯6789')≡'1.23450e¯6789':
assert(¯6⍕Q qi⊂'¯1.2345e¯6789')≡'¯1.23450e¯6789':
assert(¯4⍕Q qi⊂'1.2345e6789')≡'1.235e6789':
assert(¯4⍕Q qi⊂'¯1.2345e6789')≡'¯1.235e6789':
assert(¯4⍕Q qi⊂'1.2345e¯6789')≡'1.235e¯6789':
assert(¯4⍕Q qi⊂'¯1.2345e¯6789')≡'¯1.235e¯6789':
assert(¯1⍕Q qi 1.2)≡'1e0':
assert(¯1⍕Q qi 16)≡'2e1':
assert(¯1⍕Q qi ¯0.044)≡'¯4e¯2':
assert(¯1⍕Q qi ¯0.065)≡'¯7e¯2':
assert(1⍕Q qi 0 12.34 12.36 ¯12.34 ¯12.36)≡'0.0 12.3 12.4 ¯12.3 ¯12.4':
assert(0⍕Q qi 0 2.4 2.6 ¯2.4 ¯2.6)≡'0 2 3 ¯2 ¯3':
assert(2⍕Q qi ¯0.006 ¯0.004 0.004 0.006)≡'¯0.01 0.00 0.00 0.01':
assert(¯3⍕Q qi 10*x)≡1↓∊' 1.00e'∘,¨⍕¨x←¯5+⍳11:
x←qi'12e¯4' '666e889'(○¯1)
y←qi'¯12e1234' 0(*1)
assert(¯6⍕Q⍪x)≡↑'/'(≠⊆⊢)'/ 1.20000e¯3/ 6.66000e891/¯3.14159e0':
assert(¯6⍕Q⍪y)≡↑'/'(≠⊆⊢)'/¯1.20000e1235/ 0.00000e0/ 2.71828e0':
assert(¯6⍕Q x,⍤0⊢y)≡(¯6⍕Q⍪x),' ',¯6⍕Q⍪y:
1
}
⍕Q has a separate set of tests because implementing the functionality was ... interesting.
created:_2020-07-07
updated: 2020-08-24
- Roger|Dyalog
- Posts: 238
- Joined: Thu Jul 28, 2011 10:53 am
1 post
• Page 1 of 1
Who is online
Users browsing this forum: No registered users and 1 guest
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group