1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
| \ part 2
big_digit_pointer multiplicand big_digit_pointer multiplier
big_digit_pointer product
: big_product ( a1 a2 -- a3 ) \ a3 has the result of multiplying
\ the absolute value of the n digit operand at a1 by the absolute
\ value of the m digit operand at a2.
to_pointer multiplier
to_pointer multiplicand \ store operand addresses
here dup to_pointer product \ address of result
0 multiplier @ abs
0 multiplicand @ abs 2dup
+ , \ store product size
\ allot and clear the first
dup 0 do 0 , loop \ n digits of the product
over cells allot \ allot remaining digits of product
over 1+ 1 \ for each multiplier digit,
\ starting with the first
do
0 \ initial carry
over 1+ 1 \ for each multiplicand digit,
\ starting with the first
do
i multiplicand @ \ mulitplicand digit times
j multiplier @ \ multiplier digit
m*
i j 1- + >r \ current product digit index
r@ product @ m+ \ add previous product result
rot m+ \ add carry
d>carry \ split into digit & carry
r> product ! \ store product digit
loop
over i + product ! \ store carry
loop
2drop
dup clip ; \ if there is a high-order zero,
\ remove it
big_digit_pointer dividend big_digit_pointer divisor
big_digit_pointer quotient variable normalizer
: divisor(n) ( -- n) \ n is the high digit of the divisor
0 divisor @ abs divisor @ ;
: divisor(n-1) ( -- n) \ n is the next-to-high-order digit of the divisor
0 divisor @ abs 1- divisor @ ;
: normalize ( -- ) \ multiply dividend and divisor by a factor that
\ will guarantee that the leading "digit" of the divisor will be
\ > bigbase/2
bigbd \ big number base as double number
divisor(n) \ high order digit of divisor
1+ um/mod normalizer ! \ normalizing factor (base/(vn+1))
drop \ discard remainder
here \ this will be the address of the
\ normalized dividend
0 dividend big>here \ copy dividend to end of data space
to_pointer dividend \ new dividend address
normalizer @ 1 >
if 0 dividend
normalizer @ big*s \ normalize the dividend.
then
0 , \ append high order zero to dividend
0 dividend dup @ 0< \ negative dividend?
if -1 else 1 then
swap +! \ up the dividend digit count
here \ address of the normalized divisor
0 divisor big>here \ copy divisor to end of data space
dup to_pointer divisor
normalizer @ big*s ; \ normalize the divisor
: big. ( a --) \ "big dot" display the big number at a
here >r r@
swap big>here \ copy for nondestructive write
dup @ \ sign of the number
swap bigstring
space
r> here - allot ; \ recover space used by big>here
: big.digits ( a --) \ "big dot digits" print the digits of the
\ big number at a
dup cell+ swap dup @ abs cells + do i ? -1 cells +loop ;
: trial ( n1 -- n2) \ n2 is trial quotient digit n1
\ cr ." trial "
0 divisor @ abs + >r
r@ dividend @ bigbase um* \ u(j)*b
r@ 1- dividend @ m+ \ [u(j)*b+u(j-1)]
divisor(n) \ v(1), high digit of divisor
r@ dividend @ = \ equal to uj?
\ data stack: low[u(j)*b+u(j-1)]
\ high[u(j)*b+u(j-1)]
\ flag
if 2drop
r@ 1- dividend @
0 divisor(n) m+ \ rhat = u(j-1) + v(1)
biggest
swap
if r> drop exit then \ we have the right q
else divisor(n) um/mod \ rhat qhat
then \ ( rhat qhat) (j)
begin \ test trial quotient
2dup divisor(n-1) um* \ v(n-1)*qhat
rot bigbase um* \ rhat*b
r@ 2 - dividend @ \ u(j-2)
m+
2swap d<
while \ ( rhat qhat) (j)
1- \ decrease trial quotient
swap divisor(n) + \ adjust remainder
swap
repeat
r> drop \ clear return stack
nip ; \ drop trial remainder
: div_subtract ( quotient j -- quotient flag)
\ subtract (vn..v1)q from (u(j+n)..u(j))
\ flag is true if the result is negative
0 \ borrow
0 divisor @ abs 1+ 1
do \ ( quotient j borrow)
>r 2dup r>
rot i divisor @ m*
d>carry \ convert from double number to
\ big digits
\ ( quotient j j borrow carry digit)
rot + \ add the previous borrow to the digit
overflow?
rot dividend @ \ uj
rot - \ new uj
begin
overflow? over 0<
while swap repeat
>r \ park new borrow
over dividend ! \ store new uj
1+ \ bump j
r>
loop
over dividend @
swap - \ subtract the last borrow from the
\ next digit of u
dup
rot dividend ! \ put the result in the digit of u
0<> ; \ test for overflow
: addback ( j --) \ add (vn..v1) to (u(j+n)..u(j))
0 \ carry
0 divisor @ abs 1+ 1
do \ j carry
over dup dividend @ \ j carry j u(j)
i divisor @
+ rot + \ j j (v(i)+u(j)+carry)
carry
rot dividend !
swap 1+ swap \ increment j
loop
dup
if \ deal with the final carry (i'm not sure
\ this is strictly necessary (if you can
\ prove it one way or the other, i would be
\ interested in seeing it) but it is neater)
swap dividend +!
else 2drop
then ;
: |divide| ( a1 a2 -- a3) \ a3 contains the result of dividing
\ the absolute value of the big number at a1 by the absolute value
\ of the big number at a2. the numbers must be unequal and the
\ divisor must have at least two "digits".
to_pointer divisor
to_pointer dividend
normalize
here dup \ address of quotient
to_pointer quotient
1 \ limit for do - stop after digit 1
0 dividend @ abs \ number of digits in normalized dividend
0 divisor @ abs \ number of digits in divisor
- 1 max dup , \ number of digits in quotient
dup cells allot \ space for quotient
do
i trial \ trial quotient digit
i div_subtract
if 1- i addback then \ ( qi)
i quotient ! \ store qi
-1 +loop
dup clip ;
: divide ( a1 a2 -- a3) \ a3 contains the result of dividing
\ the absolute value of the big number at a1 by the absolute value
\ of the big number at a2.
2dup |big|< \ is the number at a1 < num at a2?
if 2drop here 1 , 0 , \ answer is 0
else 2dup |big|= \ are the numbers equal?
if 2drop here 1 , 1 , \ answer is 1
else dup @ abs 1 = \ single "digit" divisor?
if cell+ @ \ divisor on stack
here rot big>here \ dividend to here
dup rot big/mods
drop \ drop remainder
dup @ abs \ absolute value for sign of quotient
over !
else |divide|
then
then
then ;
\ finally! the words for the user.
: big ( <cccc> -- a ) \ a is the address of the big number
\ created from characters cccc in the input stream.
bl word count \ ( a u )
>r big_string r@ cmove \ move characters from input stream to buffer
big_string r> make_big_number ;
big_digit_pointer op1 big_digit_pointer op2
: big+ ( a1 a2 -- a3) \ a3 has the result of algebraically
\ adding the operand at a1 to the operand at a2.
here >r
2dup to_pointer op2
to_pointer op1
0 op1 @ 0 op2 @ xor
0< if difference \ operands are of opposite sign
0 op1 0 op2 |big|<
if 0 op2 @ \ result has the sign of operand 2
else 0 op1 0 op2 |big|=
if 1 \ result is zero, plus sign
else 0 op1 @ \ result has the sign of operand 1
then
then
else sum \ operands have same sign
0 op1 @
then
over @ \ size of result
swap 0< if negate then \ add the sign
over !
r@ reposition r> ;
: big- ( a1 a2 - a3) \ a3 has the result of algebraically
\ subtracting the operand at a2 from the operand at a1.
here >r big>here \ copy second operand
r@ bignegate \ switch its sign
r@ big+ \ add
r@ reposition r> ;
: big* ( a1 a2 -- a3) \ a3 is the address of the result of
\ multiplying the operand at a1 by the operand at a2
2dup big_product adjust_sign ;
: big/ ( a1 a2 -- a3) \ a3 contains the floored quotient of
\ the big number at a1 dvided by the big number at a2.
\ a3 is the value of here before the operation.
here >r
2dup divide
( adjust_sign)
rot @ rot @ xor 0< \ do we need an adjustment?
if dup 1 big+s dup bignegate then
r@ reposition r> ;
\ big 288,265,561,597,526,014 big 17,593,259,786,239 big/ should leave a
\ result of 16384. this tests the rare "trial divisor off by two" division
\ branch on a 16 bit system. see regener for more on this
: bigmod ( a1 a2 -- a3) \ a3 is the remainder after dividing
\ the big number at a1 by the big number at a2. a3 is the value
\ returned by here before the operation.
here >r
2dup big/ \ (a1 a2 qoutient-a)
big* big-
r@ reposition r> ;
: big/mod ( a1 a2 -- a3 a4) \ a3 is the remainder and a4
\ is the quotient after dividing the big number at a1 by the big
\ number at a2.
2dup big/
dup >r
big* big-
r> ;
base ! \ End of file; restore BASE
|