1 ;;; Rational numbers
  2 
  3 import math ; abs, divmod, gcd, pow, sign
  4 import util ; die!, neq
  5 
  6 ; maximum value of either component of a rational number before behavior is
  7 ; undefined; 2^31 by default to give interpreters without bignums a chance,
  8 ; but customizable from userland.
  9 $RAT = push 2,31 :pow
 10 
 11 ; encodes a numerator N and a denominator D as a rational number R with the
 12 ; following structure: ((N × $RAT + D) << 1) + sign (0 for negative, else 1).
 13 ; This representation is nice because, with care, it makes it so that we never
 14 ; have to involve the heap (besides divmod) just to do fairly basic arithmetic.
 15 ; [N D] => [R]
 16 ;
 17 ; [22 7]  => [R(22,7)]
 18 ; [4 8]   => [R(4,8)]  ; no implicit simplification TODO: make it configurable?
 19 ; [5 1]   => [R(5,1)]  ; no conversion to integer (duh)
 20 ; [-3 4]  => [R(-3,4)]
 21 ; [3 -4]  => [R(-3,4)] ; sign always held in the numerator
 22 ; [-1 -3] => [R(1,3)]
 23 to_r: dup jz _to_r_dbz
 24   dup :sign copy 2 :sign mul push -1 :neq
 25   copy 2 :abs copy 2 :abs
 26   $RAT copy 2 mul add push 2 mul
 27   copy 2 add slide 4 ret
 28 _to_r_dbz: push "0 in denominator!" :die!
 29 
 30 ; returns the numerator N of the rational number R, which may be negative
 31 ; [R] => [N]
 32 ;
 33 ; [R(22,7)] => [22]
 34 ; [R(-3,4)] => [-3]
 35 ; [R(3,-4)] => [-3]
 36 ratnum: push 2 :divmod push 2 mul push 1 sub swap $RAT div mul ret
 37 
 38 ; returns the denominator D of the rational number R; always positive
 39 ; [R] => [D]
 40 ;
 41 ; [R(22,7)] => [7]
 42 ; [R(-3,4)] => [4]
 43 ; [R(3,-4)] => [4]
 44 ratden: push 2 div $RAT mod ret
 45 
 46 ; decomposes the rational number R into its numerator N and denominator D
 47 ; [R] => [N D]
 48 ;
 49 ; [R(22,7)] => [22 7]
 50 ; [R(-3,4)] => [-3 4]
 51 ; [R(3,-4)] => [-3 4]
 52 ; [R(4,8)]  => [4 8]  ; no implicit simplification
 53 from_r: dup :ratnum swap :ratden ret
 54 
 55 ; fully simplifies the rational number R by dividing both of its components
 56 ; by their greatest common divisor
 57 ; [R] => [R, simplified]
 58 ;
 59 ; [R(4,8)] => [R(1,2)]
 60 ; [R(-3,12)] => [R(-1,4)]
 61 ; [R(17,-51)] => [R(-1,3)]
 62 ratsimp:
 63   push 2 :divmod swap $RAT :divmod dup copy 2 :gcd
 64   swap copy 1 div copy 2 copy 2 div swap :to_r
 65   slide 2 push 2 div push 2 mul add ret
 66 
 67 ; helper prologue shared by all but ratmul
 68 _rathead:
 69   copy 1 :ratnum copy 1 :ratden mul
 70   copy 2 :ratden copy 2 :ratnum mul ret
 71 
 72 ; helper epilogue shared by all but ratdiv
 73 _rattail:
 74   copy 2 :ratden copy 2 :ratden mul
 75   :to_r :ratsimp slide 2 ret
 76 
 77 ; returns the simplified sum of the rational numbers Ra and Rb
 78 ; [Ra Rb] => [Ra + Rb]
 79 ;
 80 ; [R(-9,-14) R(-3,19)]  => [R(129,266)]
 81 ; [R(17,30)  R(18,10)]  => [R(71,30)]
 82 ; [R(-27,14) R(15,-23)] => [R(-831,322)]
 83 ; [R(3,-9)   R(8,3)]    => [R(7,3)]
 84 ; [R(-5,27)  R(-2,-27)] => [R(-1,9)]
 85 ; [R(-27,-8) R(-15,22)] => [R(237,88)]
 86 ; [R(-9,-29) R(-27,3)]  => [R(-252,29)]
 87 ; [R(2,-21)  R(4,6)]    => [R(4,7)]
 88 ratadd: :_rathead add :_rattail ret
 89 
 90 ; returns the simplified difference of the rational numbers Ra and Rb
 91 ; [Ra Rb] => [Ra - Rb]
 92 ;
 93 ; [R(21,25)  R(-28,27)]  => [R(1267,675)]
 94 ; [R(14,7)   R(13,6)]    => [R(-1,6)]
 95 ; [R(-24,-9) R(-5,-21)]  => [R(17,7)]
 96 ; [R(-27,-2) R(-2,26)]   => [R(353,26)]
 97 ; [R(-27,3)  R(2,-22)]   => [R(-98,11)]
 98 ; [R(-4,23)  R(-9,13)]   => [R(155,299)]
 99 ; [R(-14,19) R(-18,-11)] => [R(-496,209)]
100 ; [R(-29,21) R(-15,-16)] => [R(-779,336)]
101 ratsub: :_rathead sub :_rattail ret
102 
103 ; returns the simplified product of the rational numbers Ra and Rb
104 ; [Ra Rb] => [Ra × Rb]
105 ;
106 ; [R(-24,26) R(-1,-30)]  => [R(-2,65)]
107 ; [R(19,4)   R(27,2)]    => [R(513,8)]
108 ; [R(25,27)  R(4,-11)]   => [R(-100,297)]
109 ; [R(1,18)   R(4,8)]     => [R(1,36)]
110 ; [R(1,27)   R(-8,29)]   => [R(-8,783)]
111 ; [R(25,-13) R(-6,24)]   => [R(25,52)]
112 ; [R(6,-13)  R(9,23)]    => [R(-54,299)]
113 ; [R(11,8)   R(-19,-19)] => [R(11,8)]
114 ratmul: copy 1 :ratnum copy 1 :ratnum mul :_rattail ret
115 
116 ; returns the simplified quotient of the rational numbers Ra and Rb
117 ; [Ra Rb] => [Ra / Rb]
118 ;
119 ; [R(-30,-22) R(15,12)]   => [R(12,11)]
120 ; [R(13,28)   R(-15,29)]  => [R(-377,420)]
121 ; [R(7,-30)   R(-22,-12)] => [R(-7,55)]
122 ; [R(15,4)    R(8,-8)]    => [R(-15,4)]
123 ; [R(-23,28)  R(-16,-15)] => [R(-345,448)]
124 ; [R(-18,12)  R(6,18)]    => [R(-9,2)]
125 ; [R(29,-2)   R(11,-21)]  => [R(609,22)]
126 ; [R(-23,25)  R(25,-3)]   => [R(69,625)]
127 ratdiv: :_rathead :to_r :ratsimp slide 2 ret
128 
129 ; returns the simplified modulus of the rational numbers Ra and Rb
130 ; [Ra Rb] => [Ra % Rb]
131 ;
132 ; [R(-15,-3) R(-16,-10)] => [R(1,5)]
133 ; [R(4,2)    R(21,21)]   => [R(0,1)]
134 ; [R(24,10)  R(-18,-3)]  => [R(12,5)]
135 ; [R(3,-7)   R(-2,16)]   => [R(-3,56)]
136 ; [R(4,28)   R(-29,7)]   => [R(-4,1)]
137 ; [R(7,-27)  R(10,23)]   => [R(109,621)]
138 ; [R(28,-3)  R(30,-12)]  => [R(-11,6)]
139 ; [R(-29,21) R(19,-23)]  => [R(-268,483)]
140 ratmod: :_rathead mod :_rattail ret
141 
142 ; returns the sign of the rational number R
143 ; [R] => [-1 | 0 | 1]
144 ;
145 ; [R(4,3)] => [1]
146 ; [R(4,-3)] => [-1]
147 ; [R(0,10)] => [0]
148 ; [R(-3,4)] => [-1]
149 ; [R(-3,-4)] => [1]
150 ratsign:
151   dup $RAT push 2 mul push 2 add sub jn _ratsign_zero
152   push 2 mod push 2 mul push 1 sub ret
153 _ratsign_zero: pop push 0 ret
154 
155 ; returns the simplified inverse of the rational number R
156 ; [R] => [1/R]
157 ;
158 ; [R(4,8)] => [R(2,1)]
159 ; [R(3,-5)] => [R(-5,3)]
160 ; [R(-22,-7)] => [R(7,22)]
161 ratinv:
162   push 2 :divmod push 2 mul $--
163   swap $RAT :divmod copy 2 mul
164   swap :to_r :ratsimp slide 1 ret
165 
166 ; returns the rational R raised to the power E
167 ; [R E] => [R']
168 ;
169 ; [R(4,2) 5] => [R(32,1)]
170 ; [R(-3,7) 0] => [R(1,1)]
171 ; [R(-3,14) 2] => [R(9,196)]
172 ; [R(5,1) -1] => [R(1,5)]
173 ; [R(2,5) -2] => [R(25,4)]
174 ; [R(-3,14) -3] => [R(-2744,27)]
175 ratpow: dup jn _ratpow_neg
176   swap push 2 :divmod push 2 mul $--
177   swap $RAT :divmod copy 2 mul
178   copy 3 :pow swap copy 3 :pow
179   swap :to_r :ratsimp slide 2 ret
180 _ratpow_neg: push -1 mul swap :ratinv swap :ratpow ret
181 
182 ; converts the rational number R to a "pseudo-float" with a whole (W) part
183 ; and a fractional (F) part composed of P digits after the decimal point
184 ; [R P] => [W F]
185 ;
186 ; [R(22,7) 2] => [3 14]
187 ; [R(355,113) 6] => [3 141592]
188 ; [R(8675,309) 10] => [28 744336569] TODO: leading 0 is lost (bug)
189 ; [R(2,4) 3] => [0 500]
190 to_f: ^-2 :from_r :divmod push 0 swap
191 _to_f_loop:
192   @-1 :divmod  swap
193   copy 2 push 10 mul add
194   swap push 10 mul
195   push 2 :dig
196   push -2 dup :dec load :neg? jz _to_f_loop pop ret
197 
198 ; converts the rational number R to a string S of the form
199 ; "<simplified numerator>/<simplified denominator>"
200 ; [R] => [S]
201 ;
202 ; [R(22,7)] => ["22/7"]
203 ; [R(2,4)] => ["1/2"]
204 ; [R(99,9)] => ["11/1"]
205 ratstr: :ratsimp :from_r push 2 map (:to_s) push '/' :strjoinc ret